text:
# Peptide size dependent active transport in the proteasome

## Abstract

We investigate the transport of proteins inside the proteasome and propose an active transport mechanism based on a spatially asymmetric interaction potential of peptide chains. The transport is driven by fluctuations which are always present in such systems. We compute the peptide-size dependent transport rate which is essential for the functioning of the proteasome. In agreement with recent experiments, varying temperature changes the transport mechanism qualitatively.

Eukaryotic mammalian proteasomes are fascinating molecular machines, intended to digest tagged with ubiquitin proteins 2001\_Goldberg\_sa . The proteins are cleaved into peptides which can be further used for the synthesis of new proteins. The proteasome’s function has been directly linked to a wide range of diseases, such as cancer, neurological diseases, and diseases of the immune defense system. Starting from the discovery of their activity in the late 1970’s, presently proteasomes are in the focus of current molecular biology, see HiltWolf:2001CiechanoverMasucci:2003  and refs. therein. By now there are several models to explain the proteasome’s activity 2003\_Schliwa\_nature , cleavage mechanisms 2000\_Altuvia\_jmb ; 1999\_Holzhuetter\_jmb ; 2002\_Kesmir  and to predict the cleavage results 1999\_Kisselev\_jbc , however, basic principles of the proteasome operation mechanisms are still poorly understood, mainly due to the lack of experimental results.

In this Letter, we focus on understanding of protein translocation inside the proteasome, leaving the mechanisms of cleavage, targeting, etc. beyond the scope. The main question is: given proteasomes are highly complex pipe-like structures of tenthousands atoms of almost perfect left-right symmetry with respect to the axis of the pipe 1998\_Walz , thus, operating bi-directional. To be cleaved, a protein has to enter the proteasome at one side, pass the active sites (where the cleavage occurs), which makes about 1/3 of the total length of the pipe. Then the cleaved peptides have to pass all the way through the pipe to finally reappear at the other side of the proteasome. Although there is a number of examples where protein transport in cells occurs due to diffusion, i.e. Brownian motion 2001\_Reits\_nature\_cell , diffusive transport may be excluded as the main mechanism for translocation in the proteasome because of the enormous cargo 2000\_Elston\_siam . Also other proposed transport mechanisms, such as the power stroke model of protein translocation do not suffice to explain translocation 2002\_Elston\_bio\_journal . Therefore, the question arises how the protein’s motion inside the proteasome is driven?

Since proteasomes are large multi-subunit structures consisting of proteins, the mechanism of the protein transport is directly related to protein-protein interaction. In 2001\_Brokaw\_bio\_journal  it has been noticed that if attachment and detachment rates are specified asymmetrically, the protein-protein binding interaction acts as a ratchet. Following this argumentation, active protein transport, based on the mechanism of a molecular ratchet, has been studied for transport through membranes 1990\_Vale1994\_Astumian\_prl2002\_Reimann\_pr  and has been also discussed as a mechanism for cytosolic destruction by the proteasome 1999\_Matlack\_cell . Moreover, maximum likelihood tests have shown that other models, e.g. the power stroke model of protein translocation, do not lead to better agreement with the experiment than the Brownian ratchet model 2002\_Elston\_bio\_journal . Noteworthy, in these ratchet effects transport is possible only in a certain temperature interval, and stochasticity, intrinsically present due to fluctuations in any biochemical reaction 1981\_van\_Kampen\_book ; 1999\_Barkai\_nature , provides the driving mechanism.

Here we propose a model for active protein translocation in the proteasome to explain the peptide size dependence of the transport velocity as well as its temperature dependence which possibly explains the mechanism of temperature reaction or heat shock response, regulating the proteasome activity in the case of some diseases 2000\_Kuckelkorn ; 2000\_Wyttenbach\_pnas ; 2001\_Pritts ; 2001\_Stepanenko\_mol\_biology . The results describe a system size ratchet effect that is related to similar effects which have been described recently for other noise-induced phenomena 2002\_Pikovsky\_prl . It has been shown in vitro (e.g. HiltWolf:2001CiechanoverMasucci:2003 ) that already the 20S proteasome reveals full functionality, thus, for studying the transport we disregard the enzymatic action of the proteasome caps.

Model. Following, e.g., 2001\_Brokaw\_bio\_journal  we assume that the protein-proteasome interaction is characterized by a spatially periodic asymmetric potential, motivated by asymmetric C-N binding between amino acids that enter the proteasome always in $`NC`$ direction, i.e., with $`N`$ head on 1999\_Holzhuetter\_jmb  (see Fig. 1). The folding structure of the protein is not relevant here since due to the action of the regulatory complexes the protein enters the proteasome unfold 1997\_Kopito\_cell . Moreover, even fold proteins tends to form periodic structures to maximize their amphiphilicity 1984\_Eisenberg\_pnas .

As an abstraction, here we assume that the protein-proteasome interaction potential $`U(x)`$ is periodic with the period $`L`$. In reality there is a basic asymmetry, namely the $`CN`$ asymmetry of the protein (or peptide) backbone, that is superposed by a non-periodic (in our sense irregular) part that is attributed to the amino-acid-specific residues asymmetry . The basic structure of a peptide and its model as an asymmetric ratchet potential is sketched in Fig. 1. The commensurability of the protein and the proteasome is supported by the fact that both biological macro-molecules consist of the same basic structures, namely chains of amino-acid. That commensurability has been claimed before, and was the basis of automaton-like models of proteasome digestion processes 1999\_Holzhuetter\_jmb .

The detailed form of the asymmetric periodic interaction potential is of less importance for this qualitative study. Here, we assume a saw-tooth potential as drawn in Fig. 1. The angles are smoothed, i.e., $`dU/dx`$ can be computed in each point (for details see 2001\_Landa\_Zaikin\_csf ).

The proteasome acts upon the protein by a certain number of equidistant interaction centers. The dynamics of the protein inside the proteasome is, hence, governed by $`N`$ interactions centers, where $`N`$ is the number of protein elements (amino acids or multiplicatives of it). There appear the following forces: potential force (protein-proteasome interaction) $`NU(x)/x`$, fluctuations with collective $`NF(t)`$ and individual components $`f_1(t)+\mathrm{}+f_N(t)`$, and protein friction forces $`N\beta \dot{x}`$ 2001\_Brokaw\_bio\_journal , where $`x`$ is the coordinate of the protein with respect to the proteasome and $`\beta `$ is the coefficient of friction. Due to the small mass of all protein particles, moving in the liquid cytosol, the motion occurs in the overdamped realm 2003\_Bier\_prl , hence we neglect inertia forces. Note that transport is possible only in the case of nonequilibrium fluctuations. In the simplified case, when fluctuations can be represented by a sum of a collective periodic force and individual for every protein residue thermal noise, the model is analytically tractable, predicting the velocity dependence on the peptide size. Numerically, we investigate also different kinds of fluctuations in order to be closer to reality. Normalizing all forces by friction and taking $`\beta =1`$, the translocation of a protein in the proteasome is governed then by

$`{\displaystyle \frac{x}{t}}={\displaystyle \frac{U(x)}{x}}+F(t)+{\displaystyle \frac{1}{N}}(f_1(t)+\mathrm{}+f_N(t)).`$ (1)
When the protein chain enters the proteasome, the number of interaction centers $`N`$ is increased (Fig. 2).

After cleavage, $`N`$ is abruptly decreased for the cleavage products. We believe that this crucially changes the transport velocity, providing the cleavage product leaving the proteasome. Eq.(1) models the motion for both the initial protein and the peptides after cleavage. We consider three cases of fluctuations, regulating their motion:

Case 1: We assume collective oscillations of the peptide elements 1997\_Julicher\_prl ; 2000\_Tama\_proteins , e.g., $`F(t)=A\mathrm{cos}(\omega t)`$, where $`A`$ and $`\omega `$ stand for the amplitude and frequency of these oscillations. Additionally, each interaction center undergoes local thermal fluctuations, represented by mutually uncorrelated white noise of intensity $`\sigma ^2`$: $`f_i(t)=\xi _i(t)`$, where $`\xi _i(t)\xi _j(t^{})=\sigma ^2\delta (tt^{})\delta _{ij}`$. In this case the stochastic term in Eq. (1) is white noise of intensity $`\sigma ^2/N`$. The Fokker-Planck equation for the peptide coordinate probability distribution $`w(x,t)`$ associated with Eq. (1) is

$$\frac{w}{t}=\frac{}{x}\left[\left(F(t)\frac{U}{x}\right)w(x,t)\right]+\frac{\sigma ^2}{2N}\frac{^2w(x,t)}{x^2},$$

which may be solved in quasi-stationary adiabatic approximation $`w/t=0`$ 1990\_Jung\_pra . We obtain

$$\frac{\sigma ^2}{2N}\frac{w(x,F)}{x}\left(F\frac{U}{x}\right)w(x,F)=G(F),$$
(2)
where $`G(F)`$ is the probability flux. For any periodic potential $`U(x)`$ the quasi-stationary solution of Eq. (2) is

$$\begin{array}{cc}\hfill w(x,t)=& \left[C(F)\frac{2G(F)}{\sigma ^2/N}_0^x\mathrm{exp}\left(\frac{U(x^{})Fx^{}}{\sigma ^2/2N}\right)𝑑x^{}\right]\hfill \\ & \mathrm{exp}\left(\frac{U(x)Fx}{\sigma ^2/2N}\right),\hfill \end{array}$$

where $`C(F(t))`$ and $`G(F(t))`$ are unknown functions of $`t`$. Using the periodicity condition $`w(0,t)=w(L,t)`$ and the normalization of $`w(x,t)`$ we get $`G(F)`$. If the amplitude $`A`$ meets the condition $`LA\sigma ^2/N,`$ one can expand $`G(F)`$ and obtain

$$G(F)G_{01}F+G_{02}F^2$$
(3)
with the expansion coefficients $`G_{01}=L/(I_{10}I_{20})`$,

$$\begin{array}{cc}\hfill G_{02}& =G_{01}\left(\frac{I_{11}}{I_{10}}\frac{I_{21}}{I_{20}}\frac{NL}{\sigma ^2}\left(1\frac{2I_{30}}{I_{10}I_{20}}\right)\right)\hfill \\ \hfill I_{10}& =_0^Le^{U^{}(x)}𝑑x,I_{20}=_0^Le^{\left(U^{}(x)\right)}𝑑x\hfill \\ \hfill I_{11}& =\frac{2N}{\sigma ^2}_0^Lxe^{U^{}(x)}𝑑x,I_{21}=\frac{2N}{\sigma ^2}_0^Lxe^{\left(U^{}(x)\right)}𝑑x\hfill \\ \hfill I_{30}& =_0^L_0^xe^{\left(U^{}(x^{})U^{}(x)\right)}𝑑x^{}𝑑x,U^{}(x)=\frac{2NU(x)}{\sigma ^2}.\hfill \end{array}$$

Substituting Eq. (3) into $`\overline{\dot{x}}=_0^L\overline{G(x,t)}𝑑x`$, where $`\overline{()}`$ denotes time averaging, we obtain the average protein transport velocity, as a function of the noise intensity $`\sigma ^2`$ and the peptide size $`N`$

$$\overline{\dot{x}}\frac{L^2A^2}{2I_{10}I_{20}}\left[\frac{I_{11}}{I_{10}}\frac{I_{21}}{I_{20}}\frac{NL}{\sigma ^2}\left(1\frac{2I_{30}}{I_{10}I_{20}}\right)\right].$$
(4)
The dependence of the transport velocity on the noise intensity, i.e. temperature, is shown in Fig. 3 (left) for different system size $`N`$. This plot shows that protein transport is possible only in a certain temperature interval. Absence of noise leads to vanishing transport as well as too high temperature also decreases the velocity. For given temperature, the peptide size $`N`$ strongly influences the transport velocity, i.e., with increasing $`N`$ the dependence is shifted to the right and becomes more wide (curves 1-3). This dependence differs qualitatively for small and large temperature. For rather small temperature the velocity decays with increasing $`N`$ (curve 1 in Fig. 3 (right)). Consequently, when a protein enters the proteasome, it will initially move with rather large velocity. Then, due to the increase of $`N`$ the velocity decreases, and the motion of the protein virtually terminates. On this stage the probability of cleavage is much higher than that of translocation, thus, with large probability the protein will be cleaved there. After cleavage, $`N`$ is significantly decreased (Fig. 2 (right)), and following the velocity dependence, both fragments move rapidly to leave the proteasome. Noteworthy, in agreement with our model, the blocking of uncleaved proteins inside the proteasome has been observed also experimentally 2001\_Reid\_pnas .

Surprisingly, for large temperature the velocity depends qualitatively different on the peptide size (curves 2,3 in Fig. 3 (right)), i.e., there is an optimal peptide size that corresponds to maximal velocity. In this case the protein moves accelerated while entering the proteasome until the optimal velocity is achieved. Then the velocity decreases, and again it virtually stucks. After cleavage the behavior is quite different: if the cleaved peptide is close to the optimal size, it leaves the proteasome rapidly. If the cleaved fragment is too small or too large, it moves with small velocity, thus blocks the proteasome leading to low efficiency or malfunction of the proteasome. Note that for the protein that enters the proteasome the initial small $`N`$ increases, and the velocity will increase, whereas for the cleaved fragment of smaller $`N`$ the velocity remains constant. There is experimental evidence 2001\_Reid\_pnas  that certain proteins cannot be cleaved by the proteasome, i.e., they are inert with respect to its active sites. For such proteins we can determine the penetration depth, defined as the depth $`N`$ at which the propagation velocity falls below a predefined very small threshold. The inset of Fig. 3 (right) shows that the penetration depth increases linearly with temperature. So far, there are no experimental data to verify this prediction.

Consequently, small temperature leads to preferably small size of the cleaved fragments, whereas larger temperature leads to domination of some specified peptide size in the proteasome output. As a hypothesis, we believe that the qualitatively different proteasome behavior with respect temperature may be responsible for temperature reaction or heat shock response 2000\_Kuckelkorn ; 2001\_Stepanenko\_mol\_biology ; 2000\_Wyttenbach\_pnas , providing immune defense in the case of some diseases. For large temperature (Fig. 3 (right) curve 2,3) the protein is accelerated while entering the proteasome, hence, leading to more efficient operation of the proteasome.

To confirm these results we have performed computer simulations 1992\_Kloeden\_book  whose results agree with our analytical findings (Fig. 4).

Some quantitative deviations result from taking into account only a limited number of terms in the expansion Eq.(3).

Case 2: Since collective fluctuations of the protein or peptide elements, represented by a periodic in time force, can be used only under certain assumption, next we consider the system in the absence of deterministic periodic forces, but with colored noise $`F(t)=\zeta (t)`$, where $`\zeta (t)`$ is harmonic noise with some dominating frequency (due to the molecular spectrum and oscillations of the peptide chain 2000\_Tama\_proteins ; 1997\_Julicher\_prl ). We obtain this noise as the solution of

$$\ddot{\zeta }+2\delta \dot{\zeta }+\omega _0^2\zeta =\xi ^{},$$
(5)
where $`\xi ^{}`$ is Gaussian noise $`\xi ^{}(t)\xi ^{}(t^{})=\sigma _a^2\delta (tt^{})`$, and parameters $`\delta =0.01`$, $`\omega _0=0.1`$, $`\sigma _a^2=10^4`$ regulating the width, the dominating frequency of the power spectrum, and the intensity of $`\zeta (t)`$. The results show qualitatively the same behavior as in the case of the periodic force (Fig. 4 (b,e)). Hence, if exclusively noisy fluctuations are present in the system, the generic behavior described above, persists and leads to the same temperature and peptide size dependences.

Case 3: Next we consider the case of no collective fluctuations, $`F(t)=0`$, i.e., there is no eigen collective mode in the peptide oscillations. The local noise of each interaction center is presented now by the sum of uncorrelated colored and white noise components: $`f_i(t)=\zeta _i(t)+\xi _i(t)`$. The presence of the colored noise is necessary to generate nonequilibrium fluctuations, that is a necessary condition for directed transport. In this case the motion of the protein or peptide is described by

$`\dot{x}={\displaystyle \frac{U(x)}{x}}+{\displaystyle \frac{1}{N}}(\zeta _1+\mathrm{}+\zeta _N+\xi _1+\mathrm{}\xi _N),`$ (6)
with $`\zeta _i(t)`$ and $`\xi _i(t)`$ as above. In general, even in this case the behavior is similar, but increase of the peptide size leads not only to a shift of the velocity dependence towards large noise intensity, but also to a decrease of its maximal values (Fig. 4 (c)). This is due to the fact that nonequilibrium fluctuations originating from different interaction centers act occasionally in converse direction. For small temperature (Fig. 4 (f), curve 1) the transport decreases with increasing peptide size $`N`$, whereas for large temperature the rate depends nonmonotonously on $`N`$ (curve 2). The optimal velocity in the latter case is smaller, as for the case of collective fluctuations.

In summary, we investigated a model for active protein transport inside the proteasome analytically and numerically. The size of the peptide as well as temperature influence the transport rate significantly in agreement with experimental results. Following our predictions, protein transport occurs only in the certain interval of temperature. The size dependence of the transport velocity leads to a preferred fragment (peptid) size since for such sizes the probability of cleavage exceeds the probability of further translocation significantly. Under certain conditions, uncleaved proteins may get stuck inside the proteasome, thus leading to decrease of its efficiency which is in good agreement with experimental observations 2001\_Reid\_pnas . Calculating the average velocity, we have considered only the periodic constituent of the protein-proteasome potential. Slight inhomogeneity (nonperiodicity) of this potential will lead to fluctuations of the velocity around its average leaving the qualitative behavior of dependencies on the peptide size intact.

Varying temperature, the model predicts qualitatively different transport regimes, thus, possibly explaining the mechanism of temperature reaction or heat shock response as it is observed for certain diseases. For larger temperature, the protein is accelerated while entering the proteasome, moreover, there is a preferred fragment size, that should dominate the proteasome output. We emphasize that our model is based only on two assumptions: spatial asymmetry of the protein chain and presence of nonequilibrium fluctuations, which seem to be certainly fulfilled in real proteasome machines.

We thank C. Frömmel, H. Holzhütter, and R. Preißner for helpful discussion.


format: text

text:
# The physical and chemical structure of hot molecular cores

## 1 Introduction

Hot cores are compact (diameters $`0.1`$ pc), dense (n$`10^7`$ cm-3), hot (T$`100`$ K), and dark ($`A_v100`$ mag) molecular cloud cores. At least some of them are believed to be precursors to ultracompact HII regions and in the early stage of forming massive stars because luminous infrared sources, massive bipolar outflows and/or inflow motion are observed associated with them (e.g., Henning et al. 2000; Kurtz et al. 2000; Churchwell 2002). They are observationally characterized by anomalously large abundances of some molecular species including hydrogenated, saturated, small molecules such as H2O and H2S, and large complex molecules such as HCOOCH3 and (CH3)2O, in contrast to colder dark clouds. It is thought that these molecules are produced under conditions specific to star formation: the freeze-out of gas-phase species onto dust grains and the grain surface reactions in the cold collapsing pre-protostellar phase, and the evaporation of the grain mantle material and the subsequent gas-phase reactions in the hot protostellar phase (e.g., Millar 1993, 1997). Therefore, the observed large column densities of some molecules in hot cores are also expected to be an indicator of luminous high mass star formation which yields a large mass of high temperature gas. Previously, chemical models of hot cores have been studied together with the pre-stellar phase or with some initial conditions of evaporated grain mantle species, and have been successful in reproducing their observational properties (e.g., Brown et al. 1988; Millar et al. 1991; Charnley et al. 1992; Caselli et al. 1993; Charnley 1997; Millar et al. 1997; Viti & Williams 1999; Rodgers & Charnley 2001; Doty et al. 2002).

Physical properties of hot cores have been investigated on the basis of observational studies of molecular line and dust continuum emission. Some physical models together with radiative transfer calculation and/or chemical calculation have succeeded in reproducing the observations of spectral energy distributions and/or molecular column densities (e.g., Millar et al. 1997; Kaufman et al. 1998; Osorio et al. 1999; Hatchell et al. 2000; van der Tak et al. 2000a; Doty et al. 2002). On the other hand, the dynamical processes and the formation mechanisms of massive stars are not yet established because of the difficulty of treating in a self-consistent manner young stellar evolution, the accretion flow and radiation pressure. Non-spherical accretion to protostellar cores, for example, through accretion disks, and/or coalescence of lower mass stars may be involved in massive star formation, but there are, as yet, no accepted theories (e.g., Bonnell et al. 1998; Stahler et al. 2000; Yorke & Sonnhalter 2002). The nature and timescale of the dissipation of the surrounding envelopes of young massive stars, for example, by powerful outflows and ionizing strong winds which stop the accretion, are also uncertain but important in determining the mass of forming stars in the case of the isolated star formation (e.g., Shu et al. 1987; Nakano et al. 1995).

In this paper, we have modeled the density and temperature profiles of hot cores, making use of radiative transfer calculations, and then investigated the hot core chemistry, taking into account the temperature dependent grain mantle evaporation. In the following section, we present the physical model of the density profiles constrained by observational results and the temperature profiles calculated by solving the radiative transfer equation. In Sect. 3 we introduce the chemical model along with the temperature dependent grain mantle evaporation process. These models are applied to the hot core G34.3+0.15, and the results of chemical calculation are described and compared with the observations in Sect. 4. Also, we investigate the dependence of the chemical structure of hot cores on (1) the trapping of mantle molecules in water ice and (2) the density profile generally to discuss the chemical and physical condition of massive-star-forming clumps in Sect. 5. Finally, the results are summarized in Sect. 6.

## 2 Physical model

In this paper we consider a spherically symmetric molecular clump illuminated by a central star with luminosity $`L_{}`$. The density profiles modeled are constrained by observational results (Sect. 2.1), and the temperature profiles are calculated self-consistently under the assumption of radiative equilibrium by solving the radiative transfer equation (Sect. 2.2).

### 2.1 Density profiles

As we will see in Sect. 5, the density profile of a clump is one of the crucial properties which affects its chemical structure. However, there is no established theoretical model because the dynamical process of massive star formation is still unknown due to the difficulty in comprehensive treatments of young stellar evolution and accreting flow as well as ionizing radiation, radiation pressure, stellar wind, bipolar outflow, and so forth (e.g., Stahler et al. 2000; Yorke & Sonnhalter 2002). Thus, in this paper we model the density profiles of clumps simply based on the submillimetre observations by Hatchell et al. (2000), avoiding uncertainties in theoretical models.

Hatchell et al. (2000) investigated spectral energy distributions (SEDs) and radial profiles of 450 and 850 $`\mu `$m continuum dust emission of five ultracompact HII regions to conclude that chemically rich hot molecular clumps can be modeled by a combination of compact high density cores plus $`r^{1.5}`$ density profile envelopes. Also, van der Tak et al. (2000a) studied 14 envelopes around young massive stars by means of CO, CS, and H2CO line emission maps, and SEDs and radial profiles of submm dust emission to find that the best-fit values of the power law indices for density profiles of the envelopes range from $`1.0`$ to $`1.5`$. Theoretically, an isothermal self-gravitating gas is known to have $`r^2`$ density profile via the balance between self-gravity and pressure gradient forces (e.g., Ebert 1955; Bonnor 1956; Larson 1969; Shu 1977). If the nonthermal velocity dispersion, which has been observed in line widths towards various sizes of interstellar clouds and known to have a scaling-law of $`\delta vr^{0.5}`$ (e.g., Larson 1981; Falgarone et al. 1992; Caselli & Myers 1995), affects the pressure gradient force, it may flatten the radial density profiles of clumps like the logotrope model (e.g., McLaughlin & Pudritz 1996).

Taking into account of these observational and theoretical studies, we use the following spherically symmetric hydrogen number density distribution in this paper:

$$n(r)=\{\begin{array}{cc}n_0/(1+r/r_\mathrm{c})^2+n_1/(1+r/r_{\mathrm{tr}})^{1.5}\hfill & \mathrm{if}r_\mathrm{c}<r_{\mathrm{tr}},\hfill \\ n_1/(1+r/r_\mathrm{c})^{1.5}\hfill & \mathrm{otherwise},\hfill \end{array}$$
(1)
where the radius of $`r_{\mathrm{tr}}=0.05`$ pc is adopted for the transition between isothermal and nonthermal regions, according to the observations of molecular line widths (e.g., Fuller & Myers 1992; Myers & Fuller 1992). The hydrogen column density of the clumps, $`N_\mathrm{H}2_R_{}^{r_{\mathrm{out}}}n(r)𝑑r`$, is fixed as $`N_\mathrm{H}=10^{25}`$ cm-2, based on the observations of the chemically rich hot molecular clumps by Hatchell et al. (2000). The inner radius is set at the stellar surface, $`R_{}`$, and the outer radius is set at $`r_{\mathrm{out}}=1.0`$ pc. The central densities $`n_0`$ and $`n_1`$ are determined to satisfy the relation of $`n_0/(1+r_{\mathrm{tr}}/r_\mathrm{c})^2=n_1/(1+r_{\mathrm{tr}}/r_{\mathrm{tr}})^{1.5}`$ and the column density of $`N_\mathrm{H}=10^{25}`$ cm-2.

### 2.2 Temperature profiles

The temperature profiles are calculated self-consistently, assuming that the heating source of the clump is a central star with luminosity $`L_{}`$, under the condition of local radiative equilibrium, that is, the balance between absorption and reemission of radiation by dust grains at each point in the clump:

$$4\pi _0^{\mathrm{}}𝑑\nu \kappa _\nu B_\nu [T(r)]=_0^{\mathrm{}}𝑑\nu \kappa _\nu 𝑑\mu 𝑑\varphi I_\nu (r;\mu ,\varphi ),$$
(2)
where $`T`$, $`\kappa _\nu `$, and $`B_\nu (T)`$ represent the temperature, the monochromatic opacity, and the Planck function for blackbody radiation at a frequency $`\nu `$, respectively. Local thermodynamic equilibrium, $`\eta _\nu =\kappa _\nu B_\nu (T)`$, is adopted, where $`\eta _\nu `$ is the monochromatic emissivity. The specific intensity $`I_\nu `$ is calculated by solving the radiative transfer equation,

$$I_\nu (r;\mu ,\varphi )=_0^s\kappa _\nu (r^{})\rho (r^{})B_\nu [T(r^{})]e^{\tau _\nu (r^{})}𝑑s^{},$$
(3)
where $`\tau _\nu (r^{})`$ is the specific optical depth from $`r^{}`$ to $`r`$. The short characteristic method in spherical coordinates (Dullemond & Turolla 2000) is used to compute the integration. The radiative transfer equation is solved only in the dusty region (where the dust grains are not thermally destroyed) in this calculation. The inner radius for the numerical calculation is determined as $`r_{\mathrm{in}}=(T_{\mathrm{in}}/T_{})^2R_{}`$, where $`T_{}`$ and $`R_{}`$ are the temperature and radius of the central star, assuming that the radial profile of the temperature in the dust-free region ($`r<r_{\mathrm{in}}`$) is given by $`T(r)r^{1/2}`$, which represents the profile in the optically thin region with gray opacity heated by a point source. For the dust opacity model, we adopt the frequency dependent absorption coefficient $`\kappa _\nu `$ of Adams & Shu (1985, 1986), which approximates the interstellar dust grains model by Draine & Lee (1984). The dust destruction temperature is set as $`T_{\mathrm{in}}=`$ 2300K in this model (see Nomura 2002 for the details of radiative transfer calculation). In this paper the gas temperature is assumed to be the same as the dust temperature, based on the result of our calculation that the difference between the gas and the dust temperatures is small, except for the very inner region of the clump, where the gas temperature is $`10^3`$ K due to photoelectric heating from small grains induced by the ultraviolet radiation from the central star (cf. Doty & Neufeld 1997).

## 3 Chemical model

It is generally believed that the relatively large abundances of hydrogenated, saturated, small molecular species, such as H2O and H2S, and large complex molecules, such as HCOOCH3 and (CH3)2O, observed in hot cores are related to the evolution of young stellar objects and the evaporation of grain mantle material. In the cold, pre-stellar phase of the dense molecular clumps the accretion timescale of gas-phase atoms and molecules onto grain surface, $`3\times 10^9/n`$ \[cm-3\] yr at $`T=`$10 K, is short enough since the density, $`n`$, is high via self-gravitating contraction of the clumps. Thus, the accretion of gas-phase species onto grain surfaces and subsequent surface chemistry produce grain mantles rich in saturated species. Once massive star formation occurs in the clumps, the molecules are evaporated from the grain surface and the subsequent gas-phase reactions produce the large complex molecules seen in the hot gas (see e.g., reviews by Millar 1993, 1997 and references therein), although there is some evidence that even the most complex of molecules may have been formed on grains. On the other hand, the recent detailed analyses of observations of ice absorption features in infrared spectra towards massive young stellar objects have suggested that icy grain mantles in envelopes surrounding the central stars have different composition at different temperatures along the line-of-sight, that is, each icy mantle molecule sublimates from dust grains at a specific temperature (e.g., Ehrenfreund et al. 1998). In this paper, we set the molecular species evaporated from grain mantles as an initial condition, and then calculate the subsequent time-dependent gas-phase chemical reactions, taking into account the different binding energies of the mantle molecules.

The initial fractional abundances of species with respect to total hydrogen nuclei listed in Table 1 are used in this paper. The initial gas is assumed to be mostly molecular and slightly ionized by cosmic rays. The initial mantle composition is determined so as to be consistent with the recent observations of ice absorption features in infrared spectra towards massive young stellar objects (Gerakines et al. 1999; Keane et al. 2001) and other hot core models (Charnley et al. 1992; Millar et al. 1997; Rodgers & Charnley 2001; Doty et al. 2002), and fixed so as to best fit the results of model calculations to observations towards G34.3+0.15 (see Sect. 4). Also, the total abundances of carbon, oxygen, and nitrogen are set to agree with the mean interstellar values of C/H$`=1.4\times 10^4`$, O/H$`=3.19\times 10^4`$, and N/H$`=7.5\times 10^5`$ (Cardelli et al. 1996; Meyer et al. 1997; Meyer et al. 1998).

Each surface species, $`i`$, is expected to evaporate thermally with the timescale of

$$\tau _{\mathrm{evap},i}=\nu _{0,i}^1\mathrm{exp}(E_{b,i}/kT),$$
(4)
where $`k`$ and $`T`$ are the Boltzmann constant and the temperature, respectively. The vibrational frequency $`\nu _{0,i}`$ is given by $`\nu _{0,i}=(2n_sE_{b,i}/\pi ^2m_i)^{1/2}`$, where $`n_s(7\mathrm{\AA })^1`$ is the number of sites per unit surface area and $`m_i`$ the mass of the species (e.g., Tielens & Allamandola 1987). The binding energies $`E_{b,i}`$ are taken from Hasegawa & Herbst (1993), Aikawa et al. (1997), Willacy et al. (1998), and Fraser et al. (2001) (see also references therein). Now, the temperature in the clump is a decreasing function of distance from the central star and assumed to be unchanged with time, while the evaporation rate is very sensitive to the temperature as can be seen from Eq. (4). For example, in the case of water, whose binding energy is $`E_{b,\mathrm{H}_2\mathrm{O}}/k=5773`$ K, a temperature rise from $`T=95`$K to 105K shortens the evaporation time by about 1/100. Thus, in this paper the evaporation process is simply modeled:- each species, $`i`$, evaporates instantaneously inside a radius $`r_{\mathrm{inj},i}`$ where the evaporation time $`\tau _{\mathrm{evap},i}`$ is shorter than the accretion time,

$$\tau _{\mathrm{acc},i}=(S\pi a^2d_gnv_i)^1,$$
(5)
while it is retained on the dust otherwise. That is, the gas-phase abundances of mantle molecules are initially set as those listed in Table 1 inside the injection radii $`r_{\mathrm{inj},i}`$, and zero outside these radii. Beyond a radius of 1 pc, the density becomes so low that the accretion time becomes comparable to or greater than the chemical time-scales so that mantle formation is not efficient. We do not model this extended envelope here. In Eq. (5) $`n`$, $`a`$, $`d_g`$, and $`v_i=(2kT/m_i)^{1/2}`$ are the number density of hydrogen nuclei, the grain radius, the ratio of the number density of grains to $`n`$, and the thermal velocity of the species, $`i`$, respectively. The sticking probability $`S`$ is set as $`S=0.3`$. For the grain surface area per hydrogen nucleon, $`\pi d_ga^2`$, we used the value of $`<d_ga^2>=2.2\times 10^{22}`$cm2 (see Rawlings et al. 1992 for the parameter values). The corresponding injection radii $`r_{\mathrm{inj},i}`$ in the model of G34.3+0.15 (see next section) are also listed in Table 1. The binding energies of CO and N2 are low enough and the kinetic temperature high enough that they exist in the gas-phase throughout the region modelled here.

The subsequent gas-phase chemistry is calculated time-dependently by means of the chemical network based on that of Millar et al. (1997), which consists of 209 species connected by 2191 reactions. The standard interstellar cosmic-ray ionization rate, $`\zeta =1.3\times 10^{17}`$s-1, is adopted. The adsorption process of neutral atoms and molecules $`i`$ onto grain surface is included only if $`\tau _{\mathrm{acc},i}<\tau _{\mathrm{evap},i}`$ in this calculation. The accretion rate is given by

$$k_{\mathrm{acc},i}=S\pi a^2d_gnv_in_i,$$
(6)
where $`n_i`$ is the number density of species $`i`$. This process affects the gas-phase abundances slightly in the outer region, but is not very significant because the accretion rate of each species is comparable to the chemical reaction rates.

## 4 Model of hot core G34.3+0.15

In this section, we model the physical and chemical structure of the hot molecular core G34.3+0.15 using the radiative transfer and the pseudo-time-dependent gas-phase chemistry calculations described in the previous sections together with observations of the spectral energy distribution and the molecular line emissions toward G34.3+0.15. The hot core associated with this source consists of a small clump of hot gas situated at the nose of a cometary-shaped, ultra-compact Hii region. It has been well-observed, particularly through spectral line surveys (Macdonald et al. 1996; Kim et al. 2000) at 1, 2 and 3 mm, and has a rich molecular line spectrum, with large abundances of complex, highly-saturated molecules.

### 4.1 Density and temperature profiles

The density and temperature profiles of the hot core G34.3+0.15 are modeled so that the theoretical calculation of the spectral energy distribution (SED) is consistent with that observed. We take $`r_\mathrm{c}`$ defined in Sect. 2.1 as a parameter and obtain its best-fit value by iteratively calculating the temperature profile (see Sect. 2.2) and comparing the resulting SED with observation. The luminosity of the central star is set as $`L_{}=4.7\times 10^5L_{}`$, referring to Watt & Mundy (1999). For the temperature and radius of the star we adopt $`T_{}=4.5\times 10^4`$K and $`R_{}=11R_{}`$, assuming a massive zero-age main-sequence star (e.g., Thompson 1984; see also references therein). The distance to G34.3+0.15 of 3.7 kpc and the observational data for the SED are taken from Chini et al. (1987). As a result of the calculations, we find that the model with $`r_\mathrm{c}=0.05`$ pc represents the best fit to the observational data at wavelengths from far-infrared to millimeter. The resulting density and temperature profiles, and the calculated SED are shown in Figs.1a, 1b, and 1c, respectively, with solid lines. The crosses in Fig.1c are the observational data.

The data at wavelengths from near- to mid-infrared can not be reproduced by this model probably because of the inappropriate assumption of the spherically symmetric density profile. For example, the optical depth of non-spherical cavities due to bipolar outflows would be less than that of the dense envelopes so that if we were to observe the inner region of the clumps through optically thin cavities, namely, there would be more flux in the infrared (cf. Hatchell et al. 2000; van der Tak et al. 2000a). Actually, the model with the column density of $`N_\mathrm{H}=10^{24}`$ cm-2 (1/10 of the above envelope model) and the core radius of $`r_\mathrm{c}=0.05`$ pc, which is plotted with a dotted line in Fig. 1c, better represents the observational data in the mid-infrared. Here, the bolometric luminosity ratio contributed from the envelope region with $`N_\mathrm{H}=10^{25}`$ cm-2 to the cavity region with $`N_\mathrm{H}=10^{24}`$ cm-2 is assumed to be $`20:1`$ in order to fit the observational data. This will mean that the solid angle of the cavity will be much less than $`4\pi `$ and/or the radiation flux in the clump around the central star is not spherically symmetric in fact. Its density and temperature profiles are also shown with dotted lines in Figs. 1a and 1b. In the following subsections, we only use the model with the column density of $`N=10^{25}`$ cm-2 because most of the mass of the clump will be contained in this dense envelope region rather than the less dense, cavity region.

### 4.2 Results of chemical calculation

#### 4.2.1 Time evolution of molecular column densities

Making use of the density and temperature structures obtained in the previous subsection and the chemical model in Sect. 3, we simulate the time evolution of the hot core G34.3+0.15. The resulting time evolution of molecular column densities are shown in Fig.2. In order to compare these with the observations of beam-averaged column densities, each atomic or molecular column density is calculated as

$$N_i(t)=\frac{1}{\pi \mathrm{\Theta }_{\mathrm{beam}}^2}_0^{\mathrm{}}2\pi r\stackrel{~}{N_i}(r,t)\mathrm{exp}\left(\frac{r^2}{\mathrm{\Theta }_{\mathrm{beam}}^2}\right)𝑑r,$$
(7)
where $`\stackrel{~}{N_i}(r,t)=_{\mathrm{}}^{\mathrm{}}n(s)x_i(s,t)𝑑s`$ is the column density of species $`i`$ along a line of sight at a distance $`r`$ from the clump centre. Here, $`n(s)`$ is the number density of hydrogen nuclei, $`x_i(s,t)`$ the fractional abundance of species $`i`$ at time $`t`$, and $`s`$ the depth along the line of sight. The scale size $`\mathrm{\Theta }_{\mathrm{beam}}`$ is simply taken as 0.25 pc which corresponds to the angular size of 13” at the distance of 3.7kpc, based on the half power beam width of the James Clerk Maxwell Telescope (JCMT) at 345 GHz (cf. Thompson et al. 1999). While the data of most molecular lines we use in this paper are from single radio telescope observations such as the JCMT, the NH3 data is from the Very Large Array observation whose synthesized beam size is about 1”3 (Heaton et al. 1989). Thus, we calculate the NH3 column density with $`\mathrm{\Theta }_{\mathrm{beam}}=0.025`$ pc as well as with $`\mathrm{\Theta }_{\mathrm{beam}}=0.25`$ pc for comparison.

From Fig.2 we can see that parent molecules, such as H2O, NH3, H2S, and CH3OH, which are expected to be formed on grain surfaces, are relatively stable for timescales of at least 104 yr in this model. However, they are eventually destroyed by ionized species or atomic hydrogen to produce daughter species (see Sect. 4.2.3 for the details) as has been suggested in the previous studies (e.g., review by Millar 1997). Thus, the observed large abundances of both parent and daughter molecules will mean that the ages of the chemically rich hot cores are about 104 yr. In particular, we compare the results of our calculation with observations in Fig.3, which shows that the calculated molecular column densities from $`3\times 10^3`$ yr to $`3\times 10^4`$ yr (solid lines) are in good agreement with most of those observed toward the hot core G34.3+0.15 (diamonds plus dotted arrows if they are the lower or upper limits), consistent with the result of Millar et al. (1997). Lower limits in this Table have been derived from optically thick transitions. The calculated molecular column densities, with the exception of CH3CCH, HNC, H2CS, and CH3OCH3, satisfy the observational lower limits or fit the observations to within a factor 4 at most. Discrepancies may be resolved in a number of ways. Dimethyl ether, CH3OCH3, may in fact be a mantle molecule, while the gas-phase chemistry of H2CS is not at all well understood. Unlike H2CO, for which efficient gas-phase synthetic pathways have been studied in the laboratory, the analogous reactions which might form thioformaldehyde do not occur. The chemistry of HNC at low temperatures is well understood; the discrepancy between model and observation may imply that additional high-temperature loss mechanisms which are not included in the model are at work. Finally, the methyl acetylene, CH3CCH, column density is about a factor of 100 less than that observed by Macdonald et al. (1996) from their 330–360 GHz spectral line scan. However, it is much closer to the value of 8.2 $`\times `$ 1014 cm-2 determined from the 2 and 3 mm line surveys of Kim et al. (2000), similar to the value derived by Hatchell et al.(1998a) 20′′ offset from the hot core. Finally, we note that the formation of C3H$`{}_{}{}^{+}{}_{5}{}^{}`$, the precursor ion to CH3CCH, involves ions and neutrals related to CH4, C2H2 and C2H4. An increase in the abundances of each of these parents by a factor of a few would also help resolve the discrepancy noted above. The observational data are taken from Millar et al. (1997) (see also Macdonald et al. 1996 and Hatchell et al. 1998a) and Ikeda et al. (2001).

#### 4.2.2 Radial profiles of molecular abundances

In order to see the effects of the temperature dependent injection of molecules from dust grains on hot core chemistry, we present the radial profiles of molecular abundances at $`10^4`$ yr in Fig.4. We can see from the figures that the abundances of many parent molecules injected from grain mantles, such as C2H2, CO2, H2O, H2S, and CH3OH, change dramatically at $`r_{\mathrm{inj},i}`$, inside which the molecule $`i`$ can be easily evaporated. Inside $`r_{\mathrm{inj},i}`$, most of the parents are difficult to destroy. Nevertheless, given their large initial abundances, even a slow destruction rate is large enough to produce appreciable abundances of daughter products. Thus, the radial abundance profiles of the parent molecules influence those of daughter and other related species (see Sect. 4.2.3 for details). Abundance changes of some molecules along the radial direction in the clumps around young massive stars have been suggested by observational studies. For instance, the observations of the CH3CN rotational transition lines indicate that they are emitted from the inner compact hot region of the clumps, while the observed CH3OH lines suggest that the gas in some clumps consists of at least two components, one of which is hot and abundant in CH3OH, and the other is cooler and less abundant (Hatchell et al. 1998a; Millar & Hatchell 1998; van der Tak et al. 2000b). In addition, the observational abundances of the vibration-rotation lines of C2H2, HCN, H2O, and CO2 are found to increase with the excitation temperature (Lahuis & van Dishoeck 2000; Boonmann et al. 2000). Detailed radiative transfer calculations will be needed to compare the results of our theoretical model with these observations. We note that abundance jumps similar to those predicted have been observed in the low-mass ‘hot core’ around IRAS 16293-2422 (Schöier et al. 2002).

Our results show that a different chemistry operates in the surrounding cooler envelopes from that in the central hot regions where icy mantles can be evaporated, as has been indicated (e.g., Millar et al. 1997; Doty et al. 2002). In order to see the quantitative differences comparable with observations, we calculate the molecular column densities off-centred from the clump centre by a distance $`r_{\mathrm{off}}`$, making use of Eq. (7) with the distance from the clump centre, $`r^{}=rr_{\mathrm{off}}`$ replacing $`r`$. The results at $`t=10^4`$ yr are listed in Table 2. Parameters of $`\mathrm{\Theta }_{\mathrm{beam}}=0.25`$ pc and $`r_{\mathrm{off}}=0.4`$ pc are adopted for comparing the results with observations toward the halo associated with G34.3+0.15 by Thompson et al. (1999), which are also included in Table 2. For comparison, the molecular column densities calculated around the clump centre by Eq. (7) with $`t=10^4`$ yr are also listed together with the observations which are the same as those in Fig.3. The results of the model calculations are consistent with most of the observations. The ratio of the calculated molecular column densities around the clump centre to those off-centred range roughly from a factor 2, when the radial abundances of the species are almost constant, to one order of magnitude, in the case of parent or daughter species which are hardly generated by the gas-phase reactions without grain mantle evaporation.

#### 4.2.3 Detailed chemistry

In this subsection, we discuss some of the hot core chemistry, mainly focusing on the effects of the injection of parent molecules from grain mantles on the radial abundance profiles of daughter and other related species.

Fig.2(a) and Fig.2(b) show the time evolution of column densities of the carbon- and oxygen-bearing species, while Fig.4(a) and Fig.4(b) represent the radial profiles of their abundances at $`t=10^4`$ yr. The results related to oxygen chemistry can be interpreted based on the detailed analysis by Charnley (1997). Water is the main parent molecule among oxygen-bearing species and mainly destroyed by HCO+ to produce H3O+ in this model. This process keeps the HCO+ abundance low where H2O is injected from the grain mantles and is consistent with the result of Doty et al (2002). Hydroxyl molecules are formed by the dissociative recombination of H3O+, so that the OH abundance changes discontinuously at the H2O injection radius, $`r_{\mathrm{inj},\mathrm{H}_2\mathrm{O}}`$. One of the destruction processes of carbon monoxide is the reaction with He+, which makes atomic oxygen as well as atomic carbon. On the other hand, O and OH react with molecular hydrogen to reproduce water as O $`\stackrel{\mathrm{H}_2}{}`$ OH $`\stackrel{\mathrm{H}_2}{}`$ H2O at high temperature. Thus, the abundances of O and OH are suppressed in the inner hot region. These behaviours are also consistent with the results of Doty et al (2002). The reaction between OH and N generates NO, which makes the radial abundance profile of NO similar to that of OH. Molecular oxygen is not only injected from the grain mantles but also synthesized from O and OH. HCO is principally produced by the destruction of H2CO by OH, so that the HCO abundance jump appears at $`r_{\mathrm{inj},\mathrm{H}_2\mathrm{CO}}`$. C3H4 is chiefly formed initiated from C2H2 as C2H2 $`\stackrel{\mathrm{HCO}^+}{}`$ C2H$`{}_{}{}^{+}{}_{3}{}^{}`$ $`\stackrel{\mathrm{CH}_4}{}`$ C3H$`{}_{}{}^{+}{}_{5}{}^{}`$ $`\stackrel{e}{}`$ C3H4 in this model. Therefore, the C3H4 abundance increases over a radial distance between $`r_{\mathrm{inj},\mathrm{H}_2\mathrm{O}}`$, inside which the HCO+ abundance is suppressed, and $`r_{\mathrm{inj},\mathrm{C}_2\mathrm{H}_2}`$.

Fig.2(c) and Fig.4(c) are the time evolution of column densities and the radial abundance profiles of nitrogen-bearing species, whose detailed hot core chemistry have analysed by Millar et al. (1997) and Rodgers & Charnley (2001). The nitrogen molecule is the most dominant nitrogen-bearing species in the initial gas in this model, and is converted into N and N+, through reaction with He+. Ammonia is one of the possible nitrogen-bearing mantle species and is destroyed by ionized molecules to produce NH2. Meanwhile, it is formed from gas-phase reactions initiated from N+ as N+ $`\stackrel{\mathrm{H}_2}{}`$ NH+ $`\stackrel{\mathrm{H}_2}{}`$ NH$`{}_{}{}^{+}{}_{2}{}^{}`$ $`\stackrel{\mathrm{H}_2}{}`$ NH$`{}_{}{}^{+}{}_{3}{}^{}`$ $`\stackrel{\mathrm{H}_2}{}`$ NH$`{}_{}{}^{+}{}_{4}{}^{}`$ $`\stackrel{e}{}`$ NH3. The radial abundance profile of HCN has two apparent discontinuous jumps in this model. One is at the radius where $`\tau _{\mathrm{acc},\mathrm{HCN}}=\tau _{\mathrm{evap},\mathrm{HCN}}`$ and the other is at $`r_{\mathrm{inj},\mathrm{H}_2\mathrm{CO}}`$ as H2CO is one of the main precursors of CH2, which reacts with N to produce HCN. The major production reaction of HNC is that between NH2 and C. Thus, the HNC abundance jumps at $`r_{\mathrm{inj},\mathrm{NH}_3}`$. CH3CN is mainly produced by the radiative association between HCN and CH$`{}_{}{}^{+}{}_{3}{}^{}`$, so that the radial abundance profile of CH3CN resembles that of HCN. The CH3CN abundance is suppressed inside $`r_{\mathrm{inj},\mathrm{H}_2\mathrm{O}}`$ because H$`{}_{}{}^{+}{}_{3}{}^{}`$, which helps produce CH$`{}_{}{}^{+}{}_{3}{}^{}`$ from its reaction with CH3OH, decreases in abundance as it preferentially transfers its proton to the abundant water molecule. HC3N is generated from N as N $`\stackrel{\mathrm{C}_3\mathrm{H}}{}`$ C3N $`\stackrel{\mathrm{XH}^+}{}`$ HC3N+ $`\stackrel{\mathrm{H}_2}{}`$ H2C3N+ $`\stackrel{\mathrm{e}}{}`$ HC3N. C3H is a product of C2H2, which causes the HC3N abundance jump at $`r_{\mathrm{inj},\mathrm{C}_2\mathrm{H}_2}`$. CN is also chiefly formed from H2C3N+ inside $`r_{\mathrm{inj},\mathrm{C}_2\mathrm{H}_2}`$, but the jump does not appear in the CN abundance because CN can be synthesized by neutral-neutral reactions such as C+NO outside $`r_{\mathrm{inj},\mathrm{C}_2\mathrm{H}_2}`$. The CN abundance is suppressed in the inner hot region due to the loss by a reaction with H2 which creates HCN.

Fig.2(d) and Fig.4(d) represent the time evolution of column densities and the radial abundance profiles of sulfur-bearing species, respectively. Their basic behaviour can be explained following the analyses by Millar et al. (1997) and Charnley (1997). Hydrogen sulfide is a possible sulfur-bearing grain mantle species, and can be easily destroyed by atomic hydrogen especially in the inner hot region because of the small activation energy barrier of this reaction. This drives the sulfur chemistry in hot cores as H2S $`\stackrel{\mathrm{H}}{}`$ SH $`\stackrel{\mathrm{H}}{}`$ S $`\stackrel{\mathrm{O}_2}{}`$ SO $`\stackrel{\mathrm{OH}}{}`$ SO2. Thus, the SO abundance is lowered inside $`r_{\mathrm{inj},\mathrm{H}_2\mathrm{O}}`$, where the OH abundance is high, while the SO2 abundance drops outside $`r_{\mathrm{inj},\mathrm{H}_2\mathrm{O}}`$. Also, both the SO and SO2 abundances change dramatically at $`r_{\mathrm{inj},\mathrm{O}_2}`$. NS is produced by the reaction between SH and N, so that the NS abundance changes discontinuously at $`r_{\mathrm{inj},\mathrm{H}_2\mathrm{S}}`$. CS is mainly generated as SO $`\stackrel{\mathrm{C}}{}`$ CS in the inner region, while HCS+ $`\stackrel{\mathrm{e}}{}`$ CS in the outer region. On the other hand, HCS+ is synthesized by the proton transfer process to CS in the inner region, whereas by the reaction between CS+ and H2 in the outer region. Therefore, the radial abundance profiles of HCS+ and CS resemble each other. OCS is injected from grain mantles in this model in order to represent the observed high abundance of OCS molecular lines, and based on the observations of ice absorption features by Palumbo et al. (1997) (see also discussions in Hatchell et al. 1998b).

Fig.2(e) and Fig.4(e) show the time evolution of column densities and the radial abundance profiles of large oxygen-bearing species, which can be interpreted according to the analyses by Millar et al. (1991) and Charnley et al. (1992). Methanol is one of the possible large oxygen-bearing grain mantle molecules and is protonated by H$`{}_{}{}^{+}{}_{3}{}^{}`$. The product of the protonation, CH3OH$`{}_{}{}^{+}{}_{2}{}^{}`$, leads to the formation of daughter molecules as CH3OH$`{}_{}{}^{+}{}_{2}{}^{}`$ $`\stackrel{\mathrm{H}_2\mathrm{CO}}{}`$ HCOOCH$`{}_{}{}^{+}{}_{4}{}^{}`$ $`\stackrel{\mathrm{e}}{}`$ HCOOCH3 and CH3OH$`{}_{}{}^{+}{}_{2}{}^{}`$ $`\stackrel{\mathrm{CH}_3\mathrm{OH}}{}`$ C2H6OH+ $`\stackrel{\mathrm{e}}{}`$ (CH3)2O. Thus, abundance jumps of HCOOCH3 and (CH3)2O appear at $`r_{\mathrm{inj},\mathrm{CH}_3\mathrm{OH}}`$. CH3CHO is mainly made by the reaction between atomic oxygen and C2H5 which is generated from C2H6, so that the CH3CHO abundance changes dramatically at $`r_{\mathrm{inj},\mathrm{C}_2\mathrm{H}_6}`$. CH2CO is chiefly formed from the reactions of C2H3 $`\stackrel{\mathrm{O}}{}`$ CH2CO and CH$`{}_{}{}^{+}{}_{3}{}^{}`$ $`\stackrel{\mathrm{CO}}{}`$ CH3CO+ $`\stackrel{\mathrm{e}}{}`$ CH2CO inside $`r_{\mathrm{inj},\mathrm{C}_2\mathrm{H}_4}`$, while from the latter reaction otherwise. In this model C2H5OH is assumed to be created on dust grains in order to reproduce the observations of C2H5OH molecular lines (cf. Millar et al. 1997, Rodgers & Charnley 2001, Ikeda et al. 2001).

## 5 Discussion

### 5.1 Trapping of mantle molecules in water ice

It is known that water is one of the most abundant ice component in dense clouds. Some observations of ice absorption features in infrared spectra towards young stellar objects have suggested that a certain amount of CO is frozen in H2O-rich matrices (e.g., review by Allamandola et al. 1999), and laboratory experiments have shown that solid CO is able to diffuse into the porous structure of H2O and be trapped in it (Collings et al. 2003). These studies suggest that mantle molecules other than CO also could be trapped in water ice. Thus, in this subsection we examine the effects of this trapping on the chemical structure of the hot core G34.3+0.15. In order to see the extreme case, we put the injection radius of each mantle molecule as $`r_{\mathrm{inj},i}r_{\mathrm{inj},\mathrm{H}_2\mathrm{O}}`$, and assume that it is completely trapped in H2O matrices outside $`r_{\mathrm{inj},\mathrm{H}_2\mathrm{O}}`$, except CO, for which we assume that $`x(\mathrm{CO})=0.3x(\mathrm{H}_2\mathrm{O})`$ is trapped outside $`r_{\mathrm{inj},\mathrm{H}_2\mathrm{O}}`$ based on observations (e.g., Chiar et al. 1998; cf. Doty et al. 2002). As a result of the time-dependent chemical calculation described in Sect. 3, we get the molecular column densities of the hot core consistent with most of the observations at around $`10^4`$ yr (Table 4). The initial condition listed in Table 3 is used in this calculation, which is slightly different from that in Table 1, in order to fit the results to the observations. Although there are few differences in the integrated molecular column densities, we find that very different structure appear in the radial profiles of molecular abundances at around $`10^4`$ yr between the results in Sect. 4 (Fig.4: Model A) and those obtained from this calculation (Fig.5: Model B). Especially the differences are remarkable in the following molecules: (1) the binding energies of the molecules or their parent species are low, that is, the injection radii listed in Table 1 are large, such as H2CO, HCO, H2S, and SO, or (2) daughter molecules whose production are suppressed where H2O is evaporated from grain mantle, such as CH3C2H and CH3CN. For comparison, the resulting radial abundance profiles of SO and CH3CN from Model A (dotted lines) and B (solid lines) are plotted in the linear radial scales in Fig.6. Detailed observations of radial profiles of molecular abundances may make it possible to examine what fraction of molecules is trapped in water ice in cooler region of hot cores.

### 5.2 Dependence of chemical structure on density profile – General model

The comparison between the results of chemical calculation and the observations of various kinds of molecular lines will be useful to constrain the physical condition of massive star forming regions (e.g., Viti & Williams 1999). The density profile of hot molecular cores is one of the possible tracers of the evolution process just before and after massive star formation. In this section we will investigate the dependences of the chemical structure of hot cores on the density profile to discuss the physical condition of massive-star-forming clumps.

The chemical structures of three different types of molecular clupms are investigated here. Their density profiles are set by Eq. (1) with core radii of $`r_\mathrm{c}=`$ 0.01, 0.05, and 0.1 pc, and plotted in Fig.7(a) with dotted, solid, and dashed lines, respectively. Making use of these density profiles, we calculate the temperature profiles by solving the radiative transfer equation as described in Sect. 2.2. The resulting temperature distributions are shown in Fig.7(b). We can see from the figure that the temperature profiles directly reflect the difference in the density profiles. This is because the temperature profiles of hot cores are determined by the balance between the absorption of radiation from the central star and the reemission of the radiation by dust grains in the core as described by Eq. (2). The inner regions of the cores are optically thick enough to the radiation from the central stars so that the photons from the stars diffuse via interactions with the surrounding dusty material. If we simply assume that the dust opacity depends on the frequency of radiation as $`\kappa _\nu \nu ^n`$ and the density profile of a clump is set as $`\rho r^p`$, for example, the temperature profile of the clump illuminated by a central star is derived as $`Tr^{(p+1)/(4n)}`$ in the optically thick limit from the equation of local radiative equilibrium and the radiative transfer equation (e.g., Adams & Shu 1985). Meanwhile, since the evaporation process of grain mantle species is very sensitive to the temperature as shown in Sect. 3, it is to be expected that the density profiles affect the radial abundance profiles of molecules in hot cores via the temperature profiles. In order to see the dependence on the density models clearly, the temperature profiles in the region of grain mantle evaporation are plotted with linear scales in Fig.7(c).

The luminosity of the central star is another factor that influences the temperature profile of the clump. As we can see from Fig.7(b) and (c), the temperature distribution of the clump whose core radius is $`r_\mathrm{c}=0.05`$ pc and whose central star luminosity is $`L_{}=10^3L_{}`$ (dot-dashed lines) is similar to that of the model with $`r_\mathrm{c}=0.01`$ pc and $`L_{}=10^5L_{}`$. The difference in luminosity will be, however, distinguishable by observing the spectral energy distribution. The hydrogen column density will not have very much effect on the temperature structure in the optically thick region because basically $`T(r)`$ depends not on the absolute density but on the steepness of the density profile as we have shown in the above example. The optically thick region spreads (shrinks) if the hydrogen column density is high (low) so that it will affect the temperature profile at the outer radii.

We calculate the gas-phase chemical reactions in the clumps with the above density and temperature profiles and $`L_{}=10^5L_{}`$ by means of the chemical model in Sect. 3. Fig.8 gives the resulting radial abundance profiles of some molecules at $`10^4`$ yr in the linear radial scales (while those in Fig.4 are scaled in the logarithmic manner). They are parent or daughter species and have dramatic changes at certain radii related to grain mantle evaporation. H2O, CO2, H2CO, H2S, C2H2, and CH3OH are possible parent molecules, and have abundance jumps at their injection radii. The daughter species HC3N, CH3CN, and HCOOCH3 are produced mainly from C2H2, HCN, and CH3OH, respectively, in this model, and their abundance jumps appear at the same radii as those of the parent species. Comparison between reasonably high resolution observations of abundance jumps and precise radiative transfer and physical-chemical models of hot cores may make it possible to test the formation processes of these molecules.

From Fig.8 we can see that the abundance jump of each molecule occurs at a smaller radius in the clump with a smaller core radius $`r_\mathrm{c}`$, that is, with a steeper density and therefore a steeper temperature profile. The molecules in Fig.8 can keep high fractional abundances only in the inner region of the clump ($`<0.1`$pc) if $`r_\mathrm{c}=0.01`$pc. This result suggests that the chemically rich hot cores, whose observational beam-averaged molecular column densities are large, have flat density profiles rather than steep ones. On the basis of the low and high mass star formation theories, the central density profiles of clumps are thought to be steep just after the star formation, at least if accretion of the surrounding material to the central protostellar cores is important (e.g., Shu 1977; Yorke & Bodenheimer 1999; Masunaga & Inutsuka 2000; Yorke & Sonnhalter 2002). Subsequently, dissipation by outflows will let the clumps gravitationally unbound and their central density profiles will be flattened (e.g., Shu et al. 1987; Nakano et al. 1995). In the case of massive star formation, powerful outflows with high mass loss rates as well as radiation pressure and strong winds from the stars may accelerate the dissipation process. If massive stars are formed as a result of coalescence of lower mass stars, the envelopes will have different density profiles from those of isolated star-forming clumps (e.g., Bonnell et al. 1998). Statistical observations of various kinds of molecular lines towards hot cores could constrain their density profiles and the massive star formation theories.

The effects of inflow motion of hot cores have not been considered in this paper. If the inflow has a velocity of a few km s-1 as is estimated from the observed atomic and molecular lines (e.g., Keto 2002), it will not affect the structure on the scale of about 0.25 pc, which corresponds to the observational beam size for most of the molecular lines we have used in this paper, within the timescale of $`10^4`$ yr. Comparison between model calculations and the forthcoming high spatial resolution and high sensitivity observations could reveal the dependence of the chemical structure of hot cores on the inflow motion. In addition, the temporal evolution of young massive stars in hot cores is also a factor that may possibly influence the fractional abundances of some species (cf. Viti & Williams 1999).

## 6 Summary

We have constructed the density and temperature profiles of hot cores that have embedded luminous objects, using radiative transfer calculations together with observations of the spectral energy distribution. We have used these to investigate hot core chemistry taking into account the temperature dependent grain mantle evaporation process.

As a result of applying the model to the hot core G34.3+0.15 and calculating the time-dependent gas-phase reactions in hot gas heated by a central massive star, it is confirmed that the age of around $`10^4`$ yr is preferred to reproduce the large column densities of both the parent and the daughter molecules. In fact, the calculated column densities are consistent with most of those observed toward G34.3+0.15 around 104 years after central luminous star formation with a particular initial composition for the grain mantle molecules.

Also, we find that dramatic changes appear in the radial profiles of fractional abundances of some parent, daughter and related species at the parents’ injection radii, inside which the evaporation time of the molecules from dust grains are shorter than their accretion time. The abundance jumps appear because (1) most of the parents are difficult to destroy within the timescale of 104 yr, and (2) the initial abundances of parents is large enough, so that despite their slow destruction rates inside the injection radii, they produce appreciable abundances of daughter products, while (3) these parent and daughter molecules are hardly produced by the gas-phase reactions outside the injection radii.

Finally, we discussed the influence of some chemical and physical properties on the chemical structure of hot cores. First, we investigated the effects of the trapping of mantle molecules in water ice. As a result, we find large changes in the location of jumps in the radial abundance profiles of some species and that it changes completely the radial abundance profiles of daughter molecules whose production are suppressed in regions where H2O is evaporated from grain mantle. Also, we investigated the chemical structure of hot cores with different density profiles. Our results suggest that the chemically rich hot cores have flat density profiles rather than sharp ones. This is because a steeper radial density profile leads to a steeper temperature distribution as long as the clump is optically thick to radiation from the central star. This means that the abundance jumps appears at smaller radii and that therefore beam-averaged molecular column density cannot be large. Observational possibilities of constraining the density profiles of hot cores and the massive star formation theories are also suggested.

###### Acknowledgements.

We would like to thank the referee, Dr. F.L. Schoier, for his comments which improved our paper. Astrophysics at UMIST is supported by a grant from PPARC.


text:
# Oscillatory Motion of Rising Bubbles in Wormlike Micellar Fluids with Different Microstructures

(October 15, 2003; submitted January 10, 2003)

## Abstract

Previous observations of the nontransient oscillations of rising bubbles and falling spheres in wormlike micellar fluids were limited to a single surfactant system. We present an extensive survey of rising bubbles in another system, an aqueous solution of cetylpyridinium chloride and sodium salicylate, with and without NaCl, across a range of concentrations and temperatures. Two different types of oscillation are seen in different concentration ranges, each with its own temperature dependence. Rheological data allows for the identification of these different hydrodynamic states with different fluid microstructures.

Fluids are often broadly categorized as either Newtonian or non-Newtonian, according to whether the Navier-Stokes equation does or does not accurately describe the fluid’s motion. Newtonian fluids consist of molecules small enough to be approximated by point masses (such as air and water), while non-Newtonian consist of larger structures (such as polymers). This increased size affords greater degrees of freedom which leads to macroscopic viscoelasticity bird . Fluids consisting of self-assembling cylindrical (or wormlike) micelles resemble polymeric fluids on the microscopic scale, with the added feature that their length distribution is determined by aggregation kinetics; micelles continually break and reform israelachvilibook ; gelbart ; cates90a . The molecular level physics is clearly more complicated for such fluids, and one can reasonably expect this will introduce a new set of flow properties to the class of non-Newtonian fluids.

Several novel results have been reported in both viscometric flows rehage88 ; shikata88 ; grand ; lerouge00 ; liuandpine  and hydrodynamic (non-viscometric) flows smolka01 ; ab00 ; anand01  of various wormlike micellar fluids. Among these, a falling pendant drop at the bottom of a thin filament has been observed to slow to a complete stop before the filament suddenly ruptures smolka01 ; this rupture has been confirmed in other configurations john . Also, rising air bubbles and falling solid spheres have been observed to oscillate without reaching a terminal velocity ab00 ; anand01 . While in polymer solutions a rising bubble displays a sharp non-axisymmetric cusp, which remains unchanged during a steady rise bird ; liu95 , in wormlike micellar fluids the cusp periodically extends to a sharp point, then retracts to a blunt edge (see Fig.1). It is likely that these new dynamics are a manifestation of the reversible scission reactions of the micelles, though the precise mechanism may be quite complicated. Another type of fluid is known to include scission-like reactions at its microscale - associating polymer solutions - and falling sphere oscillations are reported in that system as well mollinger ; weidman .

In this Letter we report a detailed survey of the oscillatory motion of rising bubbles in a wormlike micellar fluid. To study the role of the aggregation kinetics of the micelles, both concentration and fluid temperature were varied over a wide range, unlike previous studies ab00 ; anand01 . As these parameters change, different micellar architectures are possible - from short linear or branched micelles to crosslinked networks appell92 ; schmitt94 ; in99 . Four different dynamics were seen: Newtonian behavior at high temperatures, standard polymeric behavior, and two distinct oscillating responses occurring in different concentration ranges (Fig.2). We performed steady rheology experiments to identify the fluid microstate, and found that transitions in the equilibrium structure match transitions in bubble dynamics. Critical temperature bounds were also found, which can be interpreted as a minimum length of micelle required for oscillations to occur, so that for certain concentrations the fluid may be tuned with temperature to make bubbles either oscillate, rise with a stable cusp, or rise as bubbles in a Newtonian fluid.

Non-transient oscillating bubbles were first observed in aqueous solutions of cetyltrimethylammonium bromide (CTAB) and sodium salicylate (NaSal) ab00 ; here we use another familiar system, cetylpyridinium chloride (CPCl) and NaSal, with the fixed ratio \[NaSal\]/\[CPCl\] = 1 (except for the NaCl experiments described below). Our choice for this ratio gives a solution of flexible cylindrical micelles analogous to polymer chains rehage88 . All experimental results are for a single rising bubble, with volumes ranging from 14 mm3 to 110 mm3. The bubble is injected into the fluid at the bottom of a temperature controlled plastic cylinder (31 cm height, 5 cm diameter) using a syringe with a long stainless steel tube (inner diameter of either 1.0 mm or 1.5 mm). The ratio of the horizontal diameter of the bubble $`d`$ to the cylinder diameter $`D`$ is $`d/D0.14`$, and we have also checked that bubbles oscillate for $`d/D0.02`$. Typical sizes of the Reynolds number (inertia) are $`Re10^25`$, while the Deborah number (elasticity) ranged from $`De1500`$, similar to values seen for oscillating bubbles and spheres in CTAB/NaSal ab00 ; anand01 . Rheological data were taken with a Rheometrics RFS-III controlled rate of strain rheometer with circulating fluid temperature bath, using a stainless steel Couette geometry. Video images were made with a Kodak high speed CCD camera.

We have studied concentrations from 4 to 40 mM (weight fractions $`0.13\%\phi 1.3\%`$). From 5 to 15 mM, bubbles have a cusp which oscillates in length and changes shape (Fig.1). While rising, the cusp lengthens (frames 1-4 in Fig.1) during which the velocity as measured at the top of the bubble is nearly constant. At the apex of the extension, the tail abruptly retracts and the bubble jumps upward (frames 4-5). After this it slows to a nearly constant velocity until the cycle repeats. Typically, bubbles begin to oscillate within 10 seconds of their formation, with a similar time between jumps. These clearly visible oscillations, which we call type I, are apparently not a transient effect; we have observed their persistence for rise distances over one meter, during which there were more than 30 oscillations in $``$ 35 s (for 8 mM CPCl/NaSal), and velocities more than doubled during a jump. This is the typical type I behavior in the temperature ranges shown in Fig.2, consistent with the oscillations observed in CTAB/NaSal ab00 .

Above a critical temperature in the type I range (Fig.2), bubbles have a sharp cusp which does not change in length (“steady cusp” in Fig.2) and velocities smoothly reach a steady state; the oscillatory instability vanishes. By $`60^{}`$C, the solutions appear Newtonian (“no cusp” in Fig.2): large enough bubbles are ellipsoidal and undergo the well-known side-to-side oscillations hart57 . In this temperature range, the micelles must be predominantly spherical or short rigid rods ben-shaul .

Upon increasing the concentration from 15 mM to 16 mM ($`\phi =0.5\%`$ to $`0.53\%`$), *all oscillations cease*. Bubbles rise steadily with a stable cusp for concentrations up to 20 mM; there are no temperatures at which oscillations occur (Fig.2). It is striking that a rising bubble would be so affected by so slight a concentration change. Note the transition temperature from polymeric to Newtonian behavior reaches a maximum at this concentration (Fig.2).

More surprising is the re-emergence of oscillations for concentrations from 25 mM to 40 mM ($`0.8\%\phi 1.3\%`$). Yet here the oscillations, which we call type II, are visibly different from type I oscillations; the shape change involves the entire bubble, whereas in type I it seems restricted to the tail. A type II oscillation begins with a constriction in width near the top, while the whole bubble lengthens. This constriction then travels downward, as if the bubble were squeezing through a hoop. Defining $`w_{max}`$ to be the maximum (relaxed) width (at its waist), $`w_{min}`$ to be its most constricted width, and $`\mathrm{\Delta }w=w_{max}w_{min}`$, we found $`\mathrm{\Delta }w/w_{max}=0.13`$ in type II fluids and $`0.04`$ in type I. Length extensions, however, were $`\mathrm{\Delta }l/l_{min}=0.25`$ in type II and $`0.27`$ in type I; for bubbles rising in CTAB/NaSal, $`\mathrm{\Delta }l/l_{min}=0.26`$ ab00 . This previous study revealed only one type of oscillation, which from the shape dynamics we identify as type I ab00 .

Although wormlike micelles are comprised of surfactants, it is unlikely that surface tension plays a role in the bubble oscillations, since falling rigid spheres also oscillate anand01 . Furthermore, our observation of two oscillation types suggests that the viscoelastic character of the bulk fluid is changing with concentration. We address this with rheology measurements, controlling the shear strain rate $`\dot{\gamma }`$ and recording the effective viscosity, $`\eta _{\text{eff}}=\sigma _{xy}/\dot{\gamma }`$, where $`\sigma _{xy}`$ is the steady shear stress. Transient tests were also performed to ensure that steady state was achieved. Shown in Fig.3a is $`\eta _{\text{eff}}`$ vs. $`\dot{\gamma }`$ at $`30^{}`$C, for fluids in the oscillating and non-oscillating concentration ranges. Most fluids are shear thinning - $`\eta _{\text{eff}}`$ decreases at high enough $`\dot{\gamma }`$ \- however fluids below 5 mM shear thicken ($`\eta _{\text{eff}}`$ increases above $`\eta _0`$), typically requiring $`100`$ seconds to reach steady state liuandpine . Bubbles tested in fluids below 5 mM showed no oscillations or steady cusps. While there is a noticeable difference in the rheology upon increasing from 20 to 30 mM, 30 and 40 mM (type II region) are strikingly similar; they also have a constant viscosity over a broader $`\dot{\gamma }`$ range than the other fluids. More interestingly, the zero shear viscosity $`\eta _0`$ (defined as the plateau value at low $`\dot{\gamma }`$) decreases with concentration.

A more complete study of the dependence of $`\eta _0`$ on $`\phi `$ indicates transitions in fluid microstructure (Fig.3b). The transition at low concentration to rapid $`\eta _0`$ growth marks the overlap concentration $`\phi ^{}3`$ mM, below which is the dilute regime. Above $`\phi ^{}`$ is the semi-dilute state, in which micelles exist as individual entangled worms berret93 . Here $`\eta _0\phi ^{5.8}`$, close to the value $`5.4`$ associated with stress relaxation by reptation appell92 ; candau93 . The semi-dilute regime ends near $`\phi 1215`$ mM, followed by an extraordinary decrease of $`\eta _0`$ with concentration, continuing to $`\phi 30`$ mM, after which $`\eta _0`$ varies weakly.

The dramatic change in the concentration dependence of $`\eta _0`$ indicates a transition in the equilibrium fluid structure, coinciding with the loss of type I oscillations. It may be that the entangled micelles in the semi-dilute range have begun to fuse at entanglement points, forming a crosslinked network appell92 ; candau93 ; drye ; pincus! . The junction nodes are free to slide along the micelles in such a formation, which would account for the decrease in viscosity appell92 . If the ratio of crosslinks to entanglement points grows as $`\eta _0`$ decreases, then the new state would be fully formed where $`\eta _0`$ stabilizes near $`30`$ mM ($`1\%`$), the start of type II oscillations. Note that a crosslinked network state has also been proposed for CTAB/NaSal at $`\phi 1\%`$ lequeux97 .

Transitions in micellar morphology would naturally lead to transitions in mechanisms for stress relaxation, which should be observable in the stress field around the rising bubble. These transitions are made evident with birefringent visualization, shown in Fig.4 for the three concentration regimes. Optical birefringence is a well-known technique for visualizing stress in non-Newtonian fluids opticalannrev , especially effective for wormlike micellar fluids hu93 . Fig.4a (10 mM) shows the localization of stress in the wake of a type I oscillating bubble. This birefringent tail mirrors the dynamics of the bubble’s cusp (Fig.1), and when it retracts, the birefringent tail suddenly disperses to the sides. The flow in the bubble wake is nearly in uniaxial extension harlen90 , and the birefringent tail suggests a strong relation to a thin filament whose rupture smolka01  may be related to the mechanism for oscillation. At 20 mM a shorter tail is seen (Fig.4b), whose length remains constant (like the bubble tail). Fig.4c corresponds to type II oscillations (35 mM), in which three equally spaced birefringent bands indicate a broader distribution of stress. This pattern is essentially unaltered during an oscillation. Note that $`De`$ decreases from type I to type II, with an intermediate value for non-oscillating bubbles (Fig.4 caption), contradicting the implicit critical $`De`$ condition for oscillations in anand01 .

The length distribution of the individual worms in the type I region (semi-dilute) should depend on temperature and concentration. In equilibrium, the average length of a wormlike micelle $`L_0`$ can be described with the mean field approximation cates90a ; israelachvilibook ; ben-shaul :

$$L_0(\varphi ,T)\sqrt{\varphi }e^{E/2kT}$$
(1)
where $`\varphi `$ is the total amphiphile concentration, $`k`$ is Boltzmann’s constant, $`T`$ is temperature, and $`E`$ is the scission energy required to break one wormlike micelle into two. Estimates based on light scattering measurements for $`E`$ are on the order of $`10kT`$ cates90a ; kwon ; in99  \- much lower than the covalent bonds of polymers, yet large enough for some micelles to reach appreciable lengths against thermal fluctuations cates90a . The flow near the bubble may increase the equilibrium length, and there is evidence that shear-induced structures (SIS) much larger than individual micelles form in such flow in both CTAB/NaSal liuandpine ; bolt97 ; skeller98  and CPCl/NaSal wheeler ; fischer . We assume simply that for the type I region, $`L_0`$ is the relevant quantity determining if the bubble will oscillate. Specifically, we attribute the oscillations to a breaking instability of the elongated micelles or SIS in the wake, occurring only if $`L_0`$ exceeds some critical value. Using Eqn.1, the solid line in Fig.2 corresponds to a scission energy of $`E=1.01\times 10^{19}`$ Joules, about $`24kT`$.

Type II oscillations have a much simpler dependence on temperature; they occur only for $`TT_c=36^{}`$C (Fig.2). If crosslinks have indeed formed for $`T<T_c`$, then the mean field length is no longer an appropriate quantity. The transition temperature $`T_c`$ suggests a critical energy condition ($`E_c4.3\times 10^{21}`$ Joules), which may correspond to the energy difference between junctions and endcaps. Thus for $`T>T_c`$, endcaps would dominate and the fluid is comprised of individual micelles which can entangle like polymers. This possible explanation is consistent with both the rheology and the rising bubble dynamics.

The two different oscillations we have observed appear linked to two different microstructures. We extended our study to include the well characterized ternary system CPCl-NaSal diluted in concentrated NaCl, with a ratio \[NaSal\]/\[CPCl\] = 0.5 berret93 ; berret94 ; lerouge00 . We tested several concentrations conc  of these fluids for jumping bubbles in our apparatus, at $`T=30^{}`$C, a temperature central to both type I and II oscillations (Fig.2). We observed no oscillations or any behavior different from steadily rising bubbles in conventional polymeric fluids liu95 . This classic ternary system is known to consist of entangled wormlike micelles, which when taken with the observed scaling $`\eta _0\phi ^{3.3}`$, indicates that micellar breaking occurs on a shorter timescale than reptation cates90a ; berret93 ; in contrast for our fluids $`\eta _0\phi ^{5.8}`$. Evidently, the oscillatory motion of rising bubbles is not a characteristic of fluids in this “fast-breaking” limit cates90a ; berret97 , and we conclude that type I oscillations are in the “slow-breaking” limit.

Our study of rising bubbles in various wormlike micellar fluids indicates that while scission reactions may be necessary for oscillations, there are other conditions. The discovery of a second type of oscillation provides another example of how microscale dynamics and architecture (entangled vs. crosslinked) interact to produce macroscopic instabilities. It now seems that wormlike micellar fluids are the most generic complex fluids as far as rising bubbles are concerned, since all known material-dependent dynamics can occur at different temperatures or concentrations. More surprisingly, the rise of an air bubble has been shown to be extraordinarily sensitive to fluid microstructure.

We thank A. Jayaraman, J. T. Jacobsen, P. Olmsted, T. Podgorski, and L. M. Walker for valuable discussions, R. Geist, D. M. Henderson, and M. C. Sostarecz for experimental assistance, and the referees for constructive comments. AB acknowledges support from the A. P. Sloan Foundation and the National Science Foundation (CAREER Award DMR-0094167).


text:
# Dark Matter Constraints from the Sagittarius Dwarf and Tail System

## 1. Observational Constraints

The relatively nearby Sagittarius (Sgr) dwarf galaxy offers the opportunity to explore in exquisite detail the interaction of a satellite with its parent galaxy. Moreover, the extensive Sgr tidal tail system gives sensitive leverage on the properties of the Galactic potential, much as polar ring galaxies have been exploited to determine the properties of extragalactic systems (e.g., Sparke 2002).

Since Sgr’s discovery by Ibata, Gilmore & Irwin (1994), this archetype of a dwarf galaxy merger has remained difficult to study because its core (centered at $`[l,b]=[6,14]^{}`$), lies obscured by foreground dust and stars of the Galactic disk and bulge. However, the Two Micron All Sky Survey (2MASS) has opened a new window on Sgr because of the reduced extinction by dust in the near infrared (NIR) and because Sgr contains a significant population of M giant stars, which are bright in the NIR. To improve the contrast of Sgr, these M giants can be separated from foreground M dwarfs using the combination of 2MASS $`JHK_s`$ photometry because of differential, surface gravity-sensitive opacity effects near 1.6$`\mu `$ (Bessell & Brett 1988). Majewski et al. (2003a; MSWO hereafter) used this method to make all-sky maps that reveal Sgr to be the primary source ($`>80\%`$) of M giants in the high halo ($`|Z_{GC}|>20`$ kpc), excluding those bound in the LMC and SMC (which actually show no M giant tidal tails themselves). The Sgr M giants lie along a great circle tipped only 13.5 from a truly polar ring.

A three dimensional analysis is made possible by converting $`K_s`$ magnitudes into photometric parallaxes using the red giant branch relation defined by the Sgr core (for which we adopt $`[mM]=16.90`$). The resulting three dimensional distribution reveals a flattened, planar alignment with a $`<2`$ kpc RMS spread, and clear trailing (predominantly south of the Galactic Plane) and leading (predominantly north of the Galactic Plane) tidal arms. Even without benefit of modeling, the debris arms show clearly that Sgr is orbiting within a non-precessing plane in an elliptical (approximately 3:1 or 4:1 apo:peri-Galacticon) orbit of $``$14 kpc peri-Galacticon. Together the tidal arms wrap at least 360 around the sky, and overlap at even greater length.

To constrain the dynamics of the Sgr-Milky Way interaction further, we (e.g., Majewski et al. 2003b, Paper II hereafter) have been collecting radial velocities of Sgr M giants with the Swope 1-m, KPNO 2.1-m, Yalo 1.5-m and Bok 2.3-m telescopes. As of July 2003, our sample of M giant velocities, distributed around the Sgr plane, is approaching 900 M giants. The $`5`$ km s-1 precision velocities provide critical input to Sgr models, and help compensate for remaining vagaries in the M giant spatial distribution due to contamination by disk/bulge M giants and random (and perhaps systematic) errors in the distance scale.

## 2. N-body Simulations

To understand in detail the Sgr interaction with the Milky Way (MW), we (Law et al. 2003b, Paper III) have been using N-body simulations based on those developed by Johnston et al. (1996, 1999). Early results of this work are described in Law et al. (2003a). For now these models are constructed with a focus on matching the structure and kinematics of the Sgr tails; future efforts will aim to understand the detailed nature of the Sgr core. Our immediate goal is to use the extensive tidal tail system revealed by 2MASS M giants as dynamical probes of the MW potential — its size, overall flattening, and lumpiness — and the global character of the Sgr dwarf — its orbit, mass, and dark matter content.

Sgr is represented by 105 self-gravitating particles (for both light and dark matter) distributed according to a Plummer model. This satellite orbits in a rigid MW potential represented by a Miyamoto-Nagai (1975) disk, Hernquist (1990) spheroid, and logarithmic halo — $`\mathrm{\Phi }_{halo}=v_{halo}^2\mathrm{ln}([R^2+(z/q)^2+d^2])`$ — constrained to fit the established Galactic HI/CO rotation curve interior to the Solar Circle. The solar distance to the Galactic Center, $`R_{}`$, the distance of Sgr, $`d_{Sgr}`$, and the circular speed of the Galaxy, $`v_{halo}`$, are varied but Sgr’s radial velocity is fixed at 171 km s-1 (Ibata et al. 1997, “I97”). The velocity vector of Sgr is constrained to lie within the orbital plane established by the M giants, and the satellite is evolved through the simulated Galactic potential for five orbits using the self-consistent field code of Hernquist & Ostriker (1992).

Both full N-body simulations as well as less CPU-intensive, test particle models111A single test particle orbiting in the Galactic potential with the same dynamical constraints. (when appropriate — e.g., to explore gross orbital properties) are used to explore the parameter space of Galactic potential strength and shape, Sgr orbit, and Sgr mass; the ranges of some parameters explored are given in Table 1. Models are evaluated by their ability to reproduce a set of observed properties of the M giant data, among the most important and discriminatory including the apo-Galacticon of the leading arm, the mean positional, velocity and density trends of the trailing arm, velocity and positional spreads of the debris, as well as the amount of precession in the arms (see §3). Some degeneracy of parameter combinations yielding reasonable fits to the data is found, but fixing $`q=1`$ and $`R_{}=8.5`$ kpc yields a “best-fitting” solution given by the “adopted model” in Table 2. The latter model (subject to change with further experimentation) is characterized by a 0.75 Gyr radial period Sgr orbit with a 14 kpc peri- and 52 kpc apo-Galacticon, and present space velocity of 326 km s-1. In this model the observed M giant tidal arms correspond to debris lost on at least 2.5 orbits over the last $``$1.8 Gyr or more. For any $`R_{}`$$``$8.5, the MW mass within 50 kpc is restricted to 3.7-5.1$`\times 10^{11}`$ M$`_{}`$, which is slightly smaller than the recent determination of $`5.5\times 10^{11}`$ M$`_{}`$ by Sakamoto, Chiba & Beers (2003) from an analysis of the velocities of hundreds of random halo objects. The model timescales and masses scale with the potential, while the apo- and peri-Galactica scale with the M giant distance scale, which was based on the assumed $`d_{Sgr}=24`$ kpc.

## 3. Constraints on the Halo Flattening

Because the Sgr orbital plane is not strictly polar, it should precess in a flattened potential. Yet the observed debris plane is remarkably well collimated, suggesting a nearly spherical MW halo at the distance of Sgr (Ibata et al. 2001, MSWO). To quantify this, separate planes are fit to $`120^{}`$ sections of the leading and trailing arm corresponding to debris $`180^{}`$ out of orbital phase (but up to $`300^{}`$ separated), and the angle between the planes measured. Errors reflect the quadrature sum of the errors in fitting the two planes. Table 2 gives the results of this comparison for the observed M giant debris and for closely matching simulations varying only by the degree of flattening in the halo potential.

As may be seen, the distribution of observed M giants is fully consistent with a $`q=1`$, spherical halo potential, although a 1$`\sigma `$ error permits a slight flattening of the halo ($`q=0.95`$). The constraint on halo flattening offered by the Sgr tidal arms will strengthen as the length over which the arms are traced increases, both by M giants (possible when we verify the radial velocity membership of potential M giants at even greater separations from the core — an observing program in progress) and by use of older stellar tracers that can track the arms beyond the length that can be traced with the $``$ 2-3 Gyr old Sgr M giants. Based on these results, we have adopted a $`q=1`$ halo potential in our simulations.

## 4. Constraints on the Bound Sgr Mass and $`M/L`$

The $`M/L_V`$ of Sgr has previously been found to be very large — of order 50 (I97) to 100 (Ibata & Lewis 1998). However, MSWO showed that the M giants of the Sgr core can be fit by a King model with much larger core and limiting radii than previously adopted. Inserting these radii and the 11.4 km s-1 central velocity dispersion (I97) into the standard King (1966) $`M/L`$ methodology that assumes virial equilibrium yields a $`4.9\times 10^8`$M$`_{}`$ bound Sgr mass that drops $`(M/L_V)_{tot}`$ to 25. However, as MSWO point out, there is little reason to believe even this $`M/L`$ because it is very unlikely that the Sgr tidal radius corresponds to the measured 12.6 kpc semi-major limiting radius. For example, adopting this as a tidal radius yields an absurdly discrepant $`1.6\times 10^{11}`$ M$`_{}`$ Sgr mass from the Roche tidal limit, $`m_{Sgr}=[2M_{MW}(R_{GC})][r/R_{GC}]^3`$. Even adopting the 4.4 kpc minor axis limiting radius still yields $`m_{Sgr}=6.9\times 10^9`$ M$`_{}`$ and $`M/L_v343`$. The tidal radius must be much smaller, especially if we are to explain how 2-3 Gyr old M giants, presumably formed in the central few kpc of Sgr, could so quickly have escaped across the tidal boundary into tidal arms of similar dynamical age. With no clear physical markers of a tidal radius in the spatial distribution of Sgr core stars, we must resort to alternative means to estimate its mass.

Although we do not attempt to model the Sgr core in detail, we can use the coherence of the Sgr stream to estimate the mass of the disintegrating parent core, since larger mass bodies produce commensurately “hotter”, wider debris trails. The simulated dwarf that best fits the spatial and velocity width of the streams has a $`3\times 10^8`$M$`_{}`$ mass within a semi-major tidal radius of about $`r_{\mathrm{tide}}=3.5`$ kpc, which gives $`(M/L_V)_{tot}=21`$ (adopting the Sgr $`L_V`$ within this radius). This may be regarded as an upper limit to the Sgr core mass (see §5).

When the Sgr arms can be mapped accurately over even greater lengths, it may be possible to map the degree to which dynamical friction has acted on the Sgr core, and thereby derive yet another, independent estimate of the Sgr mass. We have explored dynamical friction models as a means to explain one apparent discrepancy between the model and M giant velocities in the nearest sections of the Sgr leading arm, however no satisfactory results have been obtained that do not imply a huge recent mass loss and extreme former Sgr core mass.

## 5. Constraints on the Lumpiness of the Milky Way Halo

Current CDM models for the formation of structure in the universe predict that MW-like galaxies should contain substantial halo substructure at current epochs as a result of the accretion of thousands of subhalos over a Hubble time (e.g., Navarro et al. 1997). Because the MW currently has only eleven known luminous satellite galaxies, it is commonly held that the bulk of the subhalos must be made up of pure dark matter (see Klypin et al. 1999, Moore et al. 1999). If so, these dark matter lumps should make their presence known through the heating of dynamically cold, luminous stellar systems, like tidal tails (Moore et al. 1999, Font et al. 2001, Johnston, Spergel & Haydn 2002, Ibata et al. 2001, 2002).

The present velocity and spatial dispersion of the Sgr tails reflect both the initial dispersion of tidally released debris as well as subsequent heating of that debris imparted by encounters with large halo masses. Attributing all of the dispersion to one or the other of these phenomena provides upper limits to the effects of each. We (Paper II) have measured a $`10.4\pm 1.5`$ km s-1 velocity dispersion over a $`>`$100 expanse of the trailing Sgr arm, a dispersion nearly equivalent that of the Sgr core. If we attribute all of this dispersion to the Sgr central mass, then we derive an upper limit to that mass of 3$`\times 10^8`$ M$`_{}`$ and $`M/L_v<21`$ (as described above). However, if any of the velocity dispersion in the tidal arms is attributable to heating by subhalos, the implied bound mass of Sgr decreases. While the velocity dispersion of the trailing Sgr arm actually is observed to increase slightly with distance from the Sgr core (Paper II), it is not yet clear whether this is the signature of subhalo heating or simply the fact that the bound Sgr mass was larger in the past.

Johnston et al. (2002, “JSH”) give a prescription for a tidal debris “scattering index”, $`B`$, that measures position and velocity perturbations of tidal arm stars induced by lumps in the halo under the assumption of an initial zero velocity dispersion tidal debris population (thus the index provides an upper limit to the scattering for real debris). When applied to trailing arm M giants with 25-90 separation from the Sgr core we obtain $`B=0.031`$, a “colder” result than the 0.037 value JSH obtained for presumed Sgr carbon stars. The new, smaller $`B`$ is consistent with JSH simulations of heating in a smooth halo containing just one LMC-like (mass and orbit) lump; however, some realizations of lumpier halos in the JSH analysis are not inconsistent with the degree of scattering observed here. Unfortunately, this reflects a vagary of this type of halo probe: Dynamical heating in CDM halos tends to be dominated by the most massive lumps.

Nevertheless, still tighter constraints — in the direction of making lumpy halos even less likely — may derive from future observations and modeling of tidal streams. For example, accurately determining the zero-age dispersion of Sgr debris will make it possible to remove this contribution from any dispersion by heating. The study of initially colder streams, e.g., from globular clusters, would place even stricter constraints on the lumpiness of the halo.

We acknowledge funding from NASA, the National Science Foundation, and the David and Lucile Packard Foundation.

## References

Bessell, M. S. & Brett, J. M. 1988, PASP, 100, 1134

Font, A. S., Navarro, J. F., Stadel, J. & Quinn, T. 2001, ApJ, 563, L1

Hernquist, L. 1990, ApJ, 356, 359

Hernquist, L. & Ostriker, J. 1992, ApJ, 386, 375

Ibata, R. A, Gilmore, G. & Irwin, M. J. 1994, Nature, 370, 194

Ibata, R. A., Lewis, G. F., Irwin, M. J. & Quinn, T. 2002, MNRAS, 332, 915

Ibata, R., Lewis, G. F., Irwin, M., Totten, E. & Quinn, T. 2001, ApJ, 551, 294

Ibata, R. A., Wyse, R. F. G., GIlmore, G., Irwin, M. J. & Suntzeff, N. B. 1997, AJ, 113, 634 (I97)

Johnston, K. V., Hernquist, L. & Bolte, M. 1996, ApJ, 465, 278

Johnston, K. V., Majewski, S. R., Siegel, M. H., Reid, I. N. & Kunkel, W. E. 1999, AJ, 118, 1719

Johnston, K. V., Spergel, D. N. & Haydn, C. 2002, ApJ, 570, 656 (JSH)

King, I. R. 1966, AJ, 71, 64

Klypin, A., Kravtsov, A. V., Valenzuela, O. & Prada, F. 1999, ApJ, 522, 82

Law, D. R., Johnston, K. V. & Majewski, S. R. 2003b, in prep. (Paper III)

Law, D. R., Majewski, S. R., Skrutskie, M. F. & Johnston, K. V. 2003a, in ASP Conf. Ser. Vol., Tidal Tails and Galactic Satellites, eds. D. Martinez-Delgado & F. Prada, in press (astro-ph/0309567)

Majewski, S. R., Skrutskie, M. F., Weinberg, M. D. & Ostheimer, J. D. 2003a, ApJ, in press (astro-ph/0304198; MSWO)

Majewski, S. R., et al. 2003b, AJ, submitted (Paper II)

Miyamoto, M. & Nagai, R. 1975, PASJ, 27, 533

Moore, B., Ghigna, S., Governato, F., Lake, G., Quinn, T., Stadel, J. & Tozzi, P. 1999, ApJ, 524, L19

Navarro, J. F., Frenk, C. S. & White S. D. M. 1997, ApJ, 490, 493

Sakamoto, T., Chiba, M. & Beers, T. C. 2003, A&A, 397, 899

Sparke, L. S. 2002, in The Shapes of Galaxies and Their Dark Matter Halos, ed. P. Natarajan, (Singapore: World Scientific), p. 178