text:
# Grover’s Quantum Search Algorithm for an Arbitrary Initial Mixed State

## I Introduction

Grover’s search algorithm Grover96 ; Grover97  provides an example of the speed-up that would be offered by quantum computers, if and when they are built. The problem solved by Grover’s algorithm is finding a sought-after (“marked”) element in an unsorted database of size $`N`$. To solve this problem, a classical computer would need $`N/2`$ database queries on average, and in the worst case it would need $`N1`$ queries. Using Grover’s algorithm, a quantum computer can find the marked state using only $`O(\sqrt{N})`$ quantum database queries. The importance of Grover’s result stems from the fact that it proves the existence of a gap (albeit a polynomial gap) between the power of quantum computers and classical computers. Moreover, the algorithm may be used to speed up the solution of many problems (such as NP-complete problems), for which no efficient classical algorithms is known.

Along this paper we assume without loss of generality that $`N=2^n`$, where $`n`$ is an integer. The algorithm requires a register of $`n`$ qubits carrying the computation. When we say it is in a state $`|x`$, we mean that its qubits are in states corresponding to the binary representation of the number $`x`$. Grover’s original quantum search algorithm consists of the following steps:

1. Initialize the register to $`H|0`$. That is, reset all the qubits to 0 and apply the Hadamard transform to each of them.
2. Repeat the following operation (named the *Grover Iterate* $`Q`$) $`T=\frac{\pi \sqrt{N}}{4}`$ times:

1. Rotate the marked state $`|k`$ by a phase of $`\pi `$ radians $`(I_k^\pi )`$.
    2. Apply the Hadamard transform to the register.
    3. Rotate the $`|0`$ state by a phase of $`\pi `$ radians $`(I_0^\pi )`$.
    4. Apply the Hadamard transform again.
3. Measure the resulting state.

Several generalizations extended the original Grover algorithm. Among these is handling multiple marked states BBHT96 , and the initialization of the algorithm in any pure state BB99 . Another generalization is the replacement of the Hadamard transform by any other unitary operation Grover98a ; Ging00 ; BHMT00 . In this way the algorithm may be used to speed-up many classical decision algorithms and heuristics. Other generalizations use arbitrary rotation angles LLZT00 , replace the $`|0`$ state from Step (c) with any other state, or combine all of the above BB01 . The rotation angles may be tweaked in order to find a marked state with certainty Hoyer00 ; LLS01 .

The original Grover Iterate is $`Q=HI_0^\pi HI_k^\pi `$. It has been generalized to $`Q=UI_s^\beta U^{}I_M^\gamma `$, where $`U`$ is an arbitrary unitary operator, $`s`$ is an arbitrary state, $`\beta `$ and $`\gamma `$ are arbitrary angles, and $`M`$ includes any number of marked states. We now observe that any unitary operation $`Q`$ has a unitary diagonalization. Therefore, it can be represented as $`Q=UI_S^\stackrel{}{\beta }U^{}I_M^\gamma `$. This is a further generalization of Grover’s algorithm, where the state $`s`$ is replaced by a set of states $`S`$, each of which may have a different rotation angle. Thus, every iterative algorithm is a generalized Grover algorithm.

In this paper we study the case where the generalized Grover Iterate of BB01  is applied to a quantum register that is initialized in an arbitrary mixed state. Our study extends and corrects a result from Bose00 .

## II Arbitrary Pure Initial State

If the abovementioned search algorithm is used as a procedure by another algorithm, it might be necessary to avoid its first step. Even if the initialization is performed, gate imperfection or external noise might cause the outcome to differ from the exact $`H|0`$ state. Rather, it may well be some general pure state $`|\psi _0`$, which is a superposition of the marked state and the unmarked states. In addition, the Iterate itslef may be imperfect: the Hadamard operation might be some other unitary operation $`U`$; the rotations of (2a) and (2c) may be in angles $`\beta `$ and $`\gamma `$ (respectively), different of $`\pi `$; and the rotated state of Step (2c) may be a non-zero $`|s`$. Finally, the set of sought-after items, $`M`$, may include multiple items.

When the parameters of the problem are known, we may follow the results of Biham et al. BB01 , and calculate the probability to measure a marked state $`P_{\psi _0}`$ as a function of the number of Grover iterations $`t`$:

$$P_{\psi _0}(t)=P_{\psi _0}\mathrm{\Delta }P_{\psi _0}\mathrm{cos}(2\omega t+2\varphi _{\psi _0}).$$
(1)
$`P_{\psi _0}`$ and $`\mathrm{\Delta }P_{\psi _0}`$ denote the average over time and the amplitude of $`P_{\psi _0}`$ respectively. The subscripts $`\psi _0`$ denote that the values depend on the initial state. However, $`\omega `$, which is defined by

$$\mathrm{cos}\omega =\underset{iM}{}\left|i|U|s\right|^2\mathrm{cos}\frac{\beta +\gamma }{2}+\underset{iM}{}\left|i|U|s\right|^2\mathrm{cos}\frac{\beta \gamma }{2},$$

is independent of the initial state. In a large search problem with $`N\mathrm{}`$, and with the original Grover Iterate, $`\omega `$ may be approximated by $`\omega =\frac{2}{\sqrt{N}}`$. With the original initial state $`H|0`$ studied by Grover, we re-obtain $`P_{H|0}=\mathrm{\Delta }P_{H|0}=\frac{1}{2}`$ and $`\varphi _{H|0}0`$.

## III Arbitrary Mixed Initial State

A mixed state arises when one cannot describe the state of a quantum system deterministically, no matter what basis one chooses. Such a condition appears very often when a quantum system is entangled with its environment, while the environment cannot be accessed or manipulated. The state of such a system may be described by a completely-positive trace-1 hermitian density matrix, denoted by $`\rho `$. An equivalent description is an ensemble $`=\{p_\mu ,|\psi _\mu \}`$ where $`_\mu p_\mu =1`$ and $`\rho =_\mu p_\mu |\psi _\mu \psi _\mu |`$. According to this description, the system is in the pure state $`|\psi _\mu `$ with probability $`p_\mu `$. When a unitary operation $`V`$ is applied to the mixed state it transforms the state into $`_\mu p_\mu V|\psi _\mu \psi _\mu |V^{}`$. The mixedness of the state does not change, and it may be thought as if $`V`$ transforms each of the components of $``$ independently of the others.

Extending the argument of Section II, the initial state of the quantum register might not be pure, due to external noise, decoherence or previous manipulations. Instead, the initial state may be some general mixed state $``$. Given the description of $``$ as an ensemble, all we can say is that the register is in the pure state $`|\psi _\mu `$ with probability $`p_\mu `$ (for all $`i`$’s).

When the Grover algorithm is applied to a register which is in a pure state $`|\psi _\mu `$, the probability to measure the marked state is $`P_\mu (t)`$. The probability for the register to be in $`|\psi _\mu `$ is $`p_\mu `$ (Considering the ensemble $``$). Thus, the total probability to measure the marked state is the weighted average

$`\stackrel{~}{P}(t)`$ $`=`$ $`{\displaystyle \underset{\mu }{}}p_\mu P_\mu (t)`$ (2)
$`=`$ $`{\displaystyle \underset{\mu }{}}p_\mu \left(P_\mu \mathrm{\Delta }P_\mu \mathrm{cos}\left(2\omega t+2\varphi _\mu \right)\right).`$
The functions

$$P_\mu (t)P_\mu =\mathrm{\Delta }P_\mu \mathrm{cos}\left(2\omega t+2\varphi _\mu \right)$$

share a sinusoidal form, differing in amplitude and phase, but not in frequency. They may be thought of as the projections of vectors rotating in frequency $`\omega `$, as exemplified by Figure 1.

Therefore, their weighted sum $`\stackrel{}{\mathrm{\Delta }\stackrel{~}{P}}`$ (the Center of Mass of the vectors in the figure) is a sinusoidal function with the same frequency:

$$\stackrel{~}{P}(t)=\stackrel{~}{P}\stackrel{~}{\mathrm{\Delta }P}\mathrm{cos}\left(2\omega t+2\stackrel{~}{\varphi }\right)$$
(3)
while

$$\stackrel{~}{P}=\underset{\mu }{}p_\mu P_\mu ,$$

$$\stackrel{~}{\mathrm{\Delta }P}=\sqrt{\left(\underset{\mu }{}p_\mu \mathrm{\Delta }P_\mu \mathrm{cos}2\varphi _\mu \right)^2+\left(\underset{\mu }{}p_\mu \mathrm{\Delta }P_\mu \mathrm{sin}2\varphi _\mu \right)^2}$$
(4)
and

$$\mathrm{tan}2\stackrel{~}{\varphi }=\frac{_\mu p_\mu \mathrm{\Delta }P_\mu \mathrm{sin}\left(2\varphi _\mu \right)}{_\mu p_\mu \mathrm{\Delta }P_\mu \mathrm{cos}\left(2\varphi _\mu \right)}.$$
(5)
The probability to measure a marked state reaches its maximum value

$$\stackrel{~}{P}_{max}=\stackrel{~}{P}+\stackrel{~}{\mathrm{\Delta }P}$$

after $`T=\frac{\pi 2\stackrel{~}{\varphi }}{2\omega }`$ iterations.

If the algorithm is repeated until success with $`T`$ iterations each time, the expected total time to measure a marked state is

$$𝒯_Q=\frac{\pi 2\stackrel{~}{\varphi }}{2\omega \stackrel{~}{P}_{max}}$$

since the number of repetition until success is distributed geometrically with parameter $`\stackrel{~}{P}_{max}`$. When the original Iterate is used, and a single item is sought after, this reduces to $`𝒯_Q=\frac{\pi 2\stackrel{~}{\varphi }}{4\stackrel{~}{P}_{max}}\sqrt{N}.`$ If this value is significantly smaller than the classical expected time $`𝒯_C=N/2`$, then the quantum algorithm has an advantage. Quantitatively, the expected number of oracle queries that the quantum algorithm requires is smaller by a factor of

$$\frac{𝒯_C}{𝒯_Q}=\frac{N\omega \stackrel{~}{P}_{max}}{\pi 2\stackrel{~}{\varphi }}=\frac{2\stackrel{~}{P}_{max}\sqrt{N}}{\pi 2\stackrel{~}{\varphi }}.$$
(6)
## IV Examples

For clarity and simplicity, our examples use the original Grover Iterate and single marked state $`|k`$, with different initial mixed states.

### IV.1 Pure Initial State

When the arbitrary mixed state is chosen to be pure, the summations are degenerated and the results of BB99  are achieved. For example, if the initial state is the original $`=\{p=1,H|0\}`$, the original Grover case is found. If $`=\{p=1,|k\}`$, then $`\stackrel{~}{P}=\stackrel{~}{\mathrm{\Delta }P}=\frac{1}{2}`$ and $`\stackrel{~}{\varphi }=\frac{\pi }{2}`$. An interesting known property of the Grover algorithm is that for all states orthogonal to both $`|k`$ and $`H|0`$, $`\stackrel{~}{P}=\stackrel{~}{\mathrm{\Delta }P}=0`$.

### IV.2 Pseudo-Pure Initial State

Ensembles where a state $`|\psi `$ appears with probability $`ϵ+\frac{1ϵ}{N}`$ and any state orthogonal to it appears with equal probabilities of $`\frac{1ϵ}{N}`$ are called pseudo-pure mixed states. They are written more conveniently as $`\rho _{ϵ\text{-}pure}=(1ϵ)\frac{I}{N}+ϵ|\psi \psi |`$. Notice that $`0ϵ1`$ is a measure of the purity of $`\rho `$: when $`ϵ=0`$ it is totally mixed, and when $`ϵ=1`$ it is totally pure. It is easy to see that in the limit of large $`N`$, $`\stackrel{~}{P}=ϵP_\psi `$, $`\stackrel{~}{\mathrm{\Delta }P}=ϵ\mathrm{\Delta }P_\psi `$ and $`\stackrel{~}{\varphi }=\varphi _\psi `$. For example, for

$$\rho _{\frac{1}{\mathrm{log}N}\text{-}pure}=(1\frac{1}{\mathrm{log}N})\frac{I}{N}+\frac{1}{\mathrm{log}N}H|00|H,$$

we obtain $`\stackrel{~}{P}=\stackrel{~}{\mathrm{\Delta }P}=\frac{1}{2\mathrm{log}N}`$ and $`\stackrel{~}{\varphi }=0`$. Notice that although $`\rho `$ is extremely mixed, the quantum advantage is of factor $`\frac{2\sqrt{N}}{\pi \mathrm{log}N}`$.

### IV.3 Initial State Where $`m`$ of the Qubits Are Mixed

Let us study the case where the register is initialized to $`\rho _{m\text{-}mix}=2^m_{i=0}^{2^m1}H|ii|H`$. This state may occur if the $`m`$ least significant qubits of the register are totally mixed before the first Hadamard transform is applied. Since all $`H|i`$ are orthogonal to $`H|0`$ (except for $`H|0`$ itself) and they are almost orthogonal to $`|k`$ (since $`\left|k\left|H\right|i\right|^2=\frac{1}{N}`$), the evolution of $`\rho _{m\text{-}mix}`$ is governed by $`\{p=2^m,H|0\}`$ and we obtain $`\stackrel{~}{P}=\stackrel{~}{\mathrm{\Delta }P}=\frac{1}{2^{m+1}}`$ and $`\stackrel{~}{\varphi }=0`$. Large $`m`$ would render the algorithm useless.

## V Algorithm Usefulness and Entropy

The von Neumann entropy of a mixed state $`\rho `$ is defined as $`S(\rho )=\text{tr}\rho \mathrm{log}\rho `$. Bose et al. Bose00  presented a new model for quantum computation and laid out a new proof for the optimality of the Grover algorithm. However, one of their results was the following: if the Grover algorithm is initiated with a mixed state $`\rho `$, such that $`S(\rho )\frac{1}{2}\mathrm{log}N`$, the algorithm would have no advantage comparing to the classical case. This is in disagreement with our findings Bose00err . A counter-example to their claim is $`\rho _{\frac{1}{\mathrm{log}N}\text{-}pure}`$ as defined above. The entropy of pseudo-pure state is

$`S(\rho _{ϵ\text{-}pure})`$ $`=`$ $`S\left({\displaystyle \frac{1ϵ}{N}}I_N+ϵ|00|\right)`$
$`=`$ $`{\displaystyle \underset{1}{\overset{N1}{}}}{\displaystyle \frac{1ϵ}{N}}\mathrm{log}{\displaystyle \frac{1ϵ}{N}}`$
$`{\displaystyle \frac{1+(N1)ϵ}{N}}\mathrm{log}{\displaystyle \frac{1+(N1)ϵ}{N}}`$
$`=`$ $`(N1){\displaystyle \frac{1ϵ}{N}}\mathrm{log}{\displaystyle \frac{1ϵ}{N}}`$
$`{\displaystyle \frac{1+(N1)ϵ}{N}}\mathrm{log}{\displaystyle \frac{1+(N1)ϵ}{N}},`$
and for large $`N`$, where $`N/(N1)1`$, (7)
$``$ $`(1ϵ)\mathrm{log}{\displaystyle \frac{1ϵ}{N}}\left({\displaystyle \frac{1}{N}}+ϵ\right)\mathrm{log}\left({\displaystyle \frac{1}{N}}+ϵ\right)`$
$`=`$ $`(1ϵ)\mathrm{log}N(1ϵ)\mathrm{log}(1ϵ)`$
$`\left({\displaystyle \frac{1}{N}}+ϵ\right)\mathrm{log}\left({\displaystyle \frac{1}{N}}+ϵ\right)`$
$`=`$ $`(1ϵ)\mathrm{log}N\mathrm{},`$
where $`\mathrm{}=(1ϵ)\mathrm{log}(1ϵ)+\left(\frac{1}{N}+ϵ\right)\mathrm{log}\left(\frac{1}{N}+ϵ\right)(1,0.8)`$ for any $`0ϵ1`$ and any $`N2`$. For $`ϵ=\frac{1}{\mathrm{log}N}`$ we obtain $`S(\rho _{\frac{1}{\mathrm{log}N}\text{-}pure})=\left(1\frac{1}{\mathrm{log}N}\right)\mathrm{log}N+O(1)=\mathrm{log}N+O(1)`$. This entropy is almost maximal. However, as noted above, the Grover algorithm outperforms any classical algorithm, even when it is initialized with this state.

Entropy is not a good measure for the usefulness of Grover’s algorithm. For practically every value of entropy,there exist states that are good initializers and states that are not. For example, $`S(\rho _{(n1)\text{-}mix})=\mathrm{log}N1=S(\rho _{\frac{1}{\mathrm{log}N}\text{-}pure})`$, but when initialized in $`\rho _{(n1)\text{-}mix}`$, the Grover algorithm is as bad as guessing the marked state. Another example may be given using the pure states $`H|00|H`$ and $`H|11|H`$. With the first, Grover arrives to the marked state with quadratic speed-up, while the second state is practically unchanged by the algorithm.

## VI Acknowledgments

We thank Tal Mor for valuable discussions and suggestions that made the compiling of this paper possible. The work was partially supported by the European Commission through the IST Programme under contract IST-1999-11234. The first author was partially supported by the fund for the promotion of research at the Technion and the Israel MoD Research and Technology Unit.


format: text

text:
# 𝐾^∗(892) Production in Au+Au and pp Collisions at √𝑠_{𝑁𝑁}=200GeV at STAR

## 1 Introduction

Resonance particles, which have very short lifetime (few fm/c) comparable to the time scale for the evolution of the hot-dense matter formed in ultra-relativistic heavy-ion collisions, can be a very unique tool to probe the dynamics and properties of the high density matter . In particular, the $`K^{}`$(892) has a lifetime of $``$4fm/c, and can be produced at the chemical freeze-out stage. The $`K^{}`$ short lifetime makes it possible for the newly-formed $`K^{}`$ to undergo a period of re-interaction in the hadronic gas phase. A portion of $`K^{}`$ may decay before the kinetic freeze-out stage and their kaon and pion daughter particles might be re-scattered by other particles in the hadron gas. This effect of the $`K^{}`$ daughter particles re-scattering can destroy part of the overall $`K^{}`$ signal. On the other hand, kaon and pion particles in the hadron gas can regenerate $`K^{}`$ through the so-called pseudo-elastic collisions . This regeneration effect can compensate for the rescattering effect. Thus the measurement of $`K^{}`$ yields and their centrality dependence in heavy-ion collisions can provide information to estimate the time between chemical and kinetic freeze-out in relativistic heavy-ion collisions .

In the strongly interacting matter at high temperature and high densities, characteristics of the short-lived $`K^{}`$(892) resonance might be modified in the medium, with shifted mass, broadened width and even significantly changed line shapes. In the hadron gas, kaon and pion particles can regenerate $`K^{}`$ signals through $`K\pi K^{}K\pi `$. This regeneration channel can also interfere with the kaon and pion elastic scattering channel through $`K\pi K\pi `$. Thus the $`K^{}`$ meson may be modified not only due to this interference but also due to the kaon and pion initial phase space distributions . In addition, dynamical interactions of the $`K^{}`$ meson with the surrounding matter may also cause the modification of the $`K^{}`$ mass, width and line shape . Even though the size of the system formed in pp collisions is smaller than in Au+Au collisions, interactions that may modify the $`K^{}`$ resonance are also expected. Thus a measurement of $`K^{}`$(892) mass, width and line shape in Au+Au and pp collisions can provide very interesting information on possible in-medium effects.

## 2 Data Analysis and Results

Preliminary measurements on $`K^0`$(892)$`K\pi `$ and $`K^\pm `$(892)$`K_S^0\pi ^\pm `$ in Au+Au and pp collisions at $`\sqrt{s_{NN}}`$=200GeV are presented. These measurements were made at the Solenoidal Tracker at RHIC (STAR) with the main detector Time Projection Chamber (TPC). In Au+Au collisions, the minimum bias trigger was defined by coincidences between two Zero Degree Calorimeters (ZDC). A scintillator Central Trigger Barrel (CTB) was used to select central collision events. In pp collisions, the minimum bias trigger was defined using coincidences between two beam counters that measured the charged particles multiplicity near beam rapidity. After requiring the collision vertex to be within $`\pm `$50cm along the beam line, about 2M top 10$`\%`$ central triggered, 2M minimum bias triggered Au+Au collision events and 6M minimum bias triggered pp collision events were used in this analysis. The events from minimum bias Au+Au collisions were divided in four centrality bins from the most central to peripheral collisions: 0$`\%`$-10$`\%`$, 10$`\%`$-30$`\%`$, 30$`\%`$-50$`\%`$ and 50$`\%`$-80$`\%`$.

Through energy loss ($`dE/dx`$) in the TPC gas, charged pions and kaons were identified. In the case of $`K^0`$, charged kaons were selected by requiring their $`dE/dx`$ to be within two standard deviations (2$`\sigma `$) of the expected value. A looser $`dE/dx`$ cut of 3$`\sigma `$ was used for charged pions. In the case of $`K^\pm `$, $`K_S^0`$ candidates were selected from their decay vertex geometries via $`K_S^0\pi ^+\pi ^{}`$. The invariant mass was then calculated for each kaon and pion pair in an event. The invariant mass distribution derived in this manner was then compared to a reference distribution calculated using uncorrelated kaons and pions from different events. In this analysis, $`K^0`$ and $`\overline{K^0}`$ are added together due to limited statistics. The term $`K^0`$ refers to the average of $`K^0`$ and $`\overline{K^0}`$ unless specified otherwise.

Figure 1 shows the $`K\pi `$ invariant mass distribution after subtraction of the mixed-event reference distribution in Au+Au and pp collisions. The invariant mass distributions are fit to a combination of a p-wave Breit-Wigner function and a linear function representing the background. In Figure 1 (a), there is a large amount of residual background in the invariant mass distribution after subtracting the reference distribution. Reference  discussed possible sources of this residual background. In addition, while mixing Au+Au collision events, different events will have different reaction planes so that the reference invariant mass distribution might have slightly different shapes from the same event distribution. After subtracting the reference distribution, this slight difference will appear as part of the residual background . In Figure 1 (b), in order to precisely measure the $`K^0`$ mass and width as a function of transverse momentum in pp collisions, we only selected the kaon candidates with momentum between 0.2 GeV/c and 0.7 GeV/c to minimize the residual background. In Figure 1 (c) and (d), $`K_S^0`$ candidates were first reconstructed via decay vertex geometries requiring $`|M_{K_S^0}M_{\pi ^+\pi ^{}}|<`$ 15MeV/c2. Then, after pairing $`K_S^0`$ candidates with charged pion candidates in same events and mixed events, invariant mass distribution for the $`K^\pm `$ were obtained after reference distribution subtraction.

The $`K^0`$ mass and width as a function of transverse momentum from Au+Au and pp collisions are shown in Figure 2. In pp collisions, the $`K^0`$ masses in low $`p_T`$ region are systematically smaller than the Monte Carlo (MC) results which account for detector effects and all kinematic cuts. The mass shift is $`p_T`$ dependent, with the $`K^0`$ mass increasing as a function of $`p_T`$. In Au+Au collisions, the $`K^0`$ mass shift is also observed at low $`p_T`$. In the case of the $`K^0`$ width, there is no significant difference between the measured results and the MC results in both Au+Au and pp collisions. As already discussed in the introduction section, there are various physical effects which may cause this $`K^{}`$ mass shift. However, a soft $`K^{}`$ resonance (low $`p_T`$) is more likely to be modified in the medium and thus a larger mass shift is expected for soft $`K^{}`$ than for hard $`K^{}`$ (high $`p_T`$). A similar mass shift for the $`\rho ^0`$(770) meson has also been observed in both pp and Au+Au collisions.

The detector acceptance and efficiency corrected $`K^0`$ transverse mass ($`m_T=\sqrt{p_T^2+m^2}`$) spectra at mid-rapidity ($`|y|<`$0.5) are shown in the left plot of Figure 3. The spectra are fit with exponential functions and $`K^0`$ yields $`dN/dy`$ and inverse slopes are extracted from the exponential fit. The $`K^0`$ $`dN/dy`$ increases from pp collisions to peripheral Au+Au collisions and to central Au+Au collisions. The inverse slopes in Au+Au collisions are systematically larger than that in pp collisions. In the right plot of Figure 3, the $`K^0`$ mean $`p_T`$ as a function of number of charged hadrons is compared to that of $`\pi ^{}`$, $`K^{}`$ and $`\overline{p}`$ . In pp collisions, the $`K^0`$ mean $`p_T`$ is comparable to the mean $`p_T`$ of $`\overline{p}`$. In Au+Au collisions, the $`K^0`$ mean $`p_T`$ are systematically larger than that in pp collisions. The larger $`K^{}`$ inverse slopes and mean $`p_T`$ in Au+Au collisions than in pp collisions might be explained by the fact that in the hadronic phase of Au+Au collisions, $`K^{}`$ with higher $`p_T`$ are more likely to escape and decay outside the fireball thus avoiding the daughter particles’ re-scattering. A low $`p_T`$ $`K^{}`$ has more chances to be destroyed by daughter particles’ re-scattering in the hadron gas medium in Au+Au collisions.

$`K^{}`$ and $`K`$ mesons of the same charge have identical quark content and only differ in mass and spin so that the $`K^{}/K`$ ratio might reflect the fireball evolution conditions from chemical to kinetic freeze-out. Figure 4 shows that the $`K^{}/K`$ ratios in Au+Au collisions are significantly reduced compared to the ratio in pp collisions. This ratio reduction may indicate that in Au+Au collisions, $`K^{}`$ daughter particles’ re-scattering might destroy more $`K^{}`$ than the amount of $`K^{}`$ which are regenerated by the $`K\pi `$ interactions in the hadron gas. The upper right panel of Figure 4 shows the $`K^{}/K`$ ratio as a function of the colliding energy from different collision systems. The $`K^{}/K`$ ratio from elementary collisions does not change noticeablly in the presented $`\sqrt{s}`$ region. However, the ratio in Au+Au collisions is significantly smaller. The $`K^{}`$ and $`\varphi `$ mesons have a small mass difference and $`\mathrm{\Delta }S`$=1. Thus, the $`\varphi /K^{}`$ ratio might also be a good signature to study the strangeness enhancement effect in ultra-relativistic heavy ion collisions. The bottom right panel of Figure 4 shows that the $`\varphi /K^{}`$ ratio increases as a function of the colliding energy. We may need more statistics to identify any possible $`\varphi /K^{}`$ ratio differences between Au+Au collisions and pp collisions at $`\sqrt{s_{NN}}`$=200GeV. In this analysis, $`K`$ results are from  and $`\varphi `$ results are from , $`K^{}`$ results for Au+Au collisions at $`\sqrt{s_{NN}}`$=130GeV are from .

## 3 Summary

In summary, we have measured the $`K^0`$(892) and $`K^\pm `$(892) production from Au+Au and pp collisions at $`\sqrt{s_{NN}}`$=200GeV at STAR. Possible in-medium dynamical effects might have modified the $`K^{}`$ mass line shape and thus systematically smaller $`K^{}`$ masses are observed at the low $`p_T`$ region in pp and Au+Au collisions. In Au+Au collisions, the $`K^{}`$ daughter particles’ re-scattering effect is dominant over the regeneration effect. This picture is consistent with our measurements that the $`K^{}`$ inverse slopes and mean $`p_T`$ in Au+Au are larger than in pp and the $`K^{}/K`$ ratios in Au+Au are smaller than in pp.


text:
# Nonlinear optics via double dark resonances

## I Introduction

Dark resonances arise from the quantum superposition states that are decoupled from coherent and dissipative interactions. They are now a well-known concept in quantum optics and laser spectroscopy, leading to concepts such as electromagnetically induced transparency (EIT) Harris (1997), lasing without inversion and increased resonant index of refraction Scully and Zubairy (1997), highly efficient resonant nonlinear optics Merriam et al. (2000); Lukin and Imamoğlu (2000); Hemmer et al. (1995), “slow” light Hau et al. (1999); Kash et al. (1999); Budker et al. (1999); Kocharovskaya et al. (2001), coherent storage of photon states Liu et al. (2001); Phillips et al. (2001); Walsworth et al. (2002), and femto-second generation Harris and Sokolov (1998); Hakuta et al. (1992).

As a rule, any perturbation of the dark state leads to undesirable decoherence. However, it was shown theoretically that “perturbing” the dark state coherently leads to multiple quantum superposition states that interact coherently. These so-called “double-dark” states Lukin et al. (1999); Ye et al. (2002) can be used to mitigate decoherence effects and thus enlarge the domain of dark state physics. The perturbation can come from different sources: microwave or optical coherent fields, DC-fields, non-adiabatic optical coupling (i.e. effective fields). In this paper we discuss theoretically and experimentally the use of the magnetic component of a radio frequency (RF) field. We observe two sharp resonances at shifted frequencies. We also demonstrate that the double dark resonances allow to enhance the nonlinear optical process. In particular, we demonstrate the efficient generation of a comb of multiple frequency components associated with double dark resonances.

## II Theory

### II.1 Dark resonances and electromagnetically induced transparency

The principle of electromagnetically induced transparency (EIT) is the destructive interference of two absorption channels, which originate from a Stark split resonance. It is based on a so-called dark resonance:

For example, in a $`\mathrm{\Lambda }`$-system like the one in Fig. 1a with the excited state $`|a`$ coupled to the metastable states $`|b`$ and $`|c`$ via strong driving and weak probe fields with Rabi frequencies $`\mathrm{\Omega }`$ and $`\alpha `$, respectively, the resonant interaction Hamiltonian is $`H=\mathrm{}\mathrm{\Omega }|ca\left|+\mathrm{}\alpha \right|ba|+\text{H.c}.`$. The superposition state $`|+\mathrm{\Omega }|c+\alpha |b`$ is coupled by both fields, $`H=\mathrm{}\mathrm{\Omega }_+|+a|+\text{H.c.}`$, with $`\mathrm{\Omega }_+^2=\mathrm{\Omega }^2+\alpha ^2`$, whereas the orthogonal state $`|\mathrm{\Omega }|b\alpha |c`$ is dark, and $`H|=0`$. This is true as long as the probe and driving fields are in two-photon resonance with the transition $`|b|c`$. If all the atoms in the medium are in the dark state the medium is transparent to both fields (see probe spectrum in Fig. 2 (broken line)).

### II.2 Double-dark resonances

The usual mechanism by which perturbing the dark state destroys the transparency has been reviewed recently Lukin et al. (1999): Here we show that a coherent perturbation, e.g., a weak monochromatic electromagnetic field with Rabi frequency $`\mathrm{\Omega }_c`$ coupling an additional state (see Fig.1b), can display interesting and highly coherent features. This case is depicted in Fig. 2 (solid line), where a sharp absorption line appears on resonance, that splits the point of transparency in two. The two transmission maxima are now detuned by $`\pm \mathrm{\Omega }_c`$. This absorption line can be very narrow, as its width scales with the ratio of the intensities $`\mathrm{\Omega }_c^2/\mathrm{\Omega }^2`$. Nevertheless, its depth shows that the resonant absorption is as strong as that of the original Stark-split levels. For a simplified explanation of this phenomenon, the dressed-state picture is very useful, where the Hamiltonian is diagonalized with respect to driving and perturbing fields. The resulting resonant state, $`|0|d+\mathrm{\Omega }_c/\mathrm{\Omega }|a`$, has a small admixture of the excited state, and thus it decays very slowly, which results in a very narrow resonance. In this picture, it is easy to prove that in the all-resonant case, detuning the field $`\mathrm{\Omega }_c`$ by $`\mathrm{\Delta }_{\mathrm{RF}}`$ implies a frequency shift of the novel resonance line by $`\mathrm{\Delta }_{\mathrm{RF}}`$.

### II.3 Model

The experiment described in this article follows a slightly different setup: The states $`|c`$ and $`|d`$ are degenerate (i.e., Zeeman sublevels), whereas the new field, a radio frequency (RF) field, has the fixed frequency $`\nu _{RF}`$. This situation is the same as if there were two novel fields that are blue and red detuned by $`\pm \mathrm{\Delta }_{\mathrm{RF}}=\nu _{\mathrm{RF}}`$. One can thus expect two novel lines with half the absorption depth, shifted by $`\pm \nu _{RF}`$ compared to the center, as in Fig. 3a (inset).

In order to understand this behavior, we investigate briefly the dressed state system (i.e., the eigenvalues/eigenstates of the atom + drive + RF - system)(cf. Fig. 4): The Hamiltonian of this setup (in the interaction picture without the probe field) is

$`H`$ $`=`$ $`\mathrm{}\mathrm{\Omega }|ac|+{\displaystyle \frac{\mathrm{}}{2}}\mathrm{\Omega }_c\left(|d_1c|+|d_2c|\right)+H.c.`$ (1)
$`+\mathrm{}\mathrm{\Delta }_{\mathrm{RF}}\left(|d_1d_1||d_2d_2|\right).`$
Here, $`|d`$ is split into $`|d_1`$ and $`|d_2`$ to include the fact that the RF field is red and blue detuned at the same time. From this Hamiltonian we find the dressed state system

$`|\pm `$ $``$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left(|a\pm |c\right)`$
with
$`\lambda _\pm `$ $``$ $`\pm \mathrm{\Omega }`$
and
$`|0\pm `$ $``$ $`|d_{1/2}+{\displaystyle \frac{\mathrm{\Omega }_c}{2\mathrm{\Omega }}}f(\mathrm{\Omega },\mathrm{\Delta }_{\mathrm{RF}})\left(\pm |a+{\displaystyle \frac{\mathrm{\Delta }_{\mathrm{RF}}}{\mathrm{\Omega }}}|c\right)`$
with
$`\lambda _{0\pm }`$ $``$ $`\pm \mathrm{\Delta }_{\mathrm{RF}}`$
Here we use $`f(\mathrm{\Omega },\mathrm{\Delta }_{\mathrm{RF}})=\mathrm{\Omega }^2/(\mathrm{\Omega }^2\mathrm{\Delta }_{\mathrm{RF}}^2)`$ and assume $`\mathrm{\Omega }\mathrm{\Delta }_{\mathrm{RF}}`$. The small admixture of the excited state makes sure that the dressed states $`|0\pm `$, replacing $`|d`$, decay with a very narrow linewidth, and the eigenenergies coincide directly with the frequencies of the absorption lines.

### II.4 Nonlinear sideband generation

This arrangement also produces efficient nonlinear modulation: sidebands appear in even multiples of the RF frequency $`\nu _{\mathrm{RF}}`$. This can be explained the following way: Two RF photons couple (via state $`|d`$) to virtual states $`|c_{2\pm }`$, that are lying exactly $`2\nu _{RF}`$ above and below the level of $`|c`$ (see Fig. 5). The driving field parametrically creates a transparency for the probe field at $`2\nu _{RF}`$ above and below its original frequency. Higher order nonlinearities then can be seen for fourth, sixth, etc., order virtual levels $`|c_{4\pm }`$, $`|c_{6\pm }`$, etc., with decreasing intensities (see section III.3).

### II.5 Doppler narrowing

The experiment is performed in a hot gas and therefore Doppler broadening is an important consideration. We will give here a brief explanation of the impact of Doppler shifts, as they pose an interesting puzzle. For each velocity component of the atomic gas, the probe and drive transitions are affected by Doppler shifts equally, since these shifts depend linearly on the frequency and probe and drive fields copropagate. The RF transition, however, is essentially unaffected because of its very low frequency. It can be expected that the transparency on resonance in the Doppler-broadened case gives rise to a transparency whose width decreases for increasing Doppler broadening. Surprisingly, there are Doppler narrowed transmission lines replacing the absorption lines of the Doppler free case. This “Doppler narrowing” can be explained using Fig. 6.

For each velocity group there are 4 eigenvalues in the system. Each eigenvalue is responsible for an absorption line. In the limit of vanishing RF field, and with the probe and drive fields on resonance for the zero velocity group $`\mathrm{\Delta }_{\mathrm{Dop}}=0`$, they are

$`\lambda _{0\pm }^{\mathrm{Dop}}`$ $`=`$ $`\pm \mathrm{\Delta }_{\mathrm{RF}},`$
$`\lambda _\pm ^{\mathrm{Dop}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }_{\mathrm{Dop}}}{2}}\pm \sqrt{{\displaystyle \frac{\mathrm{\Delta }_{\mathrm{Dop}}^2}{4}}+\mathrm{\Omega }^2},`$
where $`\mathrm{\Delta }_{\mathrm{Dop}}`$ is the Doppler frequency shift. Shown in Fig. 6 are the energy eigenvalues $`\lambda _\pm ^{\mathrm{Dop}}`$ and $`\lambda _{0\pm }^{\mathrm{Dop}}`$ of the atomic dressed states as functions of the Doppler shift $`\mathrm{\Delta }_{\mathrm{Dop}}`$. Eigenvalues $`\lambda _\pm ^{\mathrm{Dop}}`$ (shown in the figure as broken lines) are responsible for the Stark split absorption lines in usual EIT. The dotted lines show all eigenvalues without both, driving and RF fields ($`\mathrm{\Omega }=\mathrm{\Omega }_c=0`$. Note that in this case there are no avoided crossings.) With both fields present all anticrossings appear in Fig. 6 (solid lines). Performing the Doppler average over all components (which in mathematical terms means a Gaussian average over all $`\mathrm{\Delta }_{\mathrm{Dop}}`$) leads to an absorption profile where narrow transparency lines can be seen very close to the $`\pm \mathrm{\Delta }_{\mathrm{RF}}`$ eigenvalues. The reason is obvious from Fig. 6, where the regions very close to $`\pm \mathrm{\Delta }_{\mathrm{RF}}`$ and to 0 are the only ones that are not crossed by any of the eigenvalues over the width of all $`\mathrm{\Delta }_{\mathrm{Dop}}`$. Thus these lines are nearly at the frequency where there is maximum absorption in the Doppler free case, and become narrower the broader the Doppler shift is spread.

The resulting transmission curve can be seen in the simulation of Fig. 3, where the three lines, (one EIT and two double dark lines) are clearly visible. The parameters in (b) are chosen to coincide with the experimental reality: The Doppler broadening in our experiment is 530 MHz.

The sideband generation is much more pronounced in the Doppler broadened medium than in the non-broadened one, which seems surprising at first. It helps to remember that the RF novel lines are absorption resonances in the non-broadened case, whereas they are transmission resonances in the Doppler broadened case. These transmission lines coincide with the zero-order dark resonance lines of the sideband spectra (see, e.g., Fig. 9). The consequence is that the ratio of neighboring sideband intensities $`I_{4\mathrm{\Delta }_{\mathrm{RF}}}/I_{2\mathrm{\Delta }_{\mathrm{RF}}}I_{2\mathrm{\Delta }_{\mathrm{RF}}}/I_0`$ is roughly two orders of magnitude higher for Doppler broadened media, such that the three-peak structure of the fundamental mode can also be seen in the second order sidebands (Fig. 9(a,b)).

## III Experiment

### III.1 Setup

The experiment was performed on the D1 resonance line of 87Rb (nuclear spin I=3/2, wavelength 795nm). The experimental setup is shown in Fig. 7. In our experiment, external cavity diode lasers (ECDL)

were used for both the drive and the probe fields. The two collinear beams, both of which are circularly polarized and which have the same diameter of 0.1 cm are overlapped and directed into a glass cell of length $`l=4`$ cm, filled with rubidium vapor. The vapor consisted of $`4\times 10^{11}`$cm-3 ($`T=65^o`$C) of the natural abundance of 85Rb and 87Rb, plus 3 Torr of Ne buffer gas. The buffer gas is needed in order to reduce the relaxation rate of the ground state coherence, which for vapor densities of this order stems mostly from the Rb atoms entering into and leaving the interaction region with the laser beams. Thus the ground state decoherence rate depends on the diameter of the laser beams and in our case can be estimated to be of the order of 5 kHz. The cell was mounted inside a three-layered magnetic shield. Perpendicular to the laser beams the RF field could be applied, using transverse coils. The amplitude of the magnetic component of the RF field is chosen such that it is of the order 0.1 G. The cell was heated using DC current in a high-resistance wire wound on the outside of the middle shield. The atomic density can then be estimated from the temperature in the cell using well-known vapor-pressure curves, e.g. from Nemeyanov (1963). The drive laser was tuned to the $`5S_{1/2}(F=2)5P_{1/2}(F=2)`$ transition; the probe laser was tuned to the $`5S_{1/2}(F=1)5P_{1/2}(F=2)`$ transition. This configuration of drive and probe fields creates an appropriate $`\mathrm{\Lambda }`$-system within the 87Rb D1 manifold (see Fig. 1c).

There must be minimal relative laser jitter between probe and drive fields, so the probe laser was phase-locked to the driving field frequency with an offset frequency of 6.8 GHz, which is equal to the hyperfine ground state splitting between the $`5S_{1/2}(F=1)`$ and $`5S_{1/2}(F=2)`$ states. This offset is fixed by a tunable microwave frequency synthesizer. After passing through the cell the laser beams hit a fast photodiode. The beat signal from this photodiode is monitored by a microwave spectrum analyzer, which is used as a broadband filter-amplifier with a bandwidth of 3 MHz, centered at the hyperfine splitting frequency.

### III.2 Spectrum

The lineshape of the dark resonances are then determined by scanning the probe laser frequency, which was accomplished by tuning the frequency of the synthesizer. While the laser power of the driving field can be changed from 0 to 6 mW, the probe laser was held at 1% of the drive power. With these parameters a broadening of the dark resonances could be observed with a slope of 30 kHz/mW. Field broadening and density narrowing of the dark resonance were investigated in Ref. Sautenkov et al. (1999). In order to observe EIT resonance and pronounced double dark structure, the driving field had to be strong enough to make the the additional double dark structures lie well within the EIT resonance. (As was tested experimentally in this scope, the novel features are not visible otherwise.) The optimal working point thus was determined to be 2.2 mW for the drive laser.

The double dark resonance structures can be seen in Fig. 8a. These additional resonances are induced by the alternating magnetic field with a frequency of 100 kHz. An RF field with a frequency of 350 kHz is shown in Fig. 8b. The new transparency lines are well outside the EIT peak. As can be seen in the figures, the frequency interval between the dark resonance and the new resonances indeed equals the RF frequency. The width of the novel resonances is 20 kHz, which suggests a total decoherence rate between the three ground states $`|b`$, $`|c`$, and $`|d`$ of 20 kHz.

### III.3 Sideband generation

The spectral distribution of the transmitted laser light was also studied. For this purpose the spectral analyzer was switched to a narrow band regime (bandwidth 30 kHz) and the central frequency of the bandwidth was scanned. The laser frequency during the spectrum analyzer scan was fixed. In this procedure several new optical components (sidebands) can be observed, which are separated by frequency intervals equal to harmonics of the applied magnetic field frequency. Obviously these optical components are generated during the propagation of laser light through the coherently prepared Rb vapor. In order to record their dependence on the laser frequency, we present snapshots where frequency steps of 20 kHz are taken, see Fig. 9.

The fundamental frequency here is the frequency of the beat-note of the drive- and probe- laser beams without the cell. The signal at the fundamental frequency is similar to the broadband detected signal in Fig. 8b, as could be expected. One can see the different behavior of the odd and even harmonics (sidebands). The amplitude of the second order sidebands that are observed at $`\pm `$700 kHz is higher than the the amplitude of the first order sidebands that are observed at $`\pm `$ 350 kHz. The reason that they are seen at all stems from the fact that the fundamental mode is broad enough to seep into the first order sidebands. In this case, the virtual states $`|d_\pm `$, detuned by $`\pm \mathrm{\Delta }_{\mathrm{RF}}`$ from state $`|d`$, are directly coupled to $`|a`$, which leads to first order sidebands. Note that in the simulations (Fig. 3) no first order sidebands are present.

The three-peak structure in the main diagonal of the lower curves of Fig. 9(c) is the main contribution to the integrated spectrum (upper curve), and for low RF intensities is all that is seen (e.g., in (a,b)). In reality, however, this contribution makes up only the fundamental mode. Sidebands originate in the resonances with $`|c_{2\pm }`$, $`|c_{4\pm }`$, etc., and can be seen for stronger RF intensity, such as in Fig. 9(c). There, each of the diagonals, which are offset from their neighbors by $`\pm 2\nu _{\mathrm{RF}}`$, consists of the fundamental three-peak structure. Finally, in the integrated picture they add up to the multi-peak spectrum in the upper curve.

## IV Conclusion

In conclusion, we have experimentally demonstrated novel coherent resonances as predicted in Ref. Lukin et al. (1999). These novel resonances are very sharp and highly sensitive to frequency and intensity of a RF field, which perturb the dark resonance. In addition, because of higher-order coupling, efficient non-linear sideband generation was demonstrated. This scheme can be used in applications like sensitive nonlinear measurements such as Ref. Yelin and Hemmer (2002) and UV laser generation Fry et al. (2000).

We want to thank M. O. Scully for his interest in this project and discussions leading to it, and we want to thank the Office of Naval Research and the National Science Foundation for support. SFY also acknowledges a Feodor-Lynen Fellowship from the Humboldt Foundation.


text:
# Chapter 1 Dissipation of Hydromagnetic Waves on Energetic Particles: Impact on Interstellar Turbulence and Cosmic Ray Transport

## Chapter 1 Dissipation of Hydromagnetic Waves on Energetic Particles: Impact on Interstellar Turbulence and Cosmic Ray Transport

### Abstract

The diffusion of galactic cosmic rays (CR) is considered. It is anticipated that the nonlinear MHD cascade sets the power-law spectrum of turbulence from the principle scale $`100`$ pc to much smaller scales. We found that the dissipation of waves due to the resonant interaction with energetic particles may terminate the cascade at less than $`10^{13}`$ cm. The shape of CR diffusion coefficient that was found may explain the observed peaks of secondary to primary nuclei ratios at a few GeV/n.

### 1. Introduction

The galactic CR have high energy density and they can not always be treated as test particles moving in given magnetic fields. In particular, the stochastic acceleration of CR by MHD waves is accompanied by the damping of waves. The wave damping causes the change of the wave spectrum that in turn affects the particle transport. Thus in principle the study of CR diffusion may need a selfconsistent consideration. We shall see that the CR action on interstellar turbulence should be taken into account at energies below $`10`$ GeV/n. The implementation of this effect in a full scale numerical simulations of CR propagation in the Galaxy was fulfilled by Moskalenko et al. (this Conference).

### 2. Equations for Cosmic Rays and Interstellar Turbulence

The steady state transport equation that describes the CR diffusion in the interstellar medium is of the form, e.g. :

$$D\mathrm{\Psi }\frac{}{p}p^2K\frac{}{p}p^2\mathrm{\Psi }=q,$$
(1)
where $`\mathrm{\Psi }(p,𝐫,t)`$ is the particle distribution function on momentum normalized on CR number density as $`N_{cr}(p)=_p^{\mathrm{}}𝑑p\mathrm{\Psi }`$, $`D(p,𝐫)`$ is the spatial diffusion coefficient, $`K(p,𝐫)`$ is the diffusion coefficient on momentum, and $`q(p,𝐫)`$ is the source term. If needed, the supplementary terms which describe the particle energy losses, nuclear fragmentation, and radioactive decays can be added to eq. (1). It is useful to introduce the diffusion mean free path $`l=3D/v`$.

It is assumed that the particle diffusion is due to the resonant wave-particle interaction. The diffusion mean free path is determined by the equation : $`l=r_gB^2\left(4\pi k_{res}W(k_{res})\right)^1`$, where $`v`$ is the particle velocity, $`r_g=pc/(ZeB)`$ is the particle Larmor radius in the average magnetic field $`B`$, $`k_{res}=1/r_g`$ is the resonant wave number, and $`W(k)`$ is the spectral energy density of waves defined as $`𝑑kW(k)=\delta B^2/4\pi `$ ($`\delta B`$ is the random magnetic field, $`\delta BB`$). The CR diffusion coefficient has the scaling $`Dv(p/Z)^a`$ for the power law spectrum $`W(k)1/k^{2a}`$. The diffusion on momentum is roughly described by the formula $`K=p^2V_a^2\left(a(4a)(4a^2)D\right)^1`$, where $`V_a`$ is the Alfven velocity.

In spite of the great progress in magnetic hydrodynamics, we do not yet have the well developed theoretical description of interstellar turbulence that would allow one to calculate $`W(k)`$ in different astrophysical conditions, see review . It is quite possible that there are two almost independent nonlinear cascades of waves in the magnetized plasma where the thermal pressure approximately equals the magnetic field pressure, e.g. . The cascade of Alfven waves (and the slow magnetosonic waves) leads to the Kolmogorov spectrum $`W(k)k^{5/3}`$, and the cascade of fast magnetosonic waves leads to the Kraichnan spectrum $`W(k)k^{3/2}`$ in the inertial range of wave numbers where the dissipation is absent. Below we consider the interaction of CR with these two cascades separately.

A. Kolmogorov type cascade. In its simplified form, the steady state equations for waves with a nonlinear transfer in $`k`$-space can be written as

$$\frac{}{k}\left(C_Ak^2\sqrt{kW(k)\left(4\pi \rho \right)^1}W(k)\right)=2\mathrm{\Gamma }_{cr}(k)W(k)+S_A\delta (kk_L).$$
(2)
$`kk_L`$, $`\rho `$ is the gas density (see , the original theory of Kolmogorov (1941) was developed for the incompressible liquid without magnetic field). The l.h.s. of eq. (3) describes the nonlinear cascade from small $`k`$ to large $`k`$, $`C_A`$ is a constant, and approximately equal to 0.3 according to the simulations . The r.h.s. of eq. (3) includes the wave damping on CR and the source term, which works on the main scale $`1/k_L=100`$ pc and describes the generation of turbulence by supernova bursts, stellar winds, and superbubbles expansion. In the limit of negligible damping $`\mathrm{\Gamma }_{cr}=0`$, the solution of eq. (2) gives the Kolmogorov scaling $`W(k)k^{5/3}`$.

The equation for wave amplitude attenuation on CR is :

$$\mathrm{\Gamma }_{cr}(k)=\pi e^2V_a^2\left(2kc^2\right)^1_{p_{res(k)}}^{\mathrm{}}𝑑pp^1\mathrm{\Psi }(p),$$
(3)
where $`p_{res}(k)=ZeB/ck`$. The solution of eqs (2), (3) allows finding the wave spectrum and the determination of the mean free path, which is

$$l(p)=l_{Ko}(p)\left[12\pi ^{3/2}V_ap^{1/3}l_{Ko}^{1/2}(p)\left(3C_AB^2r_g^{1/2}\right)^1_p^{p_L}𝑑p_2p_2^{2/3}_{p_2}^{\mathrm{}}\frac{dp_1}{p_1}\mathrm{\Psi }(p_1)\right]^2.$$
(4)
Here $`l_{Ko}`$ is the diffusion mean free path calculated for a Kolmogorov spectrum without regard of wave damping, and $`p_L=p_{res}(k_L)`$. The second term in brackets describes the modification of the mean free path due to the damping of short waves.

B. Iroshnikov-Kraichnan cascade. The simplified equation for waves reads simular to eq. (3) but with the different l.h.s.: $`\frac{}{k}\left(C_Mk^3\left(\rho V_a\right)^1W^2(k)\right)`$, where approximately $`C_M=1`$. At $`\mathrm{\Gamma }_{cr}=0`$, this gives the spectrum $`W(k)k^{3/2}`$ first found in . Using the same procedure as in the case A, one can obtain:

$$l(p)=l_{Kr}(p)\left[1\pi V_ap^{1/2}l_{Kr}(p)\left(2C_MB^2r_g\right)^1_p^{p_L}𝑑p_2p_2^{1/2}_{p_2}^{\mathrm{}}\frac{dp_1}{p_1}\mathrm{\Psi }(p_1)\right]^1.$$
(5)
Here $`l_{Kr}`$ is the diffusion mean free path calculated for a power law Kraichnan spectrum without regard of wave damping on CR.

As the most abundant species, the CR protons mainly determines the wave dissipation. Their distribution function $`\mathrm{\Psi }(p)`$ should be used to calculate $`l(p)`$ by the simultaneous solution of eqs (1) and (4) in the case A, and eqs (1) and (5) in the case B. The diffusion mean free path for other nuclei is $`l(p/Z)`$.

Let us estimate the effect of wave damping at $`1`$ GeV where approximately $`l=1`$ pc. The CR energy density is $`1`$ eV/cm3, $`V_a=10`$ km/s, $`B=3`$ $`\mu `$G. The second term in brackets in eq.(4) equals $`5\times 10^2`$, whereas the second term in brackets in eq.(5) equals $`10`$. We conclude that the Kolmogorov type cascade is not much affected by the damping on CR. The Iroshnikov-Kraichnan cascade is significantly affected, and this should lead to the modification of CR transport at energies less than about $`10`$ GeV.

### 3. Simple Selfconsistent Solution

To demonstrate the effect of wave damping, we consider a simple case of one-dimensional diffusion with the source distribution $`q=q_0(p)\delta (z)`$, $`q_0p^{\gamma _S}`$ (that corresponds to the infinitely thin disk of CR sources located at $`z=0`$) and the flat CR halo of height $`H`$, see . Let us assume that stochastic reacceleration does not essentially change the CR spectrum during the time of CR exit from the Galaxy, i.e. one can set $`K=0`$ in eq. (1). The solution of eq. (1) in the galactic disk is then $`\mathrm{\Psi }(p)=3q_0(p)H\left(2vl(p)\right)^1`$. We consider the Iroshnikov-Kraichnan cascade and simplify eq. (5) using the approximation $`_{p_2}^{\mathrm{}}𝑑p_1p_1^1\mathrm{\Psi }(p_1)=\mathrm{\Psi }(p_2)/(\gamma _s+0.5)`$. This allows one to find the self consistent diffusion mean free path for protons:

$$l(p)=l_{Kr}(p)\mathrm{exp}\left(3\sqrt{\pi }eH\left(8(\gamma _s+0.5)C_M\sqrt{\rho }c\right)^1_p^{p_L}𝑑p_1q_0(p_1)v^1(p_1)\right).$$
(6)
The form of eq. (6) is close to that needed to explain the peaks in ratios of secondary to primary nuclei in CR at a few GeV/n. Note, that the escape length (the grammage) that determines the production of secondaries is equal to $`X=3\mu H/(2l)`$, where $`\mu `$ is the surface gas density of the galactic disk .

### 4. Conclusion

The damping on CR terminates slow Iroshnikov-Kraichnan cascade but probably has no impact on the Kolmogorov-type cascade in the interstellar medium. The estimates were made for the turbulence level that provides the empirical value of CR diffusion coefficient. This finding offers a new explanation of the peaks in the secondary/primary ratios at a few GeV/n: the amplitude of short waves is small because of the damping and thus the low energy particles rapidly exit the Galaxy and almost do not produce secondaries. Another explanation, e.g. , that the peaks are produced by CR reacceleration ($`K0`$) on an undisturbed Kolmogorov-type spectrum remains as a viable alternative. The quantitative analysis of this problem is presented by Moskalenko et al. at this Conference.

Many aspects of this topic including the structure of turbulence at the main scale (the theory predicts that $`S_AS_M`$), the anisotropy in $`k`$-space (the waves propagates at large angles to local magnetic field in the Alfvenic turbulence and this should increase the diffusion coefficient compared to (2)), and the consequences for the interpretation of data on interstellar turbulence remain beyond the scope of the present short paper.

Acknowledgments. The work was supported by NASA Astrophysical Theory Program grants. It was also supported by RFBR grant at IZMIRAN.

### 5. References

1. Berezinskii V.S. et al. 1990, Astrophysics of Cosmic Rays (North Holland, Amsterdam)

2. Cho J., Lazarian A. 2002, astro-ph/0205282

3. Goldstein M.L., Roberts D.A., Matthaeus W.H. 1995, ARA&A 33, 283

4. Iroshnikov P. 1963, Sov. Astron. 7, 566

5. Jones F.C. et al. 2001, ApJ 547, 264

6. Kraichnan R. 1965, Phys. Fluids 8, 1385

7. Landau L.D., Lifschitz E.M. 1987, Fluid Mechanics (Pergamon, Oxford)

8. Norman C.A., Ferrara A. 1996, ApJ 467, 280

9. Seo E.S., Ptuskin V.S. 1994, ApJ 431, 705

10. Verma M.K. et al. 1996, JGR 101, 21619


text:
# Gravitational Collapse of Filamentary Magnetized Molecular Clouds

## 1 INTRODUCTION

Filamentary structures are commonly seen in giant molecular clouds as sites in which are embedded star-forming clumps. These structures are threaded by magnetic fields (Heiles et al. 1993; Matthews & Wilson 2000; Falgarone et al. 2001), which can support the filaments against gravitational collapse. While the mechanism by which filaments form is not well understood, it has been suggested that these filamentary structures occur as the result of the fragmentation of the parent molecular cloud through turbulent motions (Klessen & Burkert 2000; Ostriker et al. 2001). Recent theoretical work (Fiege & Pudritz 2000a, b) suggests that helical magnetic fields may play an important role in determining the structure of a filament.

The study of the collapse of a spherical pre-stellar core has been greatly benefitted by the study of one-dimensional models (Larson 1969; Penston 1969; Shu 1977; Foster & Chevalier 1993), which have become very useful in identifying the physical origin of a particular structural phenomenon. By leaving out most of the complicated and analytically intractable physics which may occur in the full three-dimensional system, these models allow us to identify the physical processes which lead to the observed structure, and can be compared to both numerical simulations and observational data of molecular cloud cores. The spherical self-similar models feature an initial gravitational collapse phase that ends with the appearance of a singularity in the density field. This is followed by an accretion phase, wherein an expansion wave propagates outward from the newly formed core enabling the infall of the remaining envelope. Similar analyses of the process of filament formation and evolution have not been as well studied, despite the frequency of filamentary structures forming in numerical simulations (Porter et al. 1994; Klessen & Burkert 2000; Ostriker et al. 2001), and the abundance of filaments in observations of molecular clouds (e.g. Falgarone et al. 2001).

Stodółkiewicz (1963) found solutions for the equilibrium structure of an isothermal cylinder, for a fluid where the ratio of gas pressure to magnetic pressure is constant. Changing the value of this ratio has the effect of changing the scale radius of the solution, which falls off as $`(r/r_0)^4`$. Ostriker (1964) found identical solutions for the unmagnetized isothermal cylinder.

Miyama et al. (1987) derived a set of self-similar solutions for an unmagnetized isothermal cylinder. These solutions have an infall velocity which is proportional to the distance from the axis, and a density structure that has the same form as the unmagnetized equilibrium filaments of Stodółkiewicz (1963) and Ostriker (1964), but with a scale radius that decreases with time and a central density that increases with time.

Recently, Hennebelle (2002) explored a set of self-similar solutions for a rotating, magnetized filament which may undergo collapse in the axial direction, in addition to radial collapse. This set of solutions include the equations which we will derive in the next section, although Hennebelle (2002) does not explore the case where there are only purely radial motions. This is the focus of our analysis.

The goal of this paper is to extend the calculations of the self-similar collapse of an isothermal cylinder to include magnetic and rotational effects. We show that one has to choose a relationship between the magnetic field strength and the density. When this is done, we find that the density and magnetic field structures of collapsing cylinders is identical to that of the corresponding hydrostatic structure, with only the scale radius and normalization of the quantities changing in time. The important consequence of this for observations is that only careful large-scale velocity measurements can discriminate between static and evolving filamentary structure, unlike the situation for accreting spheres. We derive and solve our set of self-similar equations in Section 2. These solutions will be discussed in Section 3.

## 2 THE SELF-SIMILAR COLLAPSE EQUATIONS

We begin our analysis with the equations of ideal magnetohydrodynamics for a self-gravitating gas:

$`{\displaystyle \frac{\rho }{t}}+\left(\rho 𝐯\right)`$ $`=`$ $`0`$ (1)
$`{\displaystyle \frac{𝐁}{t}}`$ $`=`$ $`\times \left(𝐯\times 𝐁\right)`$ (2)
$`\rho {\displaystyle \frac{𝐯}{t}}+\rho \left(𝐯\right)𝐯`$ $`=`$ $`P+{\displaystyle \frac{1}{\mu }}\left[\left(𝐁\right)𝐁{\displaystyle \frac{1}{2}}\left(𝐁^2\right)\right]\rho \mathrm{\Phi }`$ (3)
$`^2\mathrm{\Phi }`$ $`=`$ $`4\pi G\rho `$ (4)
In these equations, $`\rho `$ is the density, $`𝐯`$ is the velocity, $`𝐁`$ is the magnetic field, $`P`$ is the thermal pressure, for which we adopt the isothermal equation of state $`P=c_s^2\rho `$ (with a constant speed of sound $`c_s`$), and $`\mathrm{\Phi }`$ is the gravitational potential.

The geometry of our filaments leads us to use a cylindrical coordinate system to solve the equations. Our cylindrical symmetry causes all $`/\varphi `$ and $`/z`$ derivatives to be zero, so that the resulting fields are functions of $`(r,t)`$ only. This geometry allows us to replace the gravitational potential by a mass per unit length $`m`$, defined by

$`{\displaystyle \frac{m}{r}}`$ $`=`$ $`2\pi r\rho `$ (5)
We can make the following substitutions, to change our equations from functions of ($`r,t`$) to separable functions of ($`s,t`$), for the self-similar variable $`s=r/(c_st)`$. These solutions will represent collapsing motions for $`s<0`$ (and hence $`t<0`$), and expanding motions for $`s>0`$.

$`\rho (r,t)`$ $`=`$ $`{\displaystyle \frac{1}{\pi Gt^2}}\stackrel{~}{\rho }(s)`$ (6)
$`m(r,t)`$ $`=`$ $`{\displaystyle \frac{2c_s^2}{G}}\stackrel{~}{m}(s)`$ (7)
$`𝐯(r,t)`$ $`=`$ $`c_s\stackrel{~}{𝐯}(s)`$ (8)
$`𝐁(r,t)`$ $`=`$ $`\sqrt{{\displaystyle \frac{\mu }{\pi G}}}{\displaystyle \frac{c_s}{t}}\stackrel{~}{𝐁}(s)`$ (9)
where $`\stackrel{~}{\rho }(s)`$ etc. are dimensionless functions of the self-similar variable $`s`$ alone. Substituting these into Equations (1)-(5) (see Appendix A for details),

$`\stackrel{~}{v}_r`$ $`=`$ $`s`$ (10)
$`{\displaystyle \frac{d\stackrel{~}{m}}{ds}}`$ $`=`$ $`s\stackrel{~}{\rho }`$ (11)
$`{\displaystyle \frac{d}{ds}}\left(s\stackrel{~}{B}_r\right)`$ $`=`$ $`0`$ (12)
$`{\displaystyle \frac{d}{ds}}\left(\stackrel{~}{v}_\varphi \stackrel{~}{B}_r\right)`$ $`=`$ $`0`$ (13)
$`\stackrel{~}{B}_r{\displaystyle \frac{d\stackrel{~}{B}_z}{ds}}`$ $`=`$ $`0`$ (14)
$`\stackrel{~}{B}_z`$ $`=`$ $`{\displaystyle \frac{1}{s}}{\displaystyle \frac{d}{ds}}\left(s\stackrel{~}{v}_z\stackrel{~}{B}_r\right)`$ (15)
$`{\displaystyle \frac{\stackrel{~}{\rho }\stackrel{~}{v}_\varphi ^2}{s}}`$ $`=`$ $`{\displaystyle \frac{d\stackrel{~}{\rho }}{ds}}+{\displaystyle \frac{4\stackrel{~}{m}\stackrel{~}{\rho }}{s}}+{\displaystyle \frac{\stackrel{~}{B}_\varphi }{s}}{\displaystyle \frac{d}{ds}}\left(s\stackrel{~}{B}_\varphi \right)+\stackrel{~}{B}_z{\displaystyle \frac{d\stackrel{~}{B}_z}{ds}}`$ (16)
$`\stackrel{~}{\rho }\stackrel{~}{v}_\varphi `$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{B}_r}{s}}{\displaystyle \frac{d}{ds}}\left(s\stackrel{~}{B}_\varphi \right)`$ (17)
Before we proceed to solve Equations (10)–(17), we will first recap the solutions to these equations in the case of $`𝐁=0`$. Equation (11) remains unchanged, Equation (17) shows that $`\stackrel{~}{v}_\varphi =0`$, and Equation (16) reduces to $`\mathrm{d}(\mathrm{ln}\stackrel{~}{\rho })/\mathrm{d}(\mathrm{ln}s)+4\stackrel{~}{m}=0`$. The solution to this, combined with Equation (11), is $`\stackrel{~}{\rho }=2\zeta \left(1+\zeta s^2\right)^2`$ and $`\stackrel{~}{m}=\zeta s^2\left(1+\zeta s^2\right)^2`$, where $`\zeta `$ is a constant of integration. This solution was first found by Miyama et al. (1987). It has two features worth highlighting:

1. It has the same radial density structure as the equilibrium solution of Stodółkiewicz (1963) and Ostriker (1964).
2. The total enclosed mass per unit length (i.e. as $`|s|\mathrm{}`$) is the critical mass per unit length.

The asymptotic density behaviour of this solution tends to $`s^4`$ at large values of s; we will show later that magnetic solutions tend to have this behaviour only in the limit of low magnetization.

There are two classes of solutions to Equations (10)–(17) when the effects of magnetic fields are included, which we can separate through the use of Equation (14):

1. $`\stackrel{~}{B}_r=0`$: In this case, we can show from Equation (15) that $`\stackrel{~}{B}_z=0`$ and from Equation (17) that $`\stackrel{~}{v}_\varphi =0`$ as well. This solution also has the consequence that our system of equations is underspecified, as $`\stackrel{~}{B}_\varphi `$ only appears in the radial momentum equation. We will discuss the solutions to this below.
2. $`\mathrm{d}\stackrel{~}{B}_z/\mathrm{d}s=0`$: If our $`\stackrel{~}{B}_z`$ field is constant, then from Equation (12) we can see that we have a radial field $`\stackrel{~}{B}_r=\stackrel{~}{B}_{r0}/s`$ which is singular along the cylinder axis.

The latter solution is difficult to justify physically; even though the field is still divergence-free, the field lines are not closed because of the linear singularity at the axis of symmetry. Henceforth we discard this possibility and focus our attention on Case 1.

We now proceed to solve the remaining set of equations. With the condition that $`\stackrel{~}{B}_r=0`$, we are left with

$`\stackrel{~}{v}_r`$ $`=`$ $`s`$ (18)
$`{\displaystyle \frac{d\stackrel{~}{m}}{ds}}`$ $`=`$ $`s\stackrel{~}{\rho }`$ (19)
$`0`$ $`=`$ $`{\displaystyle \frac{d\stackrel{~}{\rho }}{ds}}+{\displaystyle \frac{4\stackrel{~}{m}\stackrel{~}{\rho }}{s}}+{\displaystyle \frac{\stackrel{~}{B}_\varphi }{s}}{\displaystyle \frac{d}{ds}}\left(s\stackrel{~}{B}_\varphi \right)`$ (20)
We can integrate Equation (20) to obtain the following relationship between $`\stackrel{~}{\rho }`$, $`\stackrel{~}{m}`$ and $`\stackrel{~}{B}_\varphi `$:

$`\alpha `$ $`=`$ $`s^2\left(2\stackrel{~}{\rho }+\stackrel{~}{B}_\varphi ^2\right)+4\stackrel{~}{m}\left(\stackrel{~}{m}1\right)`$ (21)
We can constrain the value of the constant $`\alpha `$ by examining the behaviour of the system on the axis. Since $`\underset{s0}{lim}\stackrel{~}{m}=0`$, it follows that in the limit $`s0`$, $`\underset{s0}{lim}s^2(2\stackrel{~}{\rho }+\stackrel{~}{B}_\varphi ^2)=\alpha `$.

The first term in this result can be evaluated by using Equation (19). We can express the density as $`\stackrel{~}{\rho }=A|s|^n`$ as $`s0`$, as for sufficiently small $`s`$, one term in the expansion of $`\stackrel{~}{\rho }`$ will dominate the other terms. From Equation (19), we can write the mass per unit length for this density behaviour as $`\stackrel{~}{m}=A_0^{|s|}|s|^{n+1}\mathrm{d}|s|`$. For $`n2`$, $`\stackrel{~}{m}`$ will not be finite in the limit $`s0`$, so with the constraint $`n>2`$, the first term in Equation (21) has the limit $`\underset{s0}{lim}\stackrel{~}{\rho }s^2=0`$.

We can constrain the magnetic field behaviour in a similar manner, by examining the toroidal magnetic flux, $`\mathrm{\Psi }`$. We use the formulation of the magnetic flux of Fiege & Pudritz (2000a), such that $`\mathrm{\Psi }=\stackrel{~}{B}_\varphi \mathrm{da}`$, where $`\mathrm{da}=\mathrm{dr}\mathrm{dz}`$ is the area element through which the magnetic field passes (see Figure 1). If we consider the flux per unit length $`\psi =\mathrm{d}\mathrm{\Psi }/\mathrm{dz}`$, and write the magnetic field as $`\stackrel{~}{B}_\varphi =C|s|^n^{}`$ as $`s0`$, we find that $`\psi =C_0^s|s|^n^{}ds`$. For $`n^{}1`$, $`\psi `$ is not finite, so we restrict $`n^{}>1`$. We thus establish that $`\alpha =0`$ in Equation (21).

The properties of the system can be examined if we solve for the mass per unit length in Equation (21).

$`\stackrel{~}{m}_\pm `$ $`=`$ $`{\displaystyle \frac{1}{2}}\pm {\displaystyle \frac{1}{2}}\sqrt{1s^2\left(2\stackrel{~}{\rho }+\stackrel{~}{B}_\varphi ^2\right)}`$ (22)
We can identify two branches to this solution. The minus branch $`\stackrel{~}{m}_{}1/2`$ corresponds to the minus sign in Equation (22); the plus solution $`\stackrel{~}{m}_+1/2`$ is the other branch. The two branches meet at the half-mass point given by $`s_{1/2}^2\left[2\stackrel{~}{\rho }(s_{1/2})+\stackrel{~}{B}_\varphi ^2(s_{1/2})\right]=1`$. Since $`\stackrel{~}{m}(0)=0`$, we know that the critical quantity $`q_{\mathrm{crit}}=s^2(2\stackrel{~}{\rho }+\stackrel{~}{B}_\varphi ^2)0`$ as $`s0`$. Thus, we can determine that for every solution, the quantity $`q_{\mathrm{crit}}`$ must increase from $`01`$, then begin to decrease again. Thus, the minus branch begins at $`\stackrel{~}{m}_{}(0)=0`$, then increases to $`1/2`$, then decreases again; the plus branch starts at $`\stackrel{~}{m}_+(0)=1`$, decreases to $`1/2`$, then increase again.

Since the plus branch is always greater than $`1/2`$, it cannot represent the solution at $`s=0`$; similarly, at values of $`|s|>|s_{1/2}|`$, the mass per unit length must continue to increase (as a decreasing mass per unit length implies a negative density according to Equation (19), which is unphysical), and so the minus solution is unphysical in this regime. Thus the physical solution must begin on the minus branch at low $`s`$, then at the half-mass point switch to the plus branch. The maximum value of $`\stackrel{~}{m}`$ must occur when $`\stackrel{~}{\rho }=0`$. If this occurs at a finite radius $`s_{\mathrm{max}}`$, then the total mass per unit length will be

$`\stackrel{~}{m}_{\mathrm{max}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\pm {\displaystyle \frac{1}{2}}\sqrt{1s_{\mathrm{max}}^2\stackrel{~}{B}_\varphi ^2(s_{\mathrm{max}})}`$ (23)
where the plus or minus sign is determined by whether $`|s_{\mathrm{max}}|>|s_{1/2}|`$ or $`|s_{\mathrm{max}}|<|s_{1/2}|`$, respectively. This $`\stackrel{~}{m}_{\mathrm{max}}`$ is identical to the reduced critical mass per unit length $`m_{\mathrm{mag}}`$ found by Fiege & Pudritz (2000a) from the virial theorem, for a filament with a toroidal magnetic field.

Because of the indeterminacy in this solution as noted above, we need to specify a form for $`\stackrel{~}{B}_\varphi `$ in order to make further progress. We will perform this calculation for two plausible magnetic field configurations which commonly appear in star formation literature: one where the ratio of the thermal to magnetic pressure $`\beta =2\mu c_s^2\rho /B_\varphi ^2=2\stackrel{~}{\rho }/\stackrel{~}{B}_\varphi ^2`$ is constant, and one where the toroidal flux-to-mass ratio $`\mathrm{\Gamma }_\varphi =\stackrel{~}{B}_\varphi /(|s|\stackrel{~}{\rho })`$ is constant (Fiege & Pudritz 2000a). We will discuss these two cases in turn.

## 3 RESULTS

There are two models for the magnetization of the gas in molecular clouds which are used extensively in the literature: one in which the magnetic pressure is proportional to the gas pressure, and one in which the magnetic flux through a fluid element is proportional to the mass of that element. As we have the freedom in our set of equations to specify one of our variables, we choose to constrain the magnetic field through one of these relations.

As we show in the Appendix, the structure of the self-similar collapse equations (19)-(20) is identical in form to the equations for the hydrostatic solutions; the density and magnetic fields have identical radial structures. The only difference between the static and collapse solutions is that the collapse solutions are time-dependent, and have an infall velocity. It is remarkable that the radial velocity terms in the momentum, induction and continuity equations exactly cancels out the time-dependent terms. A very different behaviour is seen in spherical symmetry, for which the collapse solutions have a different density structure than the static solutions. This underlies the importance that the geometry has in determining the spatial structure of the collapsing filament.

In self-similar cylindrical coordinates, the continuity equation (1) restricts the radial velocity to a single solution, $`v_r=r/t`$. Thus, this analysis can only represent a collapsing solution for $`t<0`$. Furthermore, this implies that there is no equivalent to the accretion solution found in the self-similar collapse of a spherical isothermal fluid by Larson (1969), Penston (1969), and Shu (1977). The collapse of filaments, therefore, is very different than the collapse of spheres.

### 3.1 Case 1A: Constant $`\beta `$

A common description for magnetic fields in molecular clouds is a constant ratio of gas pressure to magnetic pressure, $`\beta `$, such that the magnetic field is $`\stackrel{~}{B}_\varphi =\left(2\beta ^1\stackrel{~}{\rho }\right)^{1/2}`$. We can use this in Equation (22) to write

$`\stackrel{~}{m}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\pm {\displaystyle \frac{1}{2}}\sqrt{12s^2(1+\beta ^1)\stackrel{~}{\rho }}`$ (24)
We can differentiate this and use Equation (19) to solve for $`\stackrel{~}{\rho }(s)`$, with the result

$`\stackrel{~}{\rho }(s)`$ $`=`$ $`{\displaystyle \frac{\mathrm{sech}^2\left[(\mathrm{ln}|s|)/(1+\beta ^1)\right]}{2s^2\left(1+\beta ^1\right)}}`$ (25)
$`\stackrel{~}{m}(s)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left\{1+\mathrm{tanh}\left[(\mathrm{ln}|s|)/(1+\beta ^1)\right]\right\}`$ (26)
$`\stackrel{~}{B}_\varphi (s)`$ $`=`$ $`{\displaystyle \frac{\mathrm{sech}\left[(\mathrm{ln}|s|)/(1+\beta ^1)\right]}{|s|\sqrt{1+\beta }}}`$ (27)
This solution displays a number of the properties mentioned previously. The half-mass point occurs when $`|s|=1`$; the value of $`s_{1/2}^2(2\stackrel{~}{\rho }(s_{1/2})+\stackrel{~}{B}_\varphi ^2(s_{1/2}))=\mathrm{sech}^2(0)=1`$, as expected. This is also the sonic point in the fluid (that is, where $`|\stackrel{~}{v}_r|=1`$). The limit $`\underset{|s|\mathrm{}}{lim}(s\stackrel{~}{B}_\varphi )=0`$, so the total mass per unit length ought to be equal to the critical value $`\stackrel{~}{m}_{\mathrm{max}}=1`$; from Equation (26) we can see that this is indeed the case. The solutions to Equations (25)–(27) are plotted in Figure 2. The behaviour of the quantity $`q_{\mathrm{crit}}=s^2\left(2\stackrel{~}{\rho }+\stackrel{~}{B}_\varphi ^2\right)`$ is plotted in Figure 3a, and displays exactly the behaviour discussed in Section (2). The branches of $`\stackrel{~}{m}`$ are shown in Figure 3b.

Since most measurements of molecular clouds have difficulty resolving the interior regions of cores and filaments, we examine the behaviour of these solutions at large $`s`$ (equivalent to large $`r`$). We can do this by expanding the $`\mathrm{sech}^2`$ term into exponentials. The net result is

$`\stackrel{~}{\rho }|s|^{2\frac{2\beta +1}{\beta +1}}`$ (28)
for $`|s|1`$. Equation (28) shows that for $`\beta 1`$ (corresponding to weak magnetic fields), $`\stackrel{~}{\rho }s^4`$. For $`\beta 1`$ (strong magnetic fields), $`\stackrel{~}{\rho }s^2`$. Note that in the limit $`s\mathrm{}`$ (equivalent to $`t0`$), $`\stackrel{~}{\rho }0`$ for all $`r0`$. This simply shows that there cannot be an accretion solution for $`t>0`$—there is no material left to be accreted.

We can similarly write an expression for the density behaviour as $`s0`$:

$`\stackrel{~}{\rho }|s|^{\frac{2}{1+\beta }}`$ (29)
For a finite value of $`\beta `$, the density will be singular at the axis. This is in contrast to the result of Stodółkiewicz (1963), who used an arbitrary boundary condition of $`\rho (r=0)=\mathrm{𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡}`$. We believe our result to be physically more self-consistent, as it satisfies the physical constraints of zero mass per unit length and toroidal flux at the origin (the solution of Stodółkiewicz (1963) cannot simultaneously satisfy the hydrodynamic equations and the restriction of $`m(0)=0`$, for a purely toroidal field). Our solution shows that the density in a $`\beta =\mathrm{constant}`$ model will be singular on the axis. A singularity in the spherical solutions can be rationalized as an approximation to the effects of a protostar in the center of the cloud. In the cylindrical case, even though collapse as a singular filament occurs, this structure itself will undergo fragmentation in a finite amount of time, as is shown in Section 3.3. Our treatment of the $`\beta =\mathrm{constant}`$ case shows that it may be of limited astrophysical interest.

The toroidal field contributes both a magnetic pressure force and a magnetic tension force to the momentum equations. It is instructive to compare the magnitudes of these two forces. The magnetic pressure is $`(1/2)\mathrm{d}(\stackrel{~}{B}_\varphi ^2)/\mathrm{d}s`$, and the magnetic tension force is $`\stackrel{~}{B}_\varphi ^2/s`$. The ratio is thus

$`R`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{\stackrel{~}{𝐁}^2}{(\stackrel{~}{𝐁})\stackrel{~}{𝐁}}}={\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{ln}(\stackrel{~}{B}_\varphi ^2)}{\mathrm{ln}|s|}}`$ (30)
For our $`\beta =\mathrm{𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡}`$ fluid, $`R=(1/2)\mathrm{d}(\mathrm{ln}\stackrel{~}{\rho })/\mathrm{d}(\mathrm{ln}|s|)=2\left(2\stackrel{~}{m}+\beta ^1\right)`$. The magnetic pressure force will dominate for $`|R|>1`$, which corresponds to $`\stackrel{~}{m}>1/4(1/2)\beta ^1`$; note that this will always be the case for a fluid with $`\beta <2`$, as then the right-hand side will be negative. Strongly magnetized filaments are therefore supported by the pressure forces associated with their toroidal fields.

### 3.2 Case 1B: Constant Flux-to-Mass Ratio

Another form for the magnetic field configuration of a filament that has been proposed (Fiege & Pudritz 2000a) is one wherein the magnetic flux through a volume is proportional to the mass in the volume, such that the flux-to-mass ratio $`\mathrm{\Gamma }_\varphi `$ is constant everywhere ($`\stackrel{~}{B}_\varphi =\mathrm{\Gamma }_\varphi |s|\rho `$). With this assumption, we can numerically integrate Equation (20) using a Runge-Kutta integrator to obtain a set of solutions. We plot these solutions for several values of $`\mathrm{\Gamma }_\varphi `$ in Figure 4. As for the solutions for the $`\beta =\mathrm{𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡}`$ plasma, the solutions for $`\mathrm{\Gamma }_\varphi =\mathrm{𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡}`$ exhibit the properties discussed in Section (2). However, we also note some differences in the behaviour of these two different systems. We plot the critical quantity $`q_{\mathrm{crit}}`$ for $`\mathrm{\Gamma }_\varphi =\mathrm{𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡}`$ in Figure 3c, which can be directly compared to the $`q_{\mathrm{crit}}`$ of the $`\beta =\mathrm{𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡}`$ gas, plotted in the panel above it. While the half-mass point for the $`\beta =\mathrm{𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡}`$ case did not depend on the value of the parameter $`\beta `$, for the $`\mathrm{\Gamma }_\varphi =\mathrm{𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡}`$ solution the half-mass point does depend on the value of $`\mathrm{\Gamma }_\varphi `$. The density and magnetic field for the $`\mathrm{\Gamma }_\varphi =\mathrm{𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡}`$ solution are also non-singular, as can be seen in Figure 4.

At low values of magnetization (i.e. low $`\mathrm{\Gamma }_\varphi `$), the density distribution is very similar to that predicted by Miyama et al. (1987) for the self-similar collapse of an unmagnetized filament (which in turn has a the same shape as the unmagnetized hydrostatic filament). As the magnetic flux increases, the density distribution at large $`|s|`$ tends towards $`\stackrel{~}{\rho }s^2`$ rather than the steeper $`s^4`$ for weak magnetic fields. This can be shown explicitly by isolating $`d\stackrel{~}{\rho }/ds`$ in Equation (20):

$`{\displaystyle \frac{d\stackrel{~}{\rho }}{ds}}`$ $`=`$ $`\left({\displaystyle \frac{2\stackrel{~}{\rho }}{s}}\right){\displaystyle \frac{2\stackrel{~}{m}+\mathrm{\Gamma }_\varphi ^2s^2\stackrel{~}{\rho }}{1+\mathrm{\Gamma }_\varphi ^2s^2\stackrel{~}{\rho }}}`$ (31)
For $`\mathrm{\Gamma }_\varphi 1`$, we see that $`d\mathrm{ln}(\stackrel{~}{\rho })/d\mathrm{ln}|s|=4\stackrel{~}{m}`$; since $`m1`$ as $`|s|\mathrm{}`$, $`\stackrel{~}{\rho }s^4`$ in the case of low magnetization. For $`\mathrm{\Gamma }_\varphi 1`$, we see that $`d\mathrm{ln}(\stackrel{~}{\rho })/d\mathrm{ln}|s|=2`$, implying that $`\stackrel{~}{\rho }s^2`$, as discussed above. These results and those in Section 3.1 show that the asymptotic behaviour at large radii for very weak or very strong magnetic fields is independent of the details of how these filaments are magnetized.

As for the constant-$`\beta `$ solution, as $`s\mathrm{}`$, $`\stackrel{~}{\rho }0`$. Thus, for solutions with a constant $`\mathrm{\Gamma }_\varphi `$, there cannot be an accretion solution for $`t>0`$, either.

### 3.3 Fragmentation

Stodółkiewicz (1963) and Fiege & Pudritz (2000b) derived criteria for the fragmentation of a magnetized isothermal cylinder, for each of the forms of the magnetization that we have used. Since our collapse solution moves towards a line singularity at t=0, which is unphysical, there exists the possibility that the cylinder will become unstable to fragmentation instabilities as the density and magnetic field strengths change with time.

Fiege & Pudritz (2000b) found that for a filament with constant toroidal flux-to-mass ratio, the fragmentation timescale is

$`\tau _{frag}=1.8\left({\displaystyle \frac{\omega _{max}^2}{0.01}}\right)^{1/2}\left({\displaystyle \frac{n_c}{10^4\mathrm{cm}^3}}\right)\mathrm{Myr}`$ (32)
where $`w_{max}^2`$ is a frequency of fragmentation which depends on the value of $`\mathrm{\Gamma }_\varphi `$. Fiege & Pudritz (2000b) found that $`\omega _{max}^20.0581`$.

As the radial density structure is the same in both the self-similar case and the hydrostatic case, we can compare the fragmentation timescale to our collapse timescale, to see which process should dominate. The collapse timescale is

$`\tau _{col}={\displaystyle \frac{1}{\sqrt{\pi G\rho _c}}}0.38\left({\displaystyle \frac{n_c}{10^4\mathrm{cm}^3}}\right)^{1/2}\mathrm{Myr}`$ (33)
and so the ratio of the fragmentation to collapse timescales is

$`{\displaystyle \frac{\tau _{frag}}{\tau _{col}}}=4.75\left({\displaystyle \frac{\omega _{max}^2}{0.01}}\right)^{1/2}`$ (34)
This ratio has a minimum at $`\mathrm{\Gamma }_\varphi =0`$ of $`\tau _{frag}/\tau _{col}=1.97`$, which increases as $`\mathrm{\Gamma }_\varphi `$ increases. Thus, the toroidal magnetic field acts to prevent the cloud from fragmenting. In the case of low magnetization, the collapse and fragmentation timescales are approximately equal, suggesting that the cloud may begin to collapse along the filament axis, in addition to its radial collapse.

We can also find the ratio of the critical wavelength for fragmentation to the radial collapse scale. Fiege & Pudritz (2000b) found that the fragmentation wavelength for a constant toroidal flux-to-mass ratio is

$`\lambda _{frag}=2.8\left({\displaystyle \frac{\sigma _c}{0.5\mathrm{kms}^{}1}}\right)\left({\displaystyle \frac{n_c}{10^4\mathrm{cm}^3}}\right)^{1/2}\left({\displaystyle \frac{k_{max}}{0.2}}\right)^1\mathrm{pc}`$ (35)
where $`k_{max}`$ is a dimensionless wavenumber corresponding to the maximum instability, and has a maximum value of 0.462 at $`\mathrm{\Gamma }_\varphi =0`$. An estimate of the radial length scale can be set by setting $`|s|=1`$:

$`\lambda _{rad}`$ $`=`$ $`c_sst`$
$`\lambda _{rad}`$ $`=`$ $`c_s{\displaystyle \frac{1}{\sqrt{\pi G\rho _c}}}`$
$`\lambda _{rad}`$ $`=`$ $`0.194\left({\displaystyle \frac{\sigma _c}{0.5\mathrm{kms}^{}1}}\right)\left({\displaystyle \frac{n_c}{10^4\mathrm{cm}^3}}\right)^{1/2}\mathrm{pc}`$ (36)
From this, we get a ratio of fragmentation to radial wavelengths of

$`{\displaystyle \frac{\lambda _{frag}}{\lambda _{rad}}}=14.43\left({\displaystyle \frac{k_{max}}{0.2}}\right)^1`$ (37)
At $`\mathrm{\Gamma }_\varphi =0`$, this has its minimum value of 6.25 (3.41 if we use the critical wavelength rather than the maximum wavelength). Again, this suggests that toroidal fields support the cylinder against collapse.

We can do a similar analysis for the $`\beta =\mathrm{𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡}`$ case. The fragmentation wavelength found by Stodółkiewicz (1963) for a $`\beta =\mathrm{𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡}`$ fluid is

$`\lambda _{frag}=1.08\left({\displaystyle \frac{\sigma _c}{0.5\mathrm{kms}^{}1}}\right)\left({\displaystyle \frac{n_c}{10^4\mathrm{cm}^3}}\right)^{1/2}\sqrt{1+{\displaystyle \frac{1}{\beta }}}\mathrm{pc}`$ (38)
while the radial wavelength is still given by Equation (36), for a ratio of

$`{\displaystyle \frac{\lambda _{frag}}{\lambda _{rad}}}=5.563\sqrt{1+{\displaystyle \frac{1}{\beta }}}`$ (39)
Again, this has a minimum value when the magnetization is zero ($`\beta =\mathrm{}`$); the presence of a toroidal field serves to stabilize the cylinder against fragmentation.

In addition to the axisymmetric instabilities discussed above, there exists the possibility that non-axisymmetric instabilities (e.g. kink modes) could develop. The presence of such instabilities would destroy the axisymmetric nature of our models, a problem that cannot be addressed with a self-similar formalism. We can, however, estimate the likelihood of such modes by comparing their growth rates with that of the axisymmetric modes. Nakamura et al. (1993) has investigated the kink mode instabilities of a hydrostatic molecular cloud with helical magnetic fields described by a constant $`\beta `$. They found that the axisymmetric instabilities, such as those described above, have faster growth rates than all non-axisymmetric modes. Thus, axisymmetric modes will dominate the fragmentation, and non-axisymmetric modes are not an important effect in the evolution of these filaments on the collapse time scales that we are considering.

## 4 DISCUSSION AND CONCLUSIONS

We have determined the structure of a collapsing self-similar magnetized isothermal cylinder for the two most common magnetic field configurations found in the literature. One of the major conclusions of our study is that the cylindrical radial structure for the gas density is indistinguishable from the hydrostatic structure with the same magnetic field configuration. Consequently, one would need to examine the profiles of molecular lines for infall signatures in order to differentiate collapsing and equilibrium structures. Because the high-density regions that emit most of the thermal radiation will have low velocities, it may be more difficult to detect a characteristic line profile in comparison with the situation for accreting spherical cores that can be identified by an inverse P-Cygni profile.

The self-similar density profile at large radii from the axis follows a power law. For a cylinder with a weak magnetic field, we find that at any particular time during the collapse, $`\rho r^4`$. By way of contrast, for a cylinder with a strong magnetic field, we find that $`\rho r^2`$. These results are independent of the model for the magnetic field strength that we use ($`\beta =\mathrm{𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡}`$ or $`\mathrm{\Gamma }_\varphi =\mathrm{𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡}`$).

Our solutions require that there be no motions along the axis, and that the magnetic field is purely toroidal. In a fully dynamic, three-dimensional collapse, we would not necessarily expect this to be the case. However, if the axial motions and poloidal fields are sufficiently weak, we could reasonably expect our solutions to be a good approximation. Furthermore, the static models of Fiege & Pudritz (2000a) indicate that for models consisting of both poloidal and toroidal fields, the toroidal magnetic field strongly dominates in the outer regions of a filament, suggesting that our solutions may be applicable to the envelopes of filaments threaded with helical fields.

The velocity behaviour of our solutions bears a close resemblence to the solution found by Penston (1969) for the self-similar collapse of an isothermal sphere, in that the infall velocity $`v_r0`$ as $`r0`$. However, due to the simple nature of the continuity equation, we only find the collapse solution; there is no analogue to the Shu (1977) accretion–expansion wave solution that follows the formation of a singular core. It is also the continuity equation which requires $`vr`$ at all radii, rather than having an asymptotic infall velocity as found by Penston (1969). Clearly, arbitrarily high velocities at large radii are not physical; there needs to be a mechanism to truncate the models at finite $`r`$. Pressure truncation, which we have not included in this analysis, could satisfy this condition.

We find, in agreement with Stodółkiewicz (1963) and Fiege & Pudritz (2000b), that the presence of a toroidal field stabilizes the cylinder against fragmentation. Even in the absence of a field, however, collapse motions should dominate over fragmentation. We thus have shown that self-similar filaments should remain relatively unfragmented as they collapse. This mechanism would allow filamentary structures to survive over long periods of time while retaining their filamentary character. If this were not true, we would not expect to see many filaments because they would fragment into cores too quickly.

Although large scale velocity fields have been observed to be associated with individual cloud cores (e.g. Mardones et al. 1997; Lee et al. 1999), we emphasize that the flows predicted here are on a yet larger scale associated with filament formation. Simulations of hydromagnetic turbulence show that filaments can form from shocks (Ostriker et al. 1999, 2001). We plan to investigate the possible link between our analytic models and time-dependent simulations of magnetized clouds.

We thank an anonymous referee for useful comments on the manuscript. D.A.T. is supported by an Ontario Graduate Scholarship. R.E.P. is supported by the Natural Sciences and Engineering Research Council of Canada.

## Appendix A THE STRUCTURE OF THE SELF-SIMILAR AND HYDROSTATIC EQUATIONS

We can show explicitly that the structure of a filament that is collapsing self-similarly is identical to that of a hydrostatic filament with the same constitutive equations for the magnetic field. In our formalism, we replace the gravitational potential with the mass per unit length (Equation 5)

$`r{\displaystyle \frac{\mathrm{\Phi }}{r}}=2Gm`$ (A1)
We write out Equations (1)-(4) in component form, and use the fact that we are considering only the axisymmetric and axially symmetric case (i.e. functions depend only on the radius r, and the time t).

$`{\displaystyle \frac{m}{t}}+v_r{\displaystyle \frac{m}{r}}`$ $`=`$ $`0`$ (A2)
$`{\displaystyle \frac{B_r}{t}}`$ $`=`$ $`0`$ (A3)
$`{\displaystyle \frac{B_\varphi }{t}}`$ $`=`$ $`{\displaystyle \frac{}{r}}\left(v_rB_\varphi v_\varphi B_r\right)`$ (A4)
$`{\displaystyle \frac{B_z}{t}}`$ $`=`$ $`{\displaystyle \frac{1}{r}}{\displaystyle \frac{}{r}}\left[r\left(v_zB_rv_rB_z\right)\right]`$ (A5)
$`\rho {\displaystyle \frac{v_r}{t}}+\rho v_r{\displaystyle \frac{v_r}{r}}\rho {\displaystyle \frac{v_\varphi ^2}{r}}`$ $`=`$ $`{\displaystyle \frac{1}{\mu }}\left[B_r{\displaystyle \frac{B_r}{r}}{\displaystyle \frac{B_\varphi ^2}{r}}{\displaystyle \frac{1}{2}}{\displaystyle \frac{}{r}}\left(B_r^2+B_\varphi ^2+B_z^2\right)\right]`$ (A6)
$`c_s^2{\displaystyle \frac{\rho }{r}}{\displaystyle \frac{2Gm\rho }{r}}`$
$`\rho {\displaystyle \frac{v_\varphi }{t}}+\rho \left(v_r{\displaystyle \frac{v_\varphi }{r}}+{\displaystyle \frac{v_\varphi v_r}{r}}\right)`$ $`=`$ $`{\displaystyle \frac{1}{\mu }}\left[B_r{\displaystyle \frac{B_\varphi }{r}}+{\displaystyle \frac{B_\varphi B_r}{r}}\right]`$ (A7)
$`\rho {\displaystyle \frac{v_z}{t}}+\rho v_r{\displaystyle \frac{v_z}{r}}`$ $`=`$ $`{\displaystyle \frac{B_r}{\mu }}{\displaystyle \frac{B_z}{r}}`$ (A8)
If we make the substitutions given in Equations (6)-(9), we obtain the self-similar structure equations (10)-(17). If we instead consider the equilibrium solutions (so that $`𝐯=\mathrm{𝟎}`$, $`/t=0`$), we get

$`0`$ $`=`$ $`{\displaystyle \frac{1}{\mu }}\left[{\displaystyle \frac{B_\varphi }{r}}{\displaystyle \frac{d}{dr}}\left(rB_\varphi \right)+B_z{\displaystyle \frac{dB_z}{dr}}\right]+c_s^2{\displaystyle \frac{d\rho }{dr}}+{\displaystyle \frac{2Gm\rho }{r}}`$ (A9)
$`0`$ $`=`$ $`B_r{\displaystyle \frac{d}{dr}}\left(rB_\varphi \right)`$ (A10)
$`0`$ $`=`$ $`B_r{\displaystyle \frac{dB_z}{dr}}`$ (A11)
$`{\displaystyle \frac{dm}{dr}}`$ $`=`$ $`2\pi r\rho `$ (A12)
From the condition that $`𝐁=0`$, we require $`B_r=0`$, for the same reasons we had in Section 2 for Equation (12). We now convert to dimensionless quantities, using (in analogy with Equations 6-9)

$`\rho `$ $`=`$ $`{\displaystyle \frac{1}{\pi G}}\overline{\rho }`$ (A13)
$`m`$ $`=`$ $`{\displaystyle \frac{2c_s^2}{G}}\overline{m}`$ (A14)
$`B_\varphi `$ $`=`$ $`c_s\sqrt{{\displaystyle \frac{\mu }{\pi G}}}\overline{B}_\varphi `$ (A15)
$`B_z`$ $`=`$ $`c_s\sqrt{{\displaystyle \frac{\mu }{\pi G}}}\overline{B}_z`$ (A16)
$`r`$ $`=c_sx`$ (A17)
This results in

$`{\displaystyle \frac{d\overline{m}}{dx}}`$ $`=`$ $`x\overline{\rho }`$ (A18)
$`0`$ $`=`$ $`{\displaystyle \frac{d\overline{\rho }}{dx}}+{\displaystyle \frac{4\overline{m}\overline{\rho }}{x}}+{\displaystyle \frac{\overline{B}_\varphi }{x}}{\displaystyle \frac{d}{dx}}\left(x\overline{B}_\varphi \right)+{\displaystyle \frac{d\overline{B}_z}{dx}}`$ (A19)
Comparing Equations (A18) and (A19) with Equations (19) and (20), we see that in the cylindrically symmetric geometry, the self-similar collapse equations and hydrostatic equilibrium equations will always give the same structure for the density (with the exception that a $`B_z`$ field in the hydrostatic case is allowed), given a constitutive relation between the toroidal magnetic field and the density.