text:
# Proton Spin Relaxation Induced by Quantum Tunneling in Fe8 Molecular Nanomagnet

## Abstract

The spin-lattice relaxation rate $`T_1^1`$ and NMR spectra of 1H in single crystal molecular magnets of Fe8 have been measured down to 15 mK. The relaxation rate $`T_1^1`$ shows a strong temperature dependence down to 400 mK. The relaxation is well explained in terms of the thermal transition of the iron state between the discreet energy levels of the total spin $`S=10`$. The relaxation time $`T_1`$ becomes temperature independent below 300 mK and is longer than 100 s. In this temperature region stepwise recovery of the 1H-NMR signal after saturation was observed depending on the return field of the sweep field. This phenomenon is attributed to resonant quantum tunneling at the fields where levels cross and is discussed in terms of the Landau-Zener transition.

Recently, molecular nanomagnets have attracted much attention to study quantum mechanical phenomena in the macroscopic system owing to their identical size, the well-defined structure and a well-characterized energy structure.  The molecular magnet \[(C6H15N3)6Fe8O2(OH)12\]Br7(H2O)Br$``$8H2O, abbreviated Fe8, is a representative compound in which quantum tunneling of magnetization (QTM) has been observed as the temperature-independent recovery of magnetization below 400 mK.

The molecular magnet Fe8 consists of eight Fe3+ ions with spin $`s=5/2`$ in each molecule. The magnetic interactions between the spins in the molecule are antiferromagnetic and their magnitudes are $`22170`$ K,  while the magnetic interactions between the molecules are negligibly small. The magnetic properties of this compound at low temperatures have been described by a total spin of $`S=10`$ for each molecule, in which six spins are parallel to each other and the remaining two spins are antiparallel to the six spins.  The spin Hamiltonian in the field $`𝑯`$ is expressed by

$`=DS_z^2+E(S_x^2S_y^2)+g\mu _\mathrm{B}𝑺𝑯,`$ (1)
where $`D`$ and $`E`$ are the easy axis and the in-plane anisotropy, respectively.  We have determined the values of $`D=0.276`$ K and $`E=0.035`$ K by magnetization measurement on a single crystal.  When there is no magnetic field, the anisotropy stabilizes the degenerate spin states of $`m=\pm 10`$ at low temperatures. These states correspond to opposite directions of the magnetization in the classical sense. The energy barrier caused by the anisotropy for reversal of the magnetization between $`m=\pm 10`$ is reported to be 25 K.

Quantum tunneling of magnetization in Fe8 has been observed only by magnetization measurements.  In order to study the spin dynamics in Fe8 from a microscopic point of view we performed NMR experiments and found a stepwise recovery of the 1H spin echo signal due to resonant quantum tunneling at the level crossing fields.

We synthesized single crystals of Fe8, following the method reported previously.  The sizes of the crystals used for the experiments were about 4 $`\times `$ 2 $`\times `$ 1 mm3. The NMR spectrum and the spin-lattice relaxation rate $`T_1^1`$ of 1H in Fe8 were measured by the coherent pulsed NMR method in external magnetic fields up to 5.4 T. The spectra were obtained by measuring the spin-echo intensity with sweeping the field at a fixed frequency. The relaxation rate was obtained by measuring the recovery of the spin echo intensity as a function of the time after saturation of the 1H spins using comb pulses. The experimental temperature was lowered down to 15 mK using a 3He-4He dilution refrigerator. The samples were sealed with 4He gas in a cell made of PCTFE plastic(polychlorotrifluoroethylene). The cell was set in the mixing chamber of the dilution refrigerator.

The NMR spectrum at room temperature is narrow with the width smaller than 50 Oe, while the spectrum becomes broader and shows structure at lower temperatures. The structure is caused by freezing of the iron magnetization and the existence of many 1H sites.

Figure 1 shows the temperature dependence of the relaxation rate $`T_1^1`$. The rate $`T_1^1`$ decreases steeply over six decades with lowering the temperature from 10 K to 400 mK, and decreases with increasing field. Below 300 mK $`T_1^1`$ is almost temperature independent and has a strong site dependence.

Figure 2 shows NMR spectra of the single crystal at 150 mK. The field was applied so as to make an angle $`\theta `$ of $`50^{}`$ from the easy axis and to be within the $`ab`$-plane. In this experiment we first saturated the 1H spin system by the comb pulses at a fixed frequency, sweeping the field up and down between 0.1 T and 1.5 T for more than three times, until there was no observation of the echo signals. Then the field was decreased to a certain field $`H_r`$ with a constant sweep rate $`dH/dt`$. Immediately after arriving at $`H_r`$, the field was increased with the same sweep rate and the spectrum was taken at a fixed repetition time. No signals were expected to be observed, because the returning duration after the saturation is planned to be enough shorter than the relaxation time that is longer than 100 s. Indeed no signals except from 19F in the material of the sample cell and 3He in the mixing chamber were observed when the field was returned at 0.05 T, as shown in the lowest spectrum in Fig. 2. However, signals were observed when the return fields were negative. The intensities of the signals were increased with decreasing return field.

Figures 3(a) and (b) show the return field dependence of the echo intensity. Figure 3(a) was obtained in the case when the field was applied with $`\theta =50^{}`$ within the $`ab`$-plane. Figure 3(b) is the case when the field was parallel to the easy axis $`(\theta =0^{})`$. The echo intensities were picked up at fields of 0.45 T and 0.60 T for Fig. 3(a) and (b), respectively. The intensity increases sharply with steps at 0, $`0.31`$,$`0.63`$ and $`1.00`$ T for Fig. 3(a), while for Fig. 3(b) the steps are at +0.02 and $`0.20`$ T. The increase in the intensity with the steps cannot be explained by recovery due to the spin-lattice relaxation time shown in Fig. 1. Moreover the intensity depends on the sweep rate of the field. Figure 4 shows the sweep rate dependence of the intensity at a field of 0.60 T when the field was applied parallel to the easy axis and the return fields were $`0.1`$ T and $`0.3`$ T. The intensity is large when the sweep rate is slow.

First we analyze the relaxation rate. The temperature dependent relaxation rates of a molecular magnet Mn12 were first discussed by Lascialfari et al.  In general $`T_1^1`$ is given by the Fourier transform of the correlation function for fluctuating transverse local fields $`h_\pm (t)`$ at nuclear sites, and is expressed as

$`{\displaystyle \frac{1}{T_1}}`$ $`=`$ $`{\displaystyle \frac{\gamma _N^2}{2}}{\displaystyle h_\pm (t)h_{}(0)e^{i\omega _Lt}𝑑t},`$ (2)
where $`\gamma _N`$ is the nuclear gyromagnetic ratio and $`\omega _L`$ is the Larmor frequency. When the correlation is assumed to be described by an exponential function with a life time $`\tau _m`$ in the $`m`$ state, the correlation can be expressed as

$`h_\pm (t)h_{}(0)={\displaystyle \underset{m=10}{\overset{+10}{}}}\mathrm{\Delta }h_\pm ^2e^{t/\tau _m}{\displaystyle \frac{e^{E_m/k_\mathrm{B}T}}{Z}},`$ (3)
where $`E_m`$ is the energy of the eigenstate $`m`$ and $`Z`$ is the partition function. The lifetime $`\tau _m`$ for the $`m`$ state is governed by the spin-phonon interaction $`_{sp}`$ and is expressed by the transition probabilities $`p_{mm1}`$ for the transition from the $`m`$ to the $`m1`$ state and $`p_{mm+1}`$ from the $`m`$ to the $`m+1`$ state, as follows

$`{\displaystyle \frac{1}{\tau _m}}=p_{mm1}+p_{mm+1}`$ (4)
$`=\{\begin{array}{cc}{\displaystyle \frac{C\mathrm{\Delta }E_m^3}{e^{\mathrm{\Delta }E_m/k_\mathrm{B}T}1}}+{\displaystyle \frac{C\mathrm{\Delta }E_{m+1}^3}{1e^{\mathrm{\Delta }E_{m+1}/k_\mathrm{B}T}}},(\mathrm{for}m>0)\hfill & \\ {\displaystyle \frac{C\mathrm{\Delta }E_{m+1}^3}{e^{\mathrm{\Delta }E_{m+1}/k_\mathrm{B}T}1}}+{\displaystyle \frac{C\mathrm{\Delta }E_m^3}{1e^{\mathrm{\Delta }E_m/k_\mathrm{B}T}}},(\mathrm{for}m<0)\hfill & \end{array}`$ (5)
where $`\mathrm{\Delta }E_m=|E_{m1}E_m|`$. The parameter $`C`$ is given by

$`C={\displaystyle \frac{3}{2\pi \rho v^5\mathrm{}^4}}|m|_{sp}|m\pm 1|^2,`$ (6)
where $`v`$ is the phonon velocity and $`\rho `$ is the specific mass. Thus the spin-lattice relaxation rate is obtained as

$`{\displaystyle \frac{1}{T_1}}`$ $`=`$ $`{\displaystyle \frac{A}{Z}}{\displaystyle \underset{m=10}{\overset{+10}{}}}{\displaystyle \frac{\tau _me^{E_m/k_\mathrm{B}T}}{1+\omega _L^2\tau _m^2}},`$ (7)
where $`A=\gamma _N^2<\mathrm{\Delta }h_\pm ^2>`$.

The relaxation rates calculated by using eq.(7) with fitting parameters $`A=4\times 10^{12}`$ rad/s2 and $`C=5\times 10^5`$ Hz/K3 are shown for several fields in Fig. 1. The experimental results of the temperature and field dependence above 400 mK are well reproduced over six decades by this calculation. This means that relaxation above 400 mK is dominated by thermal fluctuations resulting from the transitions of iron spins between neighboring states due to spin-phonon interactions. The deviations at high temperatures must be caused by the contribution from higher energy levels which are not described by the simplified spin model with $`S=10`$.

The observed $`T_1^1`$ below 300 mK deviates from the calculated values and becomes temperature independent. This may be related to the temperature independent recovery of magnetization that has been observed in SQUIDs below 400 mK at non-level crossing fields,  and may be attributed to the quantum effect.

Next we discuss the stepwise behavior of the echo intensity, shown in Fig. 3. Figure 5 shows schematically the field dependence of the energy levels of the iron spin system described by eq.(1) and the transitions of the states when the field is decreased. In the experiment, when the field is positive at low temperatures, all the spins are in the $`m=10`$ state. When the field is decreased to zero, the level of the $`m=10`$ state crosses with that of the $`m=10`$ state, and the spins in the $`m=10`$ state would change to the $`m=10`$ state by quantum tunneling with a probability $`P_{10,10}`$. Further decrease of the field causes the energy levels to cross again and the transition to occur again. The transition at the point where the levels cross has been calculated by Landau and Zener, and the probability from the $`m`$ state to the $`m^{}`$ state in this case is expressed  as

$`P_{m,m^{}}`$ $`=`$ $`1\mathrm{exp}({\displaystyle \frac{\pi \mathrm{\Delta }_{m,m^{}}^2}{2\mathrm{}g\mu _\mathrm{B}|mm^{}|dH/dt}}),`$ (8)
where $`\mathrm{\Delta }_{m,m^{}}`$ is the tunneling gap at the level crossing. These transitions would induce fluctuations of the local field at proton sites and cause extra relaxation of the nuclear spins.

Immediately after quantum tunneling, the spin state goes to the lower energy levels, following the Boltzman distribution. The life time $`\tau _m`$ is estimated to be less than $`10^7`$ sec from the measured relaxation time and eq.(7). It should be noted that these thermal transitions can occur only between states with the same sign of $`m`$, while quantum tunneling occurs between states with opposite sign of $`m`$.

The measured period of the stepwise behavior for $`\theta =50^{}`$ was 0.32 T, while that for $`\theta =0^{}`$ was 0.22 T. The period of level crossing fields is expressed approximately as $`\mathrm{\Delta }H=D/g\mu _\mathrm{B}\mathrm{cos}\theta `$ and the values are 0.31 T and 0.21 T for $`\theta =50^{}`$ and $`0^{}`$, respectively. The values coincide fairly well with the experimental results. These results clearly indicate that the sudden recovery of the 1H spins is caused by resonant quantum tunneling of the iron magnetization at the level crossing fields.

When the iron spin state in a certain molecule changes through tunneling, the 1H spins in the molecule would be relaxed immediately. The iron spins which have tunneled in zero or negative fields arrive at the $`m=10`$ state and have possibility to tunnel again in the zero or positive fields. If it is assumed that the iron spins which have experienced tunneling once at least contribute the proton relaxation, the echo intensities $`I_i`$ at the monitoring field after the field returning would be expressed as

$`I_1`$ $`=`$ $`I_{01}\{1\left(1P_{10,10}\right)^2\},`$ (9)
for $`H_{10,9}<H_r<0`$, and

$`I_2`$ $`=`$ $`I_{02}\{1\left(1P_{10,10}\right)^2\left(1P_{10,9}\right)^2\},`$ (10)
for $`H_{10,8}<H_r<H_{10,9}`$. The measured sweep rate dependence of echo intensities was fitted to eqs.(9) and (10) with eq.(8), $`\mathrm{\Delta }_{10,10}=3.52\times 10^7`$ K and $`\mathrm{\Delta }_{10,9}=9.66\times 10^7`$ K, as shown in Fig. 4. The agreement is fairly well and these values of $`\mathrm{\Delta }_{10,10}`$ and $`\mathrm{\Delta }_{10,9}`$ are close to those reported from magnetization measurements.

As is found in Fig. 3(b), the echo intensities for $`H_r<0.25`$ T at a sweep rate of 0.9 T/min remain constant and smaller than those for $`dH/dt=0.08`$ T/min. The stepwise recovery was not observed for the lower crossing fields. An intensity of the signals which were measured at the slow sweep rate of 0.08 T/min with $`H_r`$ less than $`0.6`$ T corresponds to the intensity of the signal which recovered completely. This intensity was obtained for the fast rate of 0.9 T/min by sweeping the field up and down three times between $`0.1`$ T and $`0.6`$ T. However, the intensity after sweeping up and down three times at lower fields between $`0.6`$ T and $`1.0`$ T with $`dH/dt=0.9`$ T/min remained small. This suggests that tunneling does not occur at fields lower than $`0.6`$ T, though the reason is not clear.

In conclusion, relaxation rate $`T_1^1`$ of Fe8 above 400 mK are dominated by the thermal fluctuation of the iron magnetization with spin $`S=10`$. The transition of magnetization between the states splitted by $`DS_Z^2`$ are caused by the spin-phonon interaction, and induces the fluctuation at nuclear sites. The rate becomes temperature independent below 300 mK. In this temperature region the stepwise recovery of the echo intensity caused by quantum tunneling at the level crossing fields was observed in the 1H-NMR spectrum. This indicates that the internal field at proton sites fluctuates due to resonant quantum tunneling of iron spins, which is described by the Landau-Zener transition.

We would like to thank Prof. S. Miyashita, Prof. B. Barbara and Dr. K Saito for valuable discussions. We also thank Dr. H. Miyasaka and Dr. H.-C. Chang for their helpful advice on synthesis of the samples.


format: text

text:
# Temperature- and Force-Induced 𝛽-Sheet Unfolding in an Exactly Solvable Model

## Abstract

The stability of a $`\beta `$-sheeted conformation and its transition into a random coil are studied with a 2D lattice biopolymer model. At low temperature and low external force, the polymer folds back and forth on itself and forms a $`\beta `$-sheet. Our analytical calculation and Monte Carlo simulation reveal that a co-operative $`\beta `$-sheet–random coil transition takes places when the temperature or force is increased, with a dramatic decrease in the contact number. These predictions are in good agreement with experiments on titin protein. This transition is not a real phase-transition, indicating that backbone hydrogen-bonding alone is unable to stabilize a distinct $`\beta `$-sheet phase.

preprint: MPIKGF 2001-11

Single-molecule manipulation techniques were used by Gaub and coworkers and several other groups Rief et al. (1997); Carrion-Vazquez et al. (1999) to study the unfolding of proteins in real time. A co-operative $`\mathrm{`}\mathrm{`}`$all-or-nothing” transition between globular and randomly coiled states was observed at a pulling force of $`10^2`$ piconewtons. A stretched coil segment can be modelled as a wormlike chain of persistence length $`10^1`$ nanometers Rief et al. (1997); Carrion-Vazquez et al. (1999), and several (atomic-level) molecular dynamics (MD) simulations have studied the initial stage of unfolding and the final stage of refolding (see, e.g., Lu and Schulten (2000); Paci and Karplus (1999)). Monte Carlo (MC) simulations Socci et al. (1999) of lattice models also have been used to study the unfolding kinetics. The protein used in the experiments, namely titin imunoglobulin domains, is rich in $`\beta `$-sheets Rief et al. (1997); Carrion-Vazquez et al. (1999), and numerical work demonstrated that the transition to a random coil was initiated by the rupture of only a few hydrogen bonds between adjacent $`\beta `$-strands Lu and Schulten (2000). Theoretically, however, it is still unclear why the globule-coil transition is a highly co-operative process, and why there are dramatic differences between the unfolding behavior of helix-dominated and $`\beta `$-sheet-dominated proteins Rief et al. (1999).

An issue of current interest is to study the generic statistical mechanics of biopolymer $`\beta `$-sheet states. In this paper, the deformation and unfolding of $`\beta `$-sheets are addressed within an exactly-solvable 2D model. Because of the attractive interaction between monomers along the chain, the polymer could fold back and forth on itself, thereby maximizing its contact energy. However, a $`\beta `$-sheet, which is a relatively ordered structure and is favorable in energy, is deficit in configurational entropy compared with an irregular coil configuration. Because of this competition between energy and entropy, it is demonstrated by our analytical and simulational calculations that, when the temperature or the external force approaches some critical value, a small variation in temperature or force could lead to a dramatic change in the polymer’s microscopic configuration. This behavior corresponds well both with experimental observations and MD simulations on titin. However, we suggest that this highly co-operative transition is not a phase-transition in the strict sense; the nonspecific hydrogen-bonding interactions between the backbone molecules of a protein alone are hence unlikely to be able to stabilize a distinct $`\beta `$-sheet phase.

The mathematical model is constructed as follows: Consider a biopolymer of $`N`$ identical units (monomers) located on a square lattice (Fig. 1). The distance between two consecutive monomers is $`a_0`$, and to facilitate the calculation the chain is assumed to be directed, namely the bond between monomers $`i`$ and $`i+1`$ can not point to the $`𝐳_0`$ direction (Fig. 1). If two monomers with indices $`i`$ and $`i+m`$ ($`m3`$) occupy nearest neighbor lattice sites, an attractive energy of strength $`ϵ`$ is gained and these two monomers are then $`\mathrm{`}\mathrm{`}`$in contact”. Each monomer can be in one of three states, with respectively zero, one, or two contacts. By folding back and forth on itself, the polymer can gain contact energy, and a segment of $`m`$ columns in which there is at least one pair of monomers in contact between any two adjacent columns is referred to as a $`\mathrm{`}\mathrm{`}`$$`\beta `$-sheet”. Two $`\beta `$-sheet segments are separated by a $`\mathrm{`}\mathrm{`}`$coil segment” of $`m`$ columns ($`m0`$), in which all the monomers are free of contacts (Fig. 1).

We first discuss the free energy of the $`\beta `$-sheet state. Under the action of a force $`f`$ along direction $`𝐳_0`$, the energy of a $`\beta `$-sheet of $`m`$ columns is $`ϵ_{i=1}^{m1}\mathrm{v}(l_i,l_{i+1})mfa_0`$, where $`\mathrm{v}(l_i,l_{i+1})=l_i1`$ (if $`l_il_{i+1}`$) and $`l_{i+1}1`$ otherwise; $`l_i`$ is the number of monomers in the $`i`$-th column. The partition function (PF) $`Z_\beta (n)`$ of a $`\beta `$-sheet of length $`n`$ is

$$Z_\beta (n)=\underset{m=2}{\overset{[n/2]}{}}\underset{l_1=2}{}\mathrm{}\underset{l_m=2}{}\delta _{l_1+\mathrm{}+l_m}^nb^m\underset{i=1}{\overset{m1}{}}a^{\mathrm{v}(l_i,l_{i+1})},$$
(1)
where $`\delta _i^j=1`$ ($`i=j`$) and $`0`$ (otherwise); $`a=\mathrm{exp}(\beta ϵ)`$, $`b=\mathrm{exp}(\beta fa_0)`$, and $`\beta =1/T`$, $`T`$ being the temperature (the Boltzmann constant is set to $`1`$).

There can be many-body interactions among the monomers in a $`\beta `$-sheet and a contact pair can be formed between two remote monomers along the chain, consequently, direct evaluation of $`Z_\beta (n)`$ for a large system is difficult. To circumvent this, we define an auxiliary function $`T_m(\zeta ,k)=_{l_1=2}\mathrm{}_{l_m=2}\zeta ^{kl_1}a^{\mathrm{v}(l_1,l_2)l_1}_{i=2}^{m1}(\zeta /a)^{l_i}a^{\mathrm{v}(l_i,l_{i+1})}(\zeta /a)^{l_m}`$ for $`m2`$ and $`T_1(\zeta ,k)=\zeta ^{2k}/a(a\zeta ^k)`$. One finds

$$T_m(\zeta ,k)=\frac{\zeta ^{2k}T_{m1}(\zeta ,1)}{a(1\zeta ^k)}+\frac{\zeta ^k(1a)T_{m1}(\zeta ,k+1)}{a(a\zeta ^k)(1\zeta ^k)}.$$
(2)
Define $`L_\beta (\zeta ,k)=_{m=1}^{\mathrm{}}b^mT_m(\zeta ,k)`$, and from the recursive relation Eq. (2) one obtains the following two equations: (1) $`L_\beta (\zeta ,2)=A(\zeta )+B(\zeta )L_\beta (\zeta ,1)`$ and (2) $`L_\beta (\zeta ,2)=H(\zeta )+J(\zeta )L_\beta (\zeta ,1)`$. After eliminating $`L_\beta (\zeta ,2)`$ from these two equations, one finds that

$$L_\beta (\zeta ,1)=\frac{A(\zeta )H(\zeta )}{J(\zeta )B(\zeta )}.$$
(3)
In Eq. (3), $`A(\zeta )=\zeta (1\zeta )/(a1)`$, $`B(\zeta )=(a\zeta )(b\zeta ^2+a\zeta a)/b(a1)\zeta `$, $`H=(b/a)[\zeta ^4/(a\zeta ^2)+_{i=3}^{\mathrm{}}_{j=2}^{i1}\varphi (\zeta ,j)\zeta ^{2i}/(a\zeta ^i)]`$, and $`J=(b/a)[\zeta ^4/(1\zeta ^2)+_{i=3}^{\mathrm{}}_{j=2}^{i1}\varphi (\zeta ,j)\zeta ^{2i}/(1\zeta ^i)]`$, where $`\varphi (\zeta ,j)=b(1a)\zeta ^k/a(a\zeta ^k)(1\zeta ^k)`$.

Both $`H(\zeta )`$ and $`J(\zeta )`$ are convergent when $`\zeta [0,1]`$. When $`\mathrm{\Delta }=1\zeta 0^+`$, our analytical and numerical calculations suggest that $`A(\zeta )H(\zeta )\mathrm{exp}(b/a\mathrm{\Delta })f_1(\mathrm{\Delta })`$ and $`J(\zeta )B(\zeta )\mathrm{exp}(b/a\mathrm{\Delta })f_2(\mathrm{\Delta })`$, with $`f_1(\mathrm{\Delta })/f_2(\mathrm{\Delta })b/a(a1)`$. Here $`f_1(\mathrm{\Delta })`$ and $`f_2(\mathrm{\Delta })`$ are two polynomials of $`\mathrm{\Delta }`$ which contain terms with negative powers of $`\mathrm{\Delta }`$.

The grand partition function (GPF) of the $`\beta `$-sheet, $`G_\beta (\zeta )=_{n=4}(\zeta /a)^nZ_\beta (n)`$, is nothing but $`L_\beta (\zeta ,1)b^2\zeta ^2/a(a\zeta )`$. $`G_\beta (\zeta )`$ could be calculated based on Eq. (3) at given $`\zeta `$. It keeps increasing as $`\zeta `$ approaches $`\zeta _\beta (<1)`$ from $`0`$ and is divergent at $`\zeta _\beta `$, where $`\zeta _\beta `$ is the smallest positive root of the equation $`J(\zeta )B(\zeta )=0`$ (see Fig. 2). The free energy density of a $`\beta `$-sheet segment equals to $`\rho _\beta =k_BT\mathrm{ln}(\zeta _\beta /a)`$. The divergence of $`G_\beta `$ at $`\zeta _\beta `$ indicates the absence of a distinct $`\beta `$-sheet phase.

Armed with the partition function of a $`\beta `$-sheet, we proceed to investigate its stability and the configurational transition of the whole polymer system. To this end, one needs to know the free energy of a coil segment footnote1 . The total PF of a coil segment of length $`n`$ under the action of an external force ($`f`$, incorporated in the parameter $`b`$) is expressed as

$$Z_{\mathrm{coil}}(n)=2\underset{m=1}{\overset{n}{}}Z_{\mathrm{coil}}^\mathrm{s}(n,m)b^m+\underset{m=1}{\overset{n}{}}Z_{\mathrm{coil}}^\mathrm{m}(n,m)b^m$$
(4)
and $`Z_{\mathrm{coil}}(0)=1`$. In this expression, $`Z_{\mathrm{coil}}^\mathrm{s}(n,m)`$ is by definition the total number of allowed configurations for a coil of length $`n`$ whose two ends are separated by $`m`$ lattice steps along $`𝐳_0`$ and whose last column has only a single monomer; $`Z_{\mathrm{coil}}^\mathrm{m}(n,m)`$ is defined similarly, but with the last column containing at least two monomers. The coefficient $`2`$ in Eq. (4) appears because the first column of a $`\beta `$-sheet following a coil segment that ends with a column containing a single monomer can have two possible orientations (Fig. 1).

One finds $`Z_{\mathrm{coil}}^\mathrm{s}(n,m)=Z_{\mathrm{coil}}^\mathrm{s}(n1,m1)+Z_{\mathrm{coil}}^\mathrm{m}(n1,m1),`$ and $`Z_{\mathrm{coil}}^\mathrm{m}(n,m)=Z_{\mathrm{coil}}^\mathrm{s}(n,m+1)+Z_{\mathrm{coil}}^\mathrm{s}(n2,m1),`$ with $`Z_{\mathrm{coil}}^\mathrm{s}(0,m)=\delta _m^0`$, $`Z_{\mathrm{coil}}^\mathrm{m}(0,m)=Z_{\mathrm{coil}}^\mathrm{m}(1,m)=0`$, $`Z_{\mathrm{coil}}^\mathrm{m}(2,m)=\delta _m^1`$, and $`Z_{\mathrm{coil}}^\mathrm{m}(n,1)=1\delta _n^1`$. With these preparations the GPF of the coiled state is calculated to be

$$G_{\mathrm{coil}}(\zeta )=\underset{n=0}{\overset{\mathrm{}}{}}\zeta ^na^nZ_{\mathrm{coil}}(n)=\frac{a(a\zeta )(a+b\zeta )}{a^3a^2(b+1)\zeta b^2\zeta ^3},$$
(5)
It follows that the free energy density of a coil segment equals $`\rho _{\mathrm{coil}}=k_BT\mathrm{ln}(\zeta _{\mathrm{coil}}/a)`$. $`\zeta _{\mathrm{coil}}`$ is the point where $`G_{\mathrm{coil}}(\zeta )=+\mathrm{}`$: $`\zeta _{\mathrm{coil}}=\mathrm{ln}[(a/b)((\sqrt{\lambda }+b/2)^{1/3}(\sqrt{\lambda }b/2)^{1/3})]`$ with $`\lambda =(b+1)^3/27+b^2/4`$.

In the lattice model, any configuration of the polymer is a chain of $`\beta `$-sheets and random coils, with the two conformations occurring alternately (Fig. 1). The total PF is then

$`Z(N)`$ $`=`$ $`{\displaystyle \underset{s=0}{\overset{[N/4]}{}}}{\displaystyle \underset{i_0=0}{}}{\displaystyle \underset{j_1=4}{}}{\displaystyle \underset{i_1=0}{}}\mathrm{}{\displaystyle \underset{j_s=4}{}}{\displaystyle \underset{i_s=0}{}}{\displaystyle \underset{j_0=0}{}}`$ (6)
$`\delta _{i_0+j_0+_{k=1}^s(i_k+j_k)}^N{\displaystyle \underset{k=0}{\overset{s}{}}}[Z_\beta (j_k)Z_{\mathrm{coil}}(i_k)],`$
Finally the GPF of the whole system may be calculated:

$$G_{\mathrm{total}}(\zeta )=\underset{N=0}{\overset{\mathrm{}}{}}\zeta ^Na^NZ_{\mathrm{total}}(N)=\frac{[1+G_\beta (\zeta )]G_{\mathrm{coil}}(\zeta )}{1G_\beta (\zeta )G_{\mathrm{coil}}(\zeta )}.$$
(7)
Because of the facts that $`G_\beta (\zeta _\beta )=+\mathrm{}`$ and $`G_{\mathrm{coil}}(\zeta _{\mathrm{coil}})=+\mathrm{}`$, it is always possible to find a root ($`\zeta _{\mathrm{total}}`$) of the equation

$$G_{\mathrm{coil}}(\zeta _{\mathrm{total}})G_\beta (\zeta _{\mathrm{total}})=1$$
(8)
in the range of $`0<\zeta _{\mathrm{total}}<\mathrm{min}(\zeta _\beta ,\zeta _{\mathrm{coil}})`$. Consequently, the free energy density of the whole system is neither that of the $`\beta `$-sheet state nor that of the coil state: $`\rho _{\mathrm{total}}=k_BT\mathrm{ln}(\zeta _{\mathrm{total}}/a)`$. There is no singularity in the total free energy expression, and hence no real phase-transition process. Almost forty years ago, it was revealed by intensified theoretical investigations that the protein $`\alpha `$-helix–coil transition is not a phase-transition, provided that only the main-chain hydrogen-bonding interactions are considered (see, e.g., Lifson (1964); Poland and Scheraga (1970)); we are surprised to find that the same conclusion holds for the $`\beta `$-sheet transition. Although there is no strict phase-transition between $`\beta `$-sheet and random coil in the model, our analytical and simulational results presented below indicate that, when temperature and the external force are low enough the polymer is essentially in the form of a $`\beta `$-sheet, and there is an rapid change in the polymer’s total extension and its total number of contacted pairs when the temperature or the force is increased.

The total contact number $`\sigma `$ (scaled by $`N`$) could be calculated by $`\sigma =1\mathrm{ln}(\zeta _{\mathrm{total}})/\mathrm{ln}a`$ and the total extension $`\mathrm{}`$ (scaled by $`Na_0`$) could be calculated by $`\mathrm{}=b\mathrm{ln}(\zeta _{\mathrm{total}})/b`$. When there is no external force ($`f=0`$), the temperature–contact number profile and the temperature–extension profile are shown in Fig. 3. And in Fig. 4 the relationship between force and extension and between force and contact number are shown at a fixed temperature of $`T=0.591ϵ`$. When the temperature and the force is small, $`\zeta _{\mathrm{total}}`$ approaches $`\zeta _\beta `$ (the divergent radius of $`G_\beta `$) and its value is difficult to be numerically determined in much precision through Eq. (8). Therefore, the analytical calculation in Fig. 3 is shown only for $`T0.8ϵ`$ and that in Fig. 4 for $`f0.6ϵ/a_0`$. In Figs. 3 and 4 the data obtained by a MC simulation is also shown. The Monte Carlo simulation is performed with a polymer size of $`N=1024`$ (the results for $`N=5000`$ are essentially the same as the data shown here). There are three kinds of elementary updates of configurations: (i) one monomer is moved from one column to another column or stacked on the right-hand side of the last column; (ii) the positions of two columns are exchanged; (iii) the direction of one column is changed by $`180^{}`$. Special attention is paid to make sure that ergodicity and detailed balance are obeyed. Each simulational data is obtained through $`10^9`$ MC steps (it takes about $`35`$ hours in a Compaq station). The energetically most favorable state is a $`\beta `$-sheet whose width and height both equal to $`\sqrt{N}`$. Therefore, in the numerical simulation data shown in Figs. 3 and 4 for a chain of $`N=1024`$, one has $`\sigma 0.935`$ and $`\mathrm{}0.0313`$. For an infinite system, $`\sigma 1`$ and $`\mathrm{}0`$ as $`T0`$.

Figure 3 indicates that when $`f=0`$ and $`T<0.75ϵ`$ almost all the monomers of the polymer are in the $`\beta `$-sheet state. At $`T0.75ϵ`$ there is a plateau in the $`T`$$`\sigma `$ and $`T`$$`\mathrm{}`$ profiles, suggesting that a co-operative transition into random coiled structure takes place. During this process, the monomers lose contacts and the polymer is extended. As temperature is further increased, the configuration of the polymer becomes essentially random and there are only small changes of the total extension and the total contact number. Similar co-operative $`\beta `$-sheet–coil transition is observed when the temperature is fixed but the external force keeps increasing. At $`T=0.591ϵ`$ Fig. 4 indicates the transition takes place at $`f0.5ϵ/a_0`$. In the experimental investigations of Gaub and coworkers Rief et al. (1997) a co-operative globule–extended coil transition is observed for the $`\beta `$-sheet dominated protein titin; and in a later experiment of Fernandez and coworkers Carrion-Vazquez et al. (1999) it was suggested that the force-induced denaturation of protein titin should follow the same paths as the temperature-induced denaturation of this protein. The present theoretical work seems to confirm both of these experimental observations: a co-operative transition from $`\beta `$-sheet to extended coil will occur either by increase the environmental temperature or by increase the mechanical force acting on the two ends of the polymer.

Conformational transitions in single biopolymers have attracted the interest of theorists for quite a long time, and the introduction of simple mathematical models have brought a deep understanding of the issues of the protein $`\alpha `$-helix-to-coil transitionPoland and Scheraga (1970), double-stranded DNA melting de Gennes (1969), $`\beta `$-hairpin formation in oligopeptidesGuo et al. (2000a), and RNA secondary structure denaturation by heat de Gennes (1968) or by force Zhou et al. (2001), etc.. The current work extends these efforts to a new model system, namely, $`\beta `$-sheet formation in biopolymers. It suggests that a co-operative $`\beta `$-sheet–coil transition will happen at certain temperature or mechanical stress. To make it more realistic, other effects such as the sequence heterogeneity and the interactions between adjacent layers of $`\beta `$-sheets Lu and Schulten (2000) could be considered. On the other hand, the results presented here show that even without these complications the structural transition of a $`\beta `$-sheet is $`\mathrm{`}\mathrm{`}`$two-state-like”.

Recently, the MC work in Wittkop et al. (1995) suggested that pulling a collapsed polymer in 2D leads to a co-operative transition into an extended random coil, and the MC work in Grassberger and Hsu  suggested this process actually may be a continuous transition. The present effort is also complementary to these MC simulations.

The financial support of the A. v. Humboldt Foundation is acknowledged. The author is grateful to Professor R. Lipowsky for his hospitality, and to Dr. Julian Shillcock for a careful reading of the manuscript and many good suggestions.


text:
# 𝜎 Meson In 𝐽/Ψ Decays

## 1 Introduction

Beijing Spectrometer (BES) at Beijing Electron-Positron Collider (BEPC) recently reports an evidence for the existence of the $`\sigma `$ particle in $`J/\mathrm{\Psi }`$ decays. In the $`\pi ^+\pi ^{}`$ invariant mass spectrum in the process of $`J/\mathrm{\Psi }\pi ^+\pi ^{}\omega `$ they found a low mass enhancement and the detailed analysis strongly favors that the spin-parity is $`𝒪^{++}`$ and the statistical significance for the existence of the $`\sigma `$ particle is about 18 $`\sigma `$ . The BES measured values of the $`\sigma `$ mass and width are:

$`m_\sigma `$ $`=`$ $`390_{36}^{+60}\text{MeV},`$
$`\mathrm{\Gamma }_\sigma `$ $`=`$ $`282_{50}^{+77}\text{MeV},`$ (1)
and the branching ratio is

$$\mathrm{Br}(J/\mathrm{\Psi }\sigma \omega \pi ^+\pi ^{}\omega )=(1.71\pm 3.4\pm 4.3)\times 10^3.$$
(2)
The $`\sigma `$ particle has been absent for many years in the Review of Particle Physics by the Particle Data Group (PDG), however in the recent years there has been a revival of interest in studying the light scalar-isoscalar meson, the $`\sigma `$ particle, as a broad resonanceexperimentally and theoretically. A direct experimental evidence for the $`\sigma `$ meson is reported recently by the Fermilab E791 collaboration  in the D-meson decay process, $`D^+\sigma \pi ^+3\pi `$. Theoretically, $`\sigma `$ meson can play important roles in some problems. In ref. , $`\sigma `$ can be regarded as the dominant contribution to the $`u\overline{u}+d\overline{d}`$ current. Moreover, using $`sigma`$ mass and width from E791, some people investage $`\sigma `$ contribution to $`B3\pi `$ and find it can explain the recent data from CLEO and BABAR collaboration.

In this paper, we study phenomenologically the decay process $`J/\mathrm{\Psi }\sigma \omega \pi ^+\pi ^{}\omega `$ and discuss the couplings $`g_{\sigma \pi \pi }`$ and $`g_{J/\mathrm{\Psi }\sigma \omega }`$. By taking the meson wavefunctions of $`\sigma `$ and $`\omega `$ to be similar to that of $`\pi `$ and $`\rho `$, and using the perturbative QCD technique, we calculate the decay constant $`g_{J/\mathrm{\Psi }\sigma \sigma }^{th}`$. Our theoretical prediction is shown to be within the range of experimental value $`g_{J/\mathrm{\Psi }\sigma \omega }`$.

## 2  Coupling constants $`g_{\sigma \pi \pi }`$ and $`g_{J/\mathrm{\Psi }\sigma \omega }`$

Given the data on the $`\sigma `$ mass and width measured by BES in Eqs.(1) and (2) we study the coupling constants $`g_{\sigma \pi \pi }`$ and $`g_{J/\mathrm{\Psi }\sigma \omega }`$. For a two-body decay of particle $`X`$ into final states $`X_1`$ and $`X_2`$ the decay width is given by

$$\mathrm{\Gamma }(XX_1X_2)=\frac{|𝐩|}{32\pi ^2M_X^2}||^2𝑑\mathrm{\Omega }=\frac{|𝐩|}{8\pi M_X^2}||^2,$$
(3)
where $`|𝐩|`$ is 3-momentum of the final state in the center-of-mass (c.m.) . For $`\sigma \pi ^+\pi ^{}`$, the $`\sigma `$ meson decays $`100\%`$ into $`\pi \pi `$. Furthermore isospin conservation requires for two thirds of the time the final states be the charged pions. So we have:

$$\frac{2}{3}\mathrm{\Gamma }_\sigma ^0=g_{\sigma \pi \pi }^2\frac{1}{8\pi m_\sigma ^2}\sqrt{\frac{m_\sigma ^2}{4}m_\pi ^2}.$$
(4)
Using BES data on the $`\sigma `$’s mass, width, and the experimental value of the $`\pi `$ mass, we obtain:

$`g_{\sigma \pi \pi }`$ $`=`$ $`2.0_{0.19}^{+0.30}\mathrm{GeV}.`$ (5)
This number is surprisingly consistent with the theoretical value of the linear sigma model

$$g_{\sigma \pi \pi }^{linear\sigma }=\frac{\sqrt{2}m_\sigma ^2}{f_\pi }=1.80_{0.30}^{+0.50}\mathrm{GeV},$$
(6)
where $`f_\pi `$ is the pion decay constant. In Eq.(6) the $`\sigma `$ mass is taken from the BES measurement in Eq.(1).

To get the phenomenological value of $`g_{J/\mathrm{\Psi }\sigma \omega }`$, we use the full three-body decay width (for example, see ):

$`\mathrm{\Gamma }(J/\mathrm{\Psi }\sigma \omega \pi ^+\pi ^{}\omega )={\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{2m_{J/\mathrm{\Psi }}}}g_{J/\mathrm{\Psi }\sigma \omega }^2g_{\sigma \pi \pi }^2`$
$`{\displaystyle _{4m_\omega ^2}^{(m_{J/\mathrm{\Psi }}m_\omega )^2}}{\displaystyle \frac{d\chi ^2}{2\pi }}{\displaystyle \frac{1}{8\pi }}\lambda ^{1/2}(1,{\displaystyle \frac{\chi ^2}{m_{J/\mathrm{\Psi }}^2}},{\displaystyle \frac{m_\omega ^2}{m_{J/\mathrm{\Psi }}^2}})`$
$`\times {\displaystyle \frac{1}{8\pi }}\lambda ^{1/2}(1,{\displaystyle \frac{m_\pi ^2}{\chi ^2}},{\displaystyle \frac{m_\pi ^2}{\chi ^2}}){\displaystyle \frac{1}{(\chi ^2m_\sigma ^2)^2+\mathrm{\Gamma }_\sigma (\chi )^2m_\sigma ^2}},`$ (7)
where factor $`1/(8\pi )\times \lambda ^{1/2}`$ is the phase space integral of the corresponding two-body decay subprocess and

$$\lambda (a,b,c)=a^2+b^2+c^22ab2bc2ca.$$
(8)
And in Eq. (7)

$$\mathrm{\Gamma }_\sigma (\chi )\mathrm{\Gamma }_\sigma ^0\times \left(m_\sigma /\chi \right)\left(p^{}(\chi )/p^{}(m_\sigma )\right)$$
(9)
is the co-moving resonance width where $`p^{}(\chi )=\sqrt{\chi ^2/4m_\pi ^2}`$ and $`\mathrm{\Gamma }_\sigma ^0=(370_{36}^{+60})`$ MeV the experimental value in Eq.(1).

With the experimental values of $`g_{\sigma \pi \pi }`$ in Eq.(5) and the branching ratio in Eq.(2), solving Eq. (7) gives rise to

$`g_{J/\mathrm{\Psi }\sigma \omega }^{exp}`$ $`=`$ $`7.3_{1.9}^{+2.6}\mathrm{MeV}.`$ (10)
## 3 $`g_{J/\mathrm{\Psi }\sigma \omega }^{th}`$ calculated by the perturbative QCD

In this section we follow closely the calculation of the exclusive decays of $`\mathrm{{\rm Y}}`$ in ref.  to study the decay of $`J/\mathrm{\Psi }\sigma \omega `$. The decay amplitude consists of two parts. One is the hard decay amplitude of the three-gluon modes and another is the bound-state matrix elements of outgoing mesons.

In general the decay amplitude of $`J/\mathrm{\Psi }\sigma \omega `$ can be written as

$$=\mathrm{\Psi }_{J/\mathrm{\Psi }}(0)_1^1𝑑x_1^1𝑑y\varphi _\sigma ^{}(x,M_c^2)\varphi _\omega ^{}(y,M_c^2)T_h(x,y,M_c^2),$$
(11)
where $`x=x_1x_2,y=y_1y_2`$ and $`x_i`$, $`y_i`$ are the constituent’s fractional longitudinal momenta which satisfy $`\mathrm{\Sigma }x_i=1`$ and $`\mathrm{\Sigma }y_i=1`$. In Eq.(11) $`\mathrm{\Psi }_{J/\mathrm{\Psi }}(0)`$ is a non-relativistic approximate wave function of the $`J/\mathrm{\Psi }`$ meson, $`\varphi _\sigma (x,M_c^2)`$ and $`\varphi _\omega (y,M_c^2)`$ are the distribution amplitude of the $`\sigma `$ and $`\omega `$ wave functions respectively and $`M_c`$ is the charm quark mass.

$`T_h(x,y,M_c^2)`$ is the hard decay amplitude of the charm quark pairs into two light quark-antiquark pairs, which is defined by

$$T_h(x,y,M_c^2)=\frac{d^4l}{(2\pi )^4}𝒯_h(x,y,l,M_c^2).$$
(12)
There are twelve Feynmann diagrams shown in Fig. 1 contributing to $`T_h`$, which are all explicitly shown in Ref. . Using the Landau rules, it is found that both infrared and collinear divergences exist in every diagram. Fortunately, through a use of color neutrality and collinear Ward identities, one can prove that when summing up all of the diagrams, the divergences cancel .

For quarkonia $`J/\mathrm{\Psi }`$, the non-relativistic approximate wave function $`\mathrm{\Psi }_{J/\mathrm{\Psi }}(0)`$ can be determined by the decay $`J/\mathrm{\Psi }e^+e^{}`$ ,

$$\mathrm{\Gamma }(J/\mathrm{\Psi }e^+e^{})=\frac{16\pi \alpha _s^2e^2}{M_{J/\psi }^2}|\mathrm{\Psi }(0)|^2.$$
(13)
According to the Brodsky-Huang-Lepage , the light-cone wavefunction of a hadron is essentially determined by the off-shell energy variable. Thus, for the light scalar meson, $`\sigma `$, we assume the light-cone wave function to be the same as that of the pion

$$\psi _\sigma (x_i,𝐤_{})=A_\sigma \mathrm{exp}[b_\sigma ^2(\frac{𝐤_{}^\mathrm{𝟐}+m_u^2}{x_1}+\frac{𝐤_{}^\mathrm{𝟐}+m_u^2}{x_2})],$$
(14)
where $`𝐤_{}`$ is the relative transverse momentum of the final meson, $`m_u`$ is $`u`$ quark mass and $`x_i`$ are the constituent’s fractional longitudinal momenta. $`A_\sigma `$ and $`b_\sigma `$ are two free parameters which are taken to be $`A_\sigma =A_\pi 32\mathrm{G}\mathrm{e}\mathrm{V}^1`$ and $`b_\sigma ^2=b_\pi ^20.84\mathrm{GeV}^2`$ . For the wave function of the vector meson, $`\omega `$, we assume it to be the same as that of the $`\rho `$ meson

$$\psi _\omega (y_i,𝐤_{})=A_\omega \mathrm{exp}[b_\omega ^2(\frac{𝐤_{}^\mathrm{𝟐}+m_u^2}{y_1}+\frac{𝐤_{}^\mathrm{𝟐}+m_u^2}{y_2})],$$
(15)
where $`A_\omega `$ and $`b_\omega `$ are taken to be $`A_\omega =A_\rho 30\mathrm{G}\mathrm{e}\mathrm{V}^1`$ and $`b_\omega ^2=b_\rho ^20.55\mathrm{GeV}^2`$, which can be determined from two constraints.

The distribution amplitude $`\varphi `$ in Eq.(11) is defined as

$$\varphi (x_1,x_2,M_c^2)=_0^{M_c^2}\frac{dk_{}^2}{16\pi ^2}\psi (x_1,x_2,k_{}^2).$$
(16)
Making use of the wave function of $`\sigma `$ and $`\omega `$ in Eqs.(14) and (15) we have:

$`\varphi _\sigma (x,M_c^2)={\displaystyle \frac{A_\sigma (1x^2)}{64\pi ^2b_\sigma ^2}}[\mathrm{exp}(4b_\sigma ^2{\displaystyle \frac{m_u^2}{1x^2}})\mathrm{exp}(4b_\sigma ^2{\displaystyle \frac{M_c^2+m_u^2}{1x^2}})],`$ (17)
$`\varphi _\omega (y,M_c^2)={\displaystyle \frac{A_\omega (1y^2)}{64\pi ^2b_\omega ^2}}[\mathrm{exp}(4b_\omega ^2{\displaystyle \frac{m_u^2}{1y^2}})\mathrm{exp}(4b_\omega ^2{\displaystyle \frac{M_c^2+m_u^2}{1y^2}})],`$ (18)
where we have used $`x=x_1x_2,y=y_1y_2`$, $`\mathrm{\Sigma }x_i=1`$ and $`\mathrm{\Sigma }y_i=1`$. In deriving the integration of Eq.(16), we have ignored the QCD evolution on the distribution amplitude since the charm quark $`M_c`$ is not large and the evolution effects are not significant in this energy range. Thus the distribution amplitude is essentially determined by the non-perturbative model. In Fig.2 we show the distribution amplitudes of $`\sigma `$, $`\omega `$.

With the values of the parameters $`A_\pi ,b_\pi ^2,A_\rho ,b_\rho ^2`$ given above , we obtain the coupling constant $`g_{J/\mathrm{\Psi }\sigma \omega }^{th}`$ responsible for the decay $`J/\mathrm{\Psi }\sigma \omega `$ ,

$`g_{J/\mathrm{\Psi }\sigma \omega }^{th}=10.7\mathrm{MeV}.`$ (19)
One can see that our theoretical prediction on $`g_{J/\mathrm{\Psi }\sigma \omega }^{th}`$ above agrees within the errors with the experimental value $`g_{J/\mathrm{\Psi }\sigma \omega }^{exp}`$ in Eq.(10).

## 4 Conclusion

In this paper we have used the new BES experimental evidence for the $`\sigma `$ particle and studied its properties in the process of $`J/\mathrm{\Psi }\sigma \omega \pi \pi \omega `$. We have obtained the phenomelogical values of the coupling constants $`g_{\sigma \pi \pi }`$ and $`g_{J/\mathrm{\Psi }\sigma \omega }`$, then the theoretical prediction on the $`g_{\sigma \pi \pi }^{linear\sigma }`$ in the linear $`\sigma `$ model. We show that linear $`\sigma `$ model gives rise to a reasonable description of the $`\sigma `$ decay into $`\pi `$’s. We have calculated $`g_{J/\mathrm{\Psi }\sigma \omega }`$ by the perturbative QCD and in this approach we take the wave functions of the $`\sigma `$ and $`\omega `$ to be similar to that of $`\pi `$ and $`\rho `$. Our theoretical predictions and the experimental values of the $`g_{J/\mathrm{\Psi }\sigma \omega }`$ are shown to be consistent.

## Acknowledgments

One of the authors (W. J. Huo) acknowledges support from the Chinese Postdoctoral Science Foundation and CAS K.C. Wong Postdoctoral Research Award Fund. We thank Zhi-Peng Zheng and N. Wu for useful discussions. This work is supported in part by the NSF of China.


text:
# Formation Pumping of Molecular Hydrogen in the Messier 17 Photodissociation Region

## 1 INTRODUCTION

### 1.1 Messier 17

Messier 17 is one of the nearest regions of massive star formation to us, situated $`1.3`$ kpc away (Hanson, Howarth & Conti 1997). It is particularly notable for the bright HII region which attracted Messier’s attention (and also known as the Omega Nebula). Within it lies an obscured cluster of young stars. Thirteen OB stars have been spectroscopically identified (Hanson et al. 1997), all with between 4 and 15 magnitudes of extinction at optical wavelengths and estimated to be $`10^6`$ years old. The HII region is adjacent to a molecular cloud complex, in particular the core M 17 SW, which contains more than $`10^4\mathrm{M}_{\mathrm{}}`$ of material (Lada 1976, Thronson & Lada 1983). Intense far–IR emission arises from the core, with a luminosity of $`6\times 10^6\mathrm{L}_{\mathrm{}}`$ (Harper et al. 1976, Wilson et al. 1979). While this originates from dust exposed to far–UV radiation, the core appears to be externally heated, by the stars that excite the HII region (Gatley et al. 1979), rather than internally. It forms a photodissociation region (PDR, Tielens & Hollenbach 1985). The method of excitation of the hydrogen molecules in it are the subject of this paper.

### 1.2 Molecular Hydrogen Observations in M 17

M 17 was first mapped in H2 through the 2.12$`\mu `$m 1–0 S(1) line by Meadows (1986) and Gatley & Kaifu (1987). Their maps reveal two obvious regions of emission, a bar running across the north, and another running south to west. They are known as the ‘Northern Bar’ and the ‘South Western Bar’, respectively. They closely follow two emission ridges seen in radio continuum from 1.3 to 21 cm (Felli et al. 1984). These bars are PDRs, interface regions between the ionized and the molecular gas. The Northern Bar also follows an optically visible bar of the HII region, suggesting that the extinction to it is relatively small. The South Western Bar, however, is optically obscured. PDR modelling of far–IR and sub–mm line emission features from the SW–Bar (Burton, Hollenbach & Tielens 1990, Meixner et al. 1992) indicates that it is a clumpy region, containing gas with a mix of densities, ranging from $`10^3`$ to $`10^7`$$`\mathrm{cm}^3`$, exposed to a far–UV radiation field of $`10^4`$ times the average interstellar radiation field.

Low spatial resolution (20 arcsec beam) near–IR spectroscopy of the 1–0 and 2–1 S(1) H2 lines (Tanaka et al. 1989) in the N–Bar found a line ratio of $`2`$. This is the ratio expected if the excitation is via a fluorescent cascade through the vibrational-rotational levels of the ground electronic state (Black & Dalgarno 1976). The molecule is first excited by far–UV photons to an excited electronic level, from which it decays. Chrysostomou et al. (1992, 1993) imaged the N–Bar with arcsecond spatial resolution in the H2 1–0 S(1) line, as well as in the 3.29$`\mu `$m PAH emission feature and the hydrogen Br$`\gamma `$ line at 2.17$`\mu `$m. They also obtained long slit spectroscopy across the N–Bar from 2.0–2.5$`\mu `$m, measuring 16 lines of H2, as well as several HII region emission lines. The spatial coincidence found for the PAH and H2 emission features confirmed the fluorescent origin, and also suggested that most of the gas had density $`>10^5`$$`\mathrm{cm}^3`$. They found the 1–0/2–1 S(1) ratio to be $`3`$, and to be constant along the 60 arcsec long slit. This not only suggests that the H2 is slightly thermalised (i.e. the density is about the critical density, $`10^5`$$`\mathrm{cm}^3`$), but that the physical conditions along the N–Bar are constant.

### 1.3 Molecular Hydrogen Excitation in PDRs

In dense PDRs, with a strong far–UV radiation field, self-shielding can bring the H/H2 dissociation front to optical depths $`A_v<1`$ mag from the ionized surface of the cloud. Here the temperature can reach 1,000–2,000 K, allowing the v=1 level of H2 to be thermalised. H2 line ratios, for instance the 1–0/2–1 S(1) ratio, then increase from the pure fluorescent value of 2, towards the thermal value of 10 observed in shocks and in gas at $`2000`$ K (see Sternberg & Dalgarno 1989, Burton, Hollenbach & Tielens 1990). The actual line ratio depends on the density of the gas. If this is not constant, then the ratio will likely vary with position in a source. However calculating its value depends on an accurate knowledge of the collisional excitation rate coefficients for H2, which are poorly known. Hence the quantitative application of “collisional fluorescence” models to the interpretation of H2 lines ratios in PDRs requires some care.

A component of the fluorescent emission from H2 must arise due to the formation of the hydrogen molecules themselves. For photo-excitation of H2 molecules by far–UV photons, approximately 15% lead to dissociation. However the hydrogen molecules can be reformed, through the recombination of two hydrogen atoms on the surfaces of dust grains (Hollenbach & Salpeter 1971). In the steady state, the rate of formation of H2 molecules balances its rate of photodissociation. The 4.5 eV bond energy of H2 is released on formation. Some of this is used in escaping from the surface of the dust grain, some goes into translational kinetic energy in the molecule, and the rest goes into internal energy. The new molecules can then emerge in a excited vibrational-rotational state, from which radiative decay will occur. Thus a formation component to the H2 spectrum is expected. However, the signature of this spectrum is unclear.

As Black & Dalgarno (1976) first realised, when calculating the infrared spectrum of fluorescent H2, the contribution from molecule formation may be important. In their original model they assumed equipartition of the binding energy released on formation. They divided this equally between the internal energy of the molecule, its translational energy on escape from the grain surface, and the internal energy imparted to the grain lattice. They further assumed that the 1.5 eV provided as internal energy to the molecule was spread with a Boltzmann distribution through the vibrational-rotational levels. However, other assumptions are possible. For instance, Hunter & Watson (1978) argue that H2 molecules are released in rotationally hot, vibrationally cold states (i.e. high–J $`7`$, low–v). A model by Duley & Williams (1986), on the other hand, although agreeing that the molecules should form hot, argued that they would appear in just the opposite combination of states (i.e. low–J, high–v, with v=6 likely). Le Bourlot et al. (1995) have investigated the infrared spectrum produced for pure formation pumping of H2, under a variety of formation models. For instance, as might be expected, the intensities of v=6 lines are found to be considerably greater under the Duley & Williams (1986) model than in other models. However, since this model does not include the fluorescent cascade component to the emission, which will generally dominate the intensity of the H2 lines, it is hard to use it to undertake a quantitative comparison with data.

### 1.4 Previous Reports of the Detection of Formation Pumping for H2

There have been few reports made of the signature for H2 formation being observed in a spectrum. Wagenblast (1992) considers UV absorption spectra of three nearby diffuse clouds, along lines of sight to background stars. He finds that the populations in the excited, pure-rotational levels, v=0, J=5, 6 & 7, along the sight lines, cannot be produced by UV or thermal excitation. Assuming that molecule formation takes place in two adjacent rotational levels within a particular vibrational state (one for ortho-H2, or odd–J, and the other for para-H2, or even–J), he calculates the possible pairings that could account for the observed ratios of the three lines. He found that H2 would need to be formed in a rotationally hot state ($`\mathrm{J}7`$), but with a range of vibrational states possible, up to v=11. Federman et al. (1995) explore this further with additional absorption measurements from the v=3 level in the source $`\zeta `$ Oph, though do not find any fit to the data particularly satisfactory.

Mouri & Taniguchi (1995) consider the 1.5–2.5$`\mu `$m spectrum of the starburst galaxy NGC 6240, where a number of H2 lines are evident. Based on the relative strengths of the 1–0 S(7) and S(9) lines, compared to the 1–0 S(1) line, they argue that formation pumping, via associative detachment of H and H- to form H2, provides an important contribution to the line intensity. However the data on which this is based have low spectral resolution, and suffers from blending. Moreover, the data are fit with a model which contains several components to the 1–0 S(1) line intensity; formation (10%), fluorescence (20%) and thermal excitation (70%). In addition, the v=6–4 Q(1) line at 1.64$`\mu `$m, observed by Elston & Maloney (1990) in this source, and arising from the same level as the 6–4 O(3) line which we report on in this paper, is weak compared to the 1–0 S(1) line, less than 10% its intensity. It has low signal to noise in the spectrum. Hence, claims for a formation signature in the data need to be regarded with caution.

Measurements of over 30 high-excitation H2 lines in the reflection nebula NGC 2023 by Burton et al. (1992) lead these authors to speculate on whether there was a formation component to the emission from several lines from the v=4 level. This was based on an excess in the level column density distribution for these lines compared to expectations for pure fluorescent emission (and see also McCartney et al. 1999 for an extension to v=6 in this source, for which the excess may still be apparent). However, the signal to noise is not sufficient to be certain that the excess is real, and moreover, might be mistaken for ortho-to-para ratio variations between lines in different vibrational levels.

This paper presents images of emission lines from the v=1, 2 and 6 levels in M 17, as well as in the hydrogen Br $`\gamma `$ line. This includes the first map ever obtained of an H2 emission line from the v=6 level. While the bulk of the H2 emission is clearly fluorescent in origin, we argue, based on the different morphology for the v=6–4 O(3) line to the lower excitation lines, that formation pumping of H2 provides a significant component to its flux.

## 2 OBSERVATIONS AND DATA REDUCTION

Two fields in the source Messier 17 were observed using the IRIS 1–2.5$`\mu `$m camera, in conjunction with the University of New South Wales Infrared Fabry-Pérot etalon (UNSWIRF, Ryder et al. 1998), on the 3.9-m Anglo-Australian Telescope (AAT). The data were obtained on 1997 July 21–23. The fields were chosen to centre on the H2 emission in the northern (N) and south-western (SW) ionization bars of M 17, at $`18^h20^m42.8^s`$, $`16^{}8^{}23^{\prime \prime }`$ and $`18^h20^m25.3^s`$, $`16^{}11^{}27^{\prime \prime }`$ (J2000), respectively. These are shown in Fig. 1, overlaid on a 3-colour near–IR image of M 17 (1–2.5$`\mu `$m) obtained by Lada et al. (1991). This image clearly shows the two ionization bars, with blue indicating the ionized gas (through Paschen lines of hydrogen in the 1.25$`\mu `$m band) and red mostly indicating the photodissociation region (PDR, through fluorescently-excited H2 lines in the 2.2$`\mu `$m band).

The UNSWIRF etalon has a FWHM spectral resolution of $`75\mathrm{km}\mathrm{s}^1`$, a pixel size of 0.77 arcsec and a 100 arcsec circular field of view. It is scanned through a spectral line of interest in order to obtain an emission line image with minimal contamination from any continuum radiation present in a source.

The H2 1–0 S(1) (2.1218$`\mu `$m), 2–1 S(1) (2.2233$`\mu `$m), 1–0 S(7) (1.7480$`\mu `$m) and 6–4 O(3) (1.7326$`\mu `$m) lines, and the hydrogen Br$`\gamma `$ (n=7–4, 2.1661$`\mu `$m) and Br10 (n=10–4, 1.7367$`\mu `$m) lines were imaged in the N field, and just the 2$`\mu `$m–band lines in the SW field. Each of the three 2$`\mu `$m lines was imaged, using appropriate 1% width blocking filters, by stepping the etalon through the relevant plate spacings for the line of interest. All three 1.7$`\mu `$m lines were imaged with the same 1.5% width blocking filter. For the weaker 6–4 O(3) line, which had never been imaged before, the spectral scan was defined to also include the nearby Br10 line. The hydrogen line provided a wavelength reference point to ensure the correct plate spacing was used for the 6–4 O(3) line.

The observing sequence with UNSWIRF included 3 on-line settings for the 2$`\mu `$m lines, and 5 on-line settings for the 1.7$`\mu `$m lines, equally spaced in 39 $`\text{km s}^1`$ steps about the respective line centres. An ‘off-line’ setting was chosen to provide continuum subtraction. Sky frames were also taken for each etalon spacing. An integration time of two minutes was used for each setting at 2$`\mu `$m, and five minutes at 1.7$`\mu `$m. The corresponding sky frames were taken immediately after each source frame for each separate Fabry-Pérot etalon setting (and located 5 arcmin N for the N–field and 5 arcmin W for the SW–field). One such sequence of frames was taken for the 1–0 S(1) and Br$`\gamma `$ lines, four for the 2–1 S(1) line (two in the SW field) and two for the 1.7$`\mu `$m lines. Each repeated sequence was centred on a slightly different source position, to minimise the effect of any bad pixels.

The stars BS7330 (K=4.97 mag) and BS6748 (K=4.57 mag) were used as flux standards, and were imaged at each of the etalon spacings. A diffused dome lamp provided a flat field for each etalon spacing, and an arc lamp was scanned through a free spectral range in order to wavelength calibrate the Fabry-Pérot etalon response for each pixel of the array.

Data reduction was through a custom software package using IRAF111Image Reduction and Analysis Facility (www.iraf.noao.edu). Frames are linearised, flat-fielded using a dome flat at the appropriate plate setting, sky-subtracted, shifted to align stars in each frame, smoothed and the off-line frame subtracted from each on-line frame (having been appropriately scaled to minimise residuals from the subtraction process). Stacking of the frames yields a data cube, which is then fitted pixel-by-pixel with the instrumental profile (a Lorentzian) to yield a line image. Results from different sequences on the same line were then coadded. The absolute accuracy in flux calibration is typically around 30% for each image, and line ratios determined from them can thus have a (constant) scaling error of up to a factor of two. However, the error in relative line ratios determined between any two pixels is much smaller than this, depending only on the S/N each line has been measured with, and not on the accuracy of the absolute calibration.

The line centre for each pixel is determined from the plate spacing found for peak of the fitted profile, with an accuracy of $`10\mathrm{km}\mathrm{s}^1`$ (depending on the intensity). Within these errors, the emission velocity for each line across the fields was found to be the same, as expected for a photodissociation region, where gas motions are very much less than the spectral resolution of the data.

## 3 RESULTS

### 3.1 Northern Bar

#### 3.1.1 H2 1–0 S(1), 1–0 S(7) and 2–1 S(1) Lines

Fig. 2 shows images of the H2 1–0 S(1), 2–1 S(1), 6–4 O(3) and 1–0 S(7) lines from the N field. Line fluxes and line ratios for selected regions are shown in Table 1. Fig. 3 compares, for the N field, the distribution of the H2 1–0 S(1) line (in greyscale) to that of the 2–1 S(1), 6–4 O(3), 1–0 S(7) and H Br$`\gamma `$ lines (in contours), respectively. The ionization front between the HII and H2 emitting gas can be seen, though since foreground ionized gas envelops the whole region, it renders its location hard to discern.

Giving regard to the 30% error in absolute calibration, discussed in §2, the intensity of the 1–0 S(1) line at the peak, $`5\times 10^{13}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$ in an $`18\times 15`$ arcsec aperture, compares reasonably with the flux measured in the N–Bar by both Tanaka et al. (1989) ($`3\times 10^{13}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$ in a 20 arcsec beam) and by Chrysostomou et al. (1993) ($`3.5\times 10^{13}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$ in long slit with 180 arcsec2 effective area). The total 1–0 S(1) line emission from the N–Bar field is $`5.3\times 10^{12}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$, and the brightest emission peaks at $`3.2\times 10^{15}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2\mathrm{arcsec}^2`$ at $`18^h20^m42.3^s`$, $`16^{}08^{}51^{\prime \prime }`$ (J2000).

The overall morphology of the line emission is one of clumpy filaments, with two roughly parallel features running SE–NW across the field. They are embedded in diffuse H2 emission, extending to the NE of the field. The morphology is similar in the 1–0 S(1), 2–1 S(1) and 1–0 S(7) line images. This indicates that the 1–0/2–1 S(1) ratio is reasonably constant across the field.

The line ratio variations are quantified in Table 1. The 1–0 S(1)/2–1 S(1) ratio varies from 1.3 to 2.1, with a mean of 1.6. The 1–0 S(1)/1–0 S(7) ratio varies from 3.9 to 4.9, with a mean of 4.0.

Chrysostomou et al. (1993) also found little variation in the 1–0 S(1)/2–1 S(1) ratio, with a value of $`3\pm 0.5`$ over most of the region they observed. Tanaka et al. (1989) determined 1.9 for this ratio. Our own determination is 1.6. However, given the uncertainty in absolute calibration for measurements made through Fabry-Pérot etalons, these determinations are consistent with one another. The value of 3 is more likely to be correct, though, since measurements through a long slit have better relative accuracy for line ratio determinations. This suggests our 1–0 S(1)/2–1 S(1) line ratios might need to be scaled by $`1.8`$ for the N–Bar.

The similarity of the 1–0 S(1) and 1–0 S(7) line images indicates that any differential extinction between their emitting wavelengths, 1.75 and 2.12$`\mu `$m, must be uniform (or minimal). This is consistent with the estimate of the extinction of 0.05 mag at 2.2$`\mu `$m found by Chrysostomou et al. (1993), based on hydrogen recombination lines at 1.09 and 2.16$`\mu `$m. For this value of extinction, the differential extinction between 1.7 and 2.1$`\mu `$m is indeed negligible.

#### 3.1.2 H2 v=6–4 O(3) Line

For the v=6–4 O(3) line, while the emission clearly arises from the same region as the lower excitation lines, its distribution within that region is quite different. The line is considerably brighter in the NE of the two emitting filaments, contrary to what is seen in the v=1 and v=2 line images. The number of distinct clumps is also greater in the 6–4 O(3) image.

Since the v=6–4 O(3) line is emitted at a similar frequency to the 1–0 S(7) line there will not be any differential extinction between these two lines. Thus, the difference in morphology between the 6–4 O(3) line and the 1–0 S(1) line reflects differences in their excitation. It is not the result of differential extinction between the emitting wavelengths of these two lines.

The 1–0 S(1)/6–4 O(3) ratio varies from 5 to 23 in the N–Bar, with a mean value of 11 (see Table 1). This is a factor of two variation from the mean ratio, compared to less than 30% variation measured for the ratio of the 1–0 S(7) and 2–1 S(1) lines with the 1–0 S(1) line.

### 3.2 South Western Bar

In Fig. 4 are shown the H2 1–0 and 2–1 S(1) lines and the H Br$`\gamma `$ line for the SW field, together with an overlay of the Br$`\gamma `$ line on the 1–0 S(1) line. Line fluxes are again listed in Table 1. The H Br10 line is not shown, but is similar to the Br$`\gamma `$ line in appearance. Again the H2 emission is distributed along clumpy filaments, and (giving regard to the lower S/N in the v=2 image) the v=1 and 2 lines show similar morphologies. A clear ionization front is evident between the Br$`\gamma `$ and 1–0 S(1) images, with the brightest emission filaments running parallel and separated by $`5`$ arcsec.

The total 1–0 S(1) line flux from the SW field is $`2.6\times 10^{12}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$. The flux peaks at $`1.9\times 10^{15}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2\mathrm{arcsec}^2`$ at $`18^h20^m23.1^s`$, $`16^{}11^{}16^{\prime \prime }`$ (J2000). Both these fluxes are about half the corresponding values in the N–Bar. This is probably mostly due to the extra extinction in the optically obscured SW–Bar.

The 1–0 S(1)/2–1 S(1) line ratio varies between 3 and 8, with a mean value of 4.4. This is both higher than in the N–Bar, and shows significantly more variation with position. This line ratio cannot be considered as constant over the field, unlike in the N–Bar.

## 4 DISCUSSION

### 4.1 PDR Emission and Gas Density

As has been discussed in several previous papers, the molecular hydrogen line emission in M 17 is dominated by UV fluorescence (e.g. Tanaka et al. 1989, Chrysostomou et al. 1992, Chrysostomou et al. 1993). This is evident through the morphology, its proximity to the ionization front of an HII region, the similarity to UV–excited 3.3$`\mu `$m PAH emission, and through the measurements made of H2 vibrational line ratios. Our data confirm this result. In particular, the strength of the v=2–1 S(1) line compared to the v=1–0 S(1) line is indicative of a non-thermal excitation method such as UV fluorescence (e.g. Black & Dalgarno 1976, Black & van Dishoeck 1987).

In the SW–Bar the line ratio rises above the pure fluorescent value, indicating that the density of some of the gas there is greater than critical ($`10^4`$$`10^5`$$`\mathrm{cm}^2`$, depending on uncertain values for the H2 collisional excitation rates). In such cases collisions can re-populate the v=1 level so that the ratio of the v=1–0 to the v=2–1 S(1) lines can appear thermal (i.e. $`10`$, see Sternberg & Dalgarno 1989, Burton, Hollenbach & Tielens 1990). However when such “collisional fluorescence” is seen, density variations within the emitting region invariably give rise to significant variations in the value of this line ratio within the source (see Ryder et al. 1998, Allen et al. 1999). Density variations must occur in the SW–Bar of M 17 too.

However the bulk of the gas in the SW–Bar must have density no more than critical. This can be determined from the separation between the ionization front and the excited H2 there, $`5`$ arcsec. Assuming this corresponds to an extinction of $`A_v1`$ mag (i.e. $`\mathrm{N}10^{21}`$ H2 molecules cm -2), it implies that the average H2 number density is $`\mathrm{n}1\times 10^4\mathrm{cm}^3`$ in the SW–Bar. Given the high density component that is also there, the medium must therefore be a clumpy one.

In the N–Bar we will take the line ratio to be 3, as determined by Chrysostomou et at. (1993). This is only a little higher than the pure fluorescent value, and suggests that the density here is about equal to the critical density. The constancy of the ratio across the N–Bar also suggests that there is also little variation in the density. This is in contrast to the clumpy SW–Bar.

### 4.2 Molecular Hydrogen Formation Pumping

The most interesting result from this work is the v=6–4 O(3) emission line image of the N–bar, shown in Fig. 2. This is the first time an image in such a high-excitation line of H2, 31,000 K above ground state, has been obtained. While the line arises from the same regions as the v=1 and v=2 lines of H2, its distribution within them is clearly different.

Since all the lines arise from ortho (odd–J) states, the difference between them is also unlikely to result from ortho-to-para ratio variations, unless these both vary with position in the source, as well as between vibrational levels.

It is also hard to see how the 6–4 O(3) line could have been mis-identified. It cannot be a line emitted from the HII region due to the completely different morphology from the hydrogen lines. Since the line is emitted from the same regions as the lower excitation H2 lines it also suggests that it is a PDR line, and not from the HII region. Moreover, the nearby H Br10 line provided a wavelength reference point for scanning the Fabry-Pérot etalon across the correct plate spacings for the 6–4 O(3) line.

The v=6–4 O(3) line emission almost certainly arises from a region where the H2 emission is dominated by UV–fluorescence. It cannot be shocked or X–ray excited since the high energy of the (v,J) = (6,1) upper level (31,000 K) would not be significantly populated by thermal means. The 6–4 O(3) line would be $`10^4`$ of the strength of the 1–0 S(1) line for thermal excitation at 2000 K, considerably weaker than the nearby 1–0 S(7) line, which would be about 10% of the 1–0 S(1) line strength. In fluorescent models, however, the v=6–4 O(3) line intensity is significant; for instance in Black & van Dishoeck’s Model 14 it is 31% of the value of the 1–0 S(1) line. In the same model the 2–1 S(1) line is 56% the strength, and the 1–0 S(7) line 18%, of the 1–0 S(1) line intensity (neglecting any differential extinction between the emitting wavelengths). This is broadly consistent with the data, for which the mean fluxes of the 2–1 S(1), 1–0 S(7) and 6–4 O(3) lines are $``$60%, 25% and 10% of the 1–0 S(1) line, respectively. However the specific predictions for the strength of the lines depend upon an additional excitation mechanism to the fluorescent cascade, formation pumping.

For every H2 molecule that is fluorescently excited, 15% lead to photodissociation (e.g. Draine & Bertoldi 1996). In steady state PDR models, molecule destruction is balanced by molecule formation, which is believed to occur on grain surfaces (e.g. Hollenbach & Salpeter 1971). The newly formed H2 molecules are released in an excited vibrational-rotational state. However, there is little hard evidence to suggest in what state, or distribution of states, this might be. Black & van Dishoeck (1987) consider three possible formation models. Their ‘standard’ Model 14 assumes that 1.5 eV of the binding energy is distributed in a Boltzmann distribution through the energy levels (after Black & Dalgarno 1976). A second model assumes that the new molecules all appear in v=14, J=0 or 1 (after Hunter & Watson 1978). The third formation model is based on a treatment of H2 catalysis by Duley & Williams (1986) which predicts the molecules are ejected into the v=6 level. Naturally, there are significant differences in the predictions of the intensities for lines from v=6 between the first two and the third formation model. Depending on the formation temperature of the molecule, v=6 lines are predicted to be up to a factor $`3`$ times brighter in the third than in the first formation model. This suggests that the 1–0 S(1)/6–4 O(3) ratio may fall from $`10`$ to $`3`$ in Black & van Dishoeck’s (1987) Model 14, if it were to be modified so that H2 formation occured in v=6, rather than into a Boltzmann distribution.

We believe that formation pumping into v=6, or to a nearby level, provides the best explanation for the difference in the images between the 6–4 O(3) line and the lower excitation lines. However, it cannot simply be due to the additional component that formation pumping adds to the intensity of a line, for that would not result in a different morphology for the line. Either the rate at which formation is occurring must vary with position, or the formation spectrum itself is varying with position (for instance, due to a varying formation temperature), across the N–Bar of M 17.

In the Duley & Williams model, about 1.5 eV of the 4.5 eV bond energy released is assumed to remain with the H2 molecule, as it is ejected from the surface of the dust grain where it formed. This puts it into an excited rotational-vibrational state, namely $`\mathrm{v}6`$. The actual formation process depends intimately on the nature of the surface of the dust grains. In their 1986 paper Duley & Williams considered the surfaces to be highly defected silicates, to which the molecule is moderately strongly bound. A later paper (Duley & Williams 1993) further considers H2 formation on amorphous H2O ice (where it is weakly bound) and on aromatic carbon molecules (e.g. PAHs or HACs, where it is strongly bound). If the former occurs, they suggest that the H2 will be released in highly-excited states (i.e. the binding energy virtually all goes into internal energy in the H2 molecule, so it ends up, perhaps, in v=13). In the latter case, most of the binding energy would be distributed through the many degrees of freedom of the aromatic molecule, leaving the H2 in a low-excitation state. Clearly, the observations in M 17 support the silicate grains model for that source over the other two grain models, since the excess appears to be in v=6.

The PDR models discussed above were steady state models, where the dissociation front is stationary. In them, the photodissociation rate of hydrogen molecules by far–UV radiation is balanced by their formation rate on grains. This need not be the case. The timescale for the H2/H dissociation front to reach equilibrium, $`\mathrm{t}_{\mathrm{eq}}5\times 10^8/\mathrm{n}`$ years, can be long compared to variations in the radiation field (for instance, soon after a star switches on) (Hollenbach & Natta 1995). In M 17, the data constrain $`\mathrm{t}_{\mathrm{eq}}`$ to be less than $`10^4`$ years. The exposure of previously shielded H2 molecules to the far–UV radiation field, as the dissociation front moves further into a molecular cloud, increases the column of fluoresced gas. Thus, the H2 line intensities rise from their steady state values. The 1–0 S(1) line intensity can be elevated by an order of magnitude (Hollenbach & Natta 1995). Moreover, since this emission is dominated by pure fluorescence (collisions can only play a minor role in redistributing the level populations), the 1–0/2–1 S(1) ratio will be $`2`$. Any thermal contribution to the 1–0 S(1) line has been minimised.

During this time the rate of molecule destruction must exceed that of formation. Thus the relative proportion of the formation component to the fluorescent H2 line intensities will be reduced. We now consider whether this explanation can be applied towards explaining the 6–4 O(3) line emission from the N–Bar of M 17. In it, the line is brighter in the NE of the two emission filaments, the one furthest away (at least in projected distance) from the ionizing stars in the HII region. Our hypothesis is that a formation component contributes to the 6–4 O(3) line intensity. This would imply that, in the weaker of the two filaments, a smaller proportion of the emission has been produced by formation pumping, compared to that produced by pure fluorescence, than in the brighter of the filaments. Thus, this suggests that the SW of the two filaments contains a non-steady state photodissociation region, where the dissociation front is moving rapidly into the molecular cloud.

We can apply this interpretation to provide a rough estimate of the fraction of the 6–4 O(3) line intensity that derives directly from H2 formation in a steady state PDR. We assume that when the 1–0 S(1)/6–4 O(3) ratio is at it largest (i.e. 23), then the UV cascade dominates its excitation. We also assume that when it is at its smallest value (i.e. 5) this represents the steady state PDR. This then yields a formation component that is $`80`$% of the total 6–4 O(3) line intensity in the steady state, the remaining $`20`$% coming from the fluorescent cascade. This is obviously a crude estimate. Without further data on other high excitation H2 lines we cannot say more about the formation spectrum of H2.

## 5 SUMMARY

The molecular hydrogen line emission from the N– and SW–Bars of M 17 clearly is fluorescently excited, from inside a photdissociation region. The N–Bar has approximately constant density, whereas the SW–Bar contains a clumpy molecular medium. However, the different distribution of the emission from the v=6–4 O(3) line, compared to that of the lower excitation lines from the v=1 and 2 levels, indicates that another mechanism is significant for its excitation. We argue that this is formation pumping, newly formed molecules being released from grain surfaces in an excited state, at or near the v=6 level. This therefore resembles the suggestion put forward by Duley & Williams (1986) for H2 formation on the surface of highly defected silicate grains, where the molecules are released with $`1.5`$ eV of the 4.5 eV molecule binding energy. However, time-dependent PDR models are also needed in order to quantitatively model such emission. This is because the proportion of the H2 line emission that has been produced by formation pumping, in comparison to fluorescence, probably varies with position in the source.

If this hypothesis is correct it may be tested by imaging in other highly excited emission lines. For instance, the v=6–4 Q(1) line at 1.6011$`\mu `$m, the 5–3 O(3) line at 1.6131$`\mu `$m, the 7–5 S(1) line at 1.6201$`\mu `$m and the 7–5 Q(1) line at 1.7283$`\mu `$m, are all emitted from nearby levels and in the same atmospheric window. The v=6 and v=5 lines would likely look similar, but, depending on whether any H2 is pumped into v=7 on formation, the distribution of the latter two lines may look closer to that of the lower excitation v=1 lines. Imaging Fabry-Pérot etalons are suitable instruments for these observations, because of their high spectral resolution, needed to separate out other lines in the window, and because they enable the spatial morphology to be examined. Long-slit spectroscopy, again at high spectral resolution, across the emission bars, would permit further examination of the merits of the formation hypothesis. It would help determine which lines have a formation component to their emission fluxes. Such data are needed to provide constraints on what the H2 formation spectrum actually is.

## ACKNOWLEDGEMENTS

We are grateful to Michael Ashley, Lori Allen, Jung-Kyu Lee, Stuart Ryder, John Storey and Andrew Walsh, as well as to the staff of the Anglo Australian Telescope, for all their help and support during the conduct of the observations. We also thank Michael Merrill for providing the beautiful colour image of M 17, shown in Fig. 1. We also acknowledge the support of the Australian Research Council for the funding to undertake this research.


text:
# A probable planetary companion to HD 39091 from the Anglo-Australian Planet Search

## 1 Introduction

Since the discovery of the first extra-solar planet (Mayor & Queloz 1995), planetary detections have been dominated by northern hemisphere search programmes – most prolifically by the high precision velocity programmes at Lick (e.g., Butler et al. 1996) and Keck (e.g. Vogt et al. 2000) and lower precision programmes at OHP (Baranne et al. 1996), McDonald Observatory (Cochran et al. 1997), AFOE (Noyes et al. 1997), and programmes at La Silla (Kurster et al. 2000; Queloz et al. 2000). Of these programmes, which achieve precisions of $``$ 10 m s-1, only the latter has access to the sky south of $`20\mathrm{deg}`$, and these achieve precisions of $``$ 10 m s-1. In 1998, the Anglo-Australian Planet Search (AAPS) began in the southern hemisphere enabling all-sky coverage of the brightest stars at precisions reaching 3 m s-1. In this paper we present the second set of results from this programme. Two companion papers, Tinney et al. (2002a, 2002b), present results for a further two planets and an initial investigation of the Calcium H and K activity among some of the target stars of the AAPS.

Precision Doppler surveys have found all of the known extrasolar planets around solar-type stars. Discoveries have included: the first system of multiple planets orbiting a Sun-like star (Butler et al. 1999); the first planet seen in transit (Henry et al. 2000, Charbonneau et al. 2000); the first two sub-Saturn-mass planets (Marcy, Butler & Vogt 2000); and the AAPS’ discovery of the first planet in a circular orbit outside the 0.1 au tidal-circularisation radius (Butler et al. 2001). A number of major surprises have emerged from the sample of extrasolar planets.

* The sub-stellar companion mass function for F, G, and K dwarfs rises strongly below 10 MJUP, and shows no signs of flattening toward the detection limit near 1 MJUP. And surprisingly, brown dwarf companions to solar type stars are rare (Butler et al. 2000). This is in spite of a strong selection bias in the observations that makes brown dwarfs easier to detect than planets. It should be noted that the mass is only known within the uncertainty of the projection factor represented by the sin $`i`$ of the orbit. However, the sin $`i`$ statistic is between 1 and 0.5 for 87 per cent of randomly inclined orbits, so 87 per cent of reported $`M`$ sin $`i`$ values will be within a factor of two of the true masses.
* About $``$0.75% of nearby solar-type stars stars have been found to have planets orbiting in circularised orbits inward of 0.1 au (the 51 Peg-like “hot Jupiters”). A smaller fraction of stars is now being found to have “eps Ret -like” planets orbiting in circularised orbits at 1 au or so (Tinney et al. 2002). But the dominant class of extra solar planets, (found around some 7 per cent of target stars) show highly eccentric orbits within 3.5 au. None of these classes were predicted a priori by planetary formation theories. Such theories have received enormous impetus from these observations, leading to models for planet formation and evolution which now include the effects of dynamical friction, disk-planet and planet-planet interactions (e.g., Rasio & Ford 1996; Weidenschilling & Marzari 1996; Artymowicz 2000; Boss 2000).
* The majority of extrasolar planets that have been found so far occur around stars which are mildy metal-rich ($``$+0.2 dex with considerable scatter) compared to the Sun (Laughlin 2000; Santos, Israelian & Mayor 2002).

This paper reports the discovery of a new planet candidate from the AAPS. Section 2 describes this precision programme. The stellar properties and Keplerian orbital fits for the new planet candidate are presented in Section 3. Section 4 provides a discussion of the new object.

## 2 The Anglo-Australian planet search

The Anglo-Australian planet search (AAPS) is carried out on the 3.92m Anglo-Australian Telescope using the University College London Echelle Spectrograph (UCLES), operated in its 31 lines/mm mode together with an I2 absorption cell before the slit. UCLES used to be operated with an MIT/LL 2048$`\times `$4096 15$`\mu `$m pixel CCD. Since August 2001 the AAPS has changed to the AAO’s EEV 2048$`\times `$4096 13.5$`\mu `$m pixel CCD. This CCD provides 30% better quantum efficiency across the 5000–6200 Å region containing I2 absorption lines. Furthermore, the EEVs smaller pixels result in improved resolution and reduced charge diffusion. This CCD change has not required any significant changes of observational or reduction procedure and has resulted in significantly increased signal-to-noise. In the future we envisage further improvements in observing efficiency and S/N through the use of a (1) an exposure meter, (2) a motorised system for moving the I2 cell in and out of the beam and (3) a prime focus fibre feed for UCLES.

Doppler shifts are measured by observing through an I2 cell mounted behind the UCLES slit. The resulting superimposed iodine lines provide a fiducial wavelength scale against which to measure radial velocity shifts. The shapes of the iodine lines convey the PSF of the spectrometer for changes in optics and illumination on all time scales. We synthesize the echelle spectrum of each observation on a sub-pixel grid using a high resolution reference template, and fit for spectrograph characteristics (the wavelength scale, scattered light and the spectrograph PSF) and Doppler shift. This analysis obtains velocities from multiple epoch observations, which are measured against this reference template. This reference template is an observation at the highest available resolution (using a small 0.5 arcsec slit) and high S$`/`$N = 200 to 300, without the I2 cell present. Such measurements can only be efficently obtained in good seeing and take about 4 times as long to acquire as a standard epoch (I2 and a 1 arcsec slit) observation. We show in Fig. 1 the residuals about zero velocity for a sample of apparently stable stars reported by Butler et al. (2001), as well as the residuals about a Keplerian fit for several detected planets. We achieve 3 m s-1 precisions down to the V = 7.5 magnitude limit of the survey. The fundamental limit to the precision that can be achieved for our sample is set by a combination of S$`/`$N (which is dependent on good seeing and weather conditions), and the intrinsic velocity stability of our target stars often called “jitter”, induced by activity and/or star spots (Saar et al. 1998; Saar & Fischer 2000), rather than our observing technique. There is currently no way to tell whether a residual scatter of larger than 3 m s-1 is due to a small-amplitude planet (either short- or long-period – the detection of the latter is one of our primary goals, as these are Jupiter-like signals), or jitter induced by star spots and/or activity. Only observations over a long enough period to allow the search for long-term periodicites can reveal the presence of such relatively small-amplitude long-period signals such as Jupiter. It is therefore vital to monitor all our targets for the lifetime of the survey not just those that appear to be good planet candidates.

Our target sample includes 178 FGK stars with declinations below $`20\mathrm{deg}`$ and V$`<`$7.5, and a further 23 metal-rich stars with V$`<`$11.5. Where age/activity information is available from log $`R_{\mathrm{HK}}`$ indices (Henry et al. 1996; Tinney et al. 2002a) we require target stars to have log $`R_{\mathrm{HK}}`$ $`<`$ –4.5 corresponding to ages greater than 3 Gyr. Stars with known stellar companions within 2 arcsec are removed from the observing list, as it is operationally difficult to get an uncontaminated spectrum of a star with a nearby companion. Spectroscopic binaries discovered during the programme have also been removed and will be the subject of a separate paper (Blundell et al. 2002, in preparation). Otherwise there is no bias against observing multiple stars. The programme is not expected to have any bias against brown dwarf companions. The observing and data processing procedures follows that described in Butler et al. (1996,2001). Our first observing run was in January 1998, and the last run for which observations are reported here was in October 2001.

## 3 Stellar Characteristics and Orbital Solution for HD 39091

A total of 26 observations of HD 39091 (GJ 9189, HIP 26394, HR 2022) have been made by HIPPARCOS yielding a distance of 18.2 pc and a V magnitude of 5.65 (ESA 1997). The resulting absolute magnitude is $`M_\mathrm{V}`$ = 4.32. The star is photometrically stable within the HIPPARCOS measurement error, with a photometric scatter of 0.006 magnitudes. The Bright Star Catalogue (Hoffleit 1982) assigns a spectral type of G1V, in reasonable agreement with the HIPPARCOS spectral type of G3IV. The star is chromospherically inactive with log $`R_{\mathrm{HK}}`$ = –4.97 (Henry et al. 1996). Based on the spectroscopic analysis of Santos, Israelian & Mayor (2002), \[Fe/H\] = 0.09 and Mass = 1.1$`\pm `$0.02 M. Decin et al. (2000) find no evidence for a 60 $`\mu `$m excess in deep ISO data for HD 39091. It is interesting to note that Queloz et al. (2000) report HD 39091 as a non-variable star, however, the Queloz et al. survey would not yet have been sensitive to the long period of HD 39091.

The 26 Doppler velocity measurements of HD 39091, obtained between November 1998 and October 2001, are listed in Table 1 and shown graphically in Fig. 2. The third column labelled error is the velocity uncertainty produced by our least-squares fitting. This fit simultaneously determines the Doppler shift and the spectrograph point-spread function for each observation made through the iodine cell, given an iodine absorption spectrum and an iodine-free template of the object (Butler et al. 1996). This uncertainty includes the effects of photon-counting uncertainties, residual errors in the spectrograph PSF model, and variation in the underlying spectrum between the template and iodine epochs. All velocities are measured relative to the zero-point defined by the template observation. Only observations where the uncertainty is less than twice the median uncertainty are listed. The best-fit Keplerian curve yields an orbital period of 2083 d, a velocity amplitude of 196 m s-1, and an eccentricity of 0.62. The minimum (M sin $`i`$) mass of the planet is 10.3 MJUP, and the semi-major axis is 3.34 au. The RMS to the Keplerian fit is 6.6 m s-1, yielding a reduced chi-squared of 1.1. The properties of the extra-solar planet in orbit around HD 39091 are summarised in Table 2. Whilst we do not yet have a full orbit for HD 39091, we do clearly have a 2nd inflection, which means that the orbit is starting to become clear. However, the parameters of the orbit are relatively poorly constrained and any follow-up observations should take this into account. We are announcing this object at this somewhat early stage in order that Adaptive Optics and astrometric follow-up on this relative long-period object can be started as soon as possible.

We investigated the HIPPARCOS dataset for HD 39091 using the set of tests devised by Pourbaix & Arenou (2001). These tests combine radial velocity measurements with those from HIPPARCOS to derive upper mass limits for companions. HD 39091 fails all the Pourbaix & Arenou tests apart from one. The Thiele-Innes orbit does improve the fit. However, the derived elements are not significantly different from zero and thus no constraint can be placed on the orbit. The reason for this is the long period. It exceeds the mission duration by a factor two so the size of the orbit is not constrained at all. Neither does HIPPARCOS give a lower limit for the inclination. Such a limit is based on the amplitude (the astrometric semi-major axis of the orbit) corresponding to the spectroscopic period. Once again, the method can not be applied to HD 39091 since the period exceeds the mission duration. Hence, unfortunately, HIPPARCOS cannot provide any additional information.

## 4 Discussion

It is interesting to note the mass of 10.3 MJUP places HD 39091’s companion rather close to the canonical deuterium fusion limit of around 13 MJUP. Most definitions of a planet posit that they should not be undergoing any nuclear fusion. Unfortunately distinguishing the presence of deuterium and thus the lack of fusion is still beyond our current observational capabilities. Once we have (1) the capability to make high resolution and high signal-to-noise infrared spectroscopic measurements of such faint close-in companions and (2) sufficiently high quality molecular line lists of HDO and H2O, HD 39091’s mass and long orbit make it a prime target.

The AAT, Keck and Lick precision velocity surveys are currently surveying a total of 1200 stars. All three of these programmes use the Iodine cell technique, and all three have demonstrated long-term precision of 3m s-1 (Marcy & Butler 1998; Vogt et al. 2000; Butler et al. 2001). The ellipticities for the substellar candidates from these surveys are shown in Fig. 3. The extra-solar planets show a wide range of eccentricity similar to the eccentricity distribution of binary systems with stellar companions (e.g., Duquennoy & Mayor 1991). This is not understood in terms of a global planetary formation model that also allows for the almost circular orbits in our Solar System. Given that we have only just started to explore the potential parameter space for extra-solar planets, history should remind us not to envisage the Solar System as a special case. It is important to note that our studies have not yet run for long enough to be sensitive to a true Solar System analogue.

The star HD 39091 is established as metal-rich star relative to the Sun which is itself metal-rich compared to the average of stars in the solar neighbourhood. The companion to HD 39091 reported here thus further confirms the statistical findings of Santos et al. (2002) that stars with planets are metal rich compared to stars without planet detections. There are two general classes of explanation for this: intrinsic metallicity bias and accretion of metal-rich material. After a variety of different spectral analyses and claims a straightforward explanation remains: the higher metallicity of planet-harbouring stars arises because high metallicity environments have a higher probability of planet formation.

## 5 Conclusions

We present data showing evidence for an eccentric-orbit extra-solar planet around the metal-rich star HD 39091. The detection of this long period object gives added impetus for the continuation of these searches to longer periods. We now must endeavour to continue to improve the precision and stability of the AAPS to be sensitive to the 10+ year periods where analogues of the gas giants in our own Solar System may become detectable around other stars.

## Acknowledgments

The Anglo-Australian Planet Search team would like to thank the support of the Director of the AAO, Dr Brian Boyle, and the superb technical support which has been received throughout the programme from AAT staff – in particular E.Penny, R.Patterson, D.Stafford, F.Freeman, S.Lee, J.Pogson and G.Schafer. We gratefully acknowledge the UK and Australian government support of the Anglo-Australian Telescope through their PPARC and DETYA funding (HRAJ, CGT, AJP); NASA grant NAG5-8299 & NSF grant AST95-20443 (GWM); NSF grant AST-9988087 (RPB); and Sun Microsystems. DP is a post-doctoral researcher at the National Fund for Scientific Research, Belgium. This research made use of the Simbad database, operated at CDS, Strasbourg, France.