text:
# DNA hybridization to mismatched templates: a chip study

## Abstract

High-density oligonucleotide arrays are among the most rapidly expanding technologies in biology today. In the GeneChip system, the reconstruction of the target concentration depends upon the differential signal generated from hybridizing the target RNA to two nearly identical templates: a perfect match (PM) and a single mismatch (MM) probe. It has been observed that a large fraction of MM probes repeatably bind targets better than the PMs, against the usual expectation from sequence-specific hybridization; this is difficult to interpret in terms of the underlying physics. We examine this problem via a statistical analysis of a large set of microarray experiments. We classify the probes according to their signal to noise ($`S/N`$) ratio, defined as the eccentricity of a (PM, MM) pair’s ‘trajectory’ across many experiments. Of those probes having large $`S/N`$ ($`>3`$) only a fraction behave consistently with the commonly assumed hybridization model. Our results imply that the physics of DNA hybridization in microarrays is more complex than expected, and they suggest new ways of constructing estimators for the target RNA concentration.

Interest in the detailed physics of DNA hybridization is rooted in both purely theoretical and practical reasons. Studies of the denaturing transition started with models of perfectly homogeneous DNA PB1 , soon followed by studies of sequence-specific disorder HWA ; ZHANG ; NELSON . The specificity with which DNA binds to its exact complement as opposed to a mismatched copy (a “defect”) has been studied experimentally DEFECT ; BONNET  and theoretically SALERNO ; SINGH ; WATTIS . In this context it has been found that a fair fraction of the energetics of DNA hybridization is related to stacking interactions between first-neighbor bases, in addition to the obvious strand-strand contact STACK ; STACK2 . In this paper we present a study of mismatch hybridization stemming from a very practical problem, hybridization in DNA microarrays. We shall show experimental evidence that the system displays behavior which appears to be hard to account for on the basis of the extant view of hybridization specificities.

DNA microarrays provide an experimental technique for measuring thousands of individual mRNA concentrations present in a given target mixture. They are made by depositing DNA oligonucleotide sequences (probes) at specific locations on solid substrates. The probes can be either pre-made sequences as in cDNA spotted arrays, or they can be grown *in situ*, letter by letter, as in high-density oligonucleotide arrays GC1 . The target mRNA is amplified (into either cDNA or cRNA depending on the protocol) and the product labeled fluorecently before being hybridized onto the array. The spatial distribution of fluorescence is then measured using a laser, providing estimates for the target concentrations. In GeneChip arrays, the synthesis of probe sequences by photolithographic techniques requires a number of different masks and deposition processes per added base, making it impractical to grow more than a few dozen nucleotides. For such lengths, hybridization specificity is not expected to be high enough. To solve this conundrum, GeneChip technology is based on a two-fold approach, involving redundancy and differential signal GC1 ; GC2 ; GC3 . First, several different sequence snippets (each 25 bases long) are used to probe a single transcript; and second, each of these probes comes in two flavors. The perfect match (PM) is perfectly complementary to a portion of the target sequence whereas the single mismatch (MM) carries a substitution to the complementary base at its middle (13th) position. The rationale behind MM sequences is that they are expected to probe for non-specific hybridization in a manner that we shall detail below.

In current incarnations of the chips, each gene is probed by 14-20 (PM,MM) pairs (a probeset), and the task is therefore to reconstruct a single number (the target concentration) from these 28-40 measurements. There are conceivably many ways in which this can be done, with various degrees of noise rejection. The standard algorithm provided in the software suite AFFYSOFT  offers one method. However, as independent experimental techniques for measuring mRNA concentrations (like northern blots) provided clues that the analysis process should be improved upon, many researchers attempted to do so; it was then discovered that a fair number of MM probes consistently report higher fluorescence signal than their PM counterpart US . This observation is most intriguing because it violates the standard hybridization model as outlined below. Thus, the notion that the specific binding signal alone can be obtained as a differential of the PM and MM signals appears to fail in a subset of the probes.

We shall show below that it is not a matter of a few stray probes, by carefully examining the statistics of PM-MM pairs. These statistics show that *most* of the probes misbehave to various degrees. Only a fraction of them, having $`MM>PM`$, exhibit a flagrant violation of the basic assumptions, however, these just point most obviously at the symptoms of a deeper problem which affects all probes. Given the number of laboratories who are currently carrying out such hybridization experiments, squeezing out even a meager extra bit of signal to noise ratio from the data would be very valuable. It has become clear that this shall not happen in the absence of a better understanding of DNA hybridization to slightly mismatched templates. We shall now attempt the first step toward this goal, which is to characterize the problem.

The rationale behind the use of MM probes is contained in the standard hybridization model WONG :

$`PM`$ $`=`$ $`I_S+I_{NS}+B`$
$`MM`$ $`=`$ $`(1\alpha )I_S+I_{NS}+B`$
$`PMMM`$ $`=`$ $`\alpha I_S`$
Here $`PM`$ ($`MM`$) are the measured brightness of the PM (MM) probe, $`I_S`$ the contribution from specific complementary binding, $`I_{NS}`$ the amount from nonspecific binding assumed to be insensitive to the substitution, and $`B`$ a background of physical origin, i.e. the photodetector dark current or light reflections from the scanning process. Then $`\alpha `$ is the reduction of specific binding due to the single mismatch. These brightnesses are related to the quantity of interest (the RNA concentration in the sample) through:

$`I_S`$ $`=`$ $`k[S]`$
$`I_{NS}`$ $`=`$ $`h[NS]`$
where $`[S]`$ denotes the concentration of target RNA, $`[NS]`$ the concentration of whatever mixture contributes to the nonspecific hybridization. $`k`$ and $`h`$ are probe dependent specific and nonspecific susceptibilities (possibly concentration dependent) and include effects such as the areal density of probe, various affinities, transcript length dependent effects (longer transcripts are likely to carry more fluorophors depending on the labeling technique).

While it is no secret that the physics of hybridization is way more complex than this simplistic model, one could still hope that it would essentially provide a correct picture of GeneChip hybridizations. To summarize, let us outline the basic assumptions made so far: (i) non-specific binding is identical in PM and MM, meaning that $`I_{NS}`$ does not see the letter change; (ii) $`\alpha >0`$; (iii) $`k`$ and $`h`$ identical for PM and MM; (iv) $`k`$, $`h`$ and $`\alpha `$ are reasonably uniform numbers across a probe set. Notice that (i)+(ii) imply that $`PM>MM`$ always (see below).

If $`PMMM`$ is not used as such, the background $`B`$ needs to be subtracted from the intensities, which can be done in a statistically proper way as described in US .

According to the basic tenets of the standard model, it follows that $`PM>MM`$ for all probe pairs if the target RNA extract contains no sequences matching exactly the MM. In reality, one observes a vast number of probe pairs for which this assumption is violated; this behavior repeats consistently for a broad range of conditions. Our experience is that most people in the know think of this problem in terms of an imperfect adherence to the standard model, or a bothersome deviation from an otherwise properly behaving norm. In other words, the way this problem is usually characterized is ”there’s a number of probe pairs that don’t work and we don’t understand why”. We shall show now that this is not so: the $`MM>PM`$ pairs are so abundant that we like to propose the alternate view that the model is simply inadequate for describing what actually happens, and that we do not understand the basic physics of MM hybridization. Table 1 summarizes the statistics for various chip series.

The human HG-U95A chip series, for instance, has ~ 400K probes for 12K different probesets. Across a wide variety of conditions, we have observed approximately $`30\%`$ of all probe pairs have $`MM>PM`$. This figure, by itself enormous, would be easy to dismiss if most of them were in the low intensity range, where noise is expected to be relatively higher and could conceivably be bigger than $`|PMMM|`$, or if they were clustered in a small set of problematic probesets. Neither is true: $`91\%`$ of all probesets have at least 1 probe pair with $`MM>PM`$, and still $`60\%`$ of probesets have 5 such probe pairs out of 16. In addition, the $`MM>PM`$ pairs are fairly distributed with respect to brightness (cf. Fig. 1).

What could conceivably be the source for observing $`MM>PM`$? A perplexing extra bit of information lies in a simple statistic, the joint probability distribution $`P(\mathrm{log}PM,\mathrm{log}MM)`$. According to the standard model,

$$\frac{PM}{MM}=\frac{I_S+I_{NS}+B}{(1\alpha )I_S+I_{NS}+B}$$

If $`I_s`$ dominates over $`I_{NS}+B`$ then $`PM/MM1/(1\alpha )`$, while if $`I_s`$ vanishes (as when the transcript is just not present in the sample) then $`PM/MM1`$.

Thus we expect

$$1\frac{PM}{MM}\frac{1}{1\alpha }$$

Thus, the standard model predicts that $`P(\mathrm{log}PM,\mathrm{log}MM)`$ should be supported in a band, with lower limit corresponding to the diagonal $`PM=MM`$ when cross-hybridization dominates, and with an upper limit given by $`MM=(1\alpha )PM`$ for fully specific binding. Naively one would further assume that for low brightness most of the signal comes from nonspecific binding, while most would come from specific binding for high brightness. Fig. 1 shows something quite otherwise: as brightness increases, the joint probability distribution forks into two branches. The crest of the lower one lies fully below the $`MM=PM`$ diagonal.

The characteristic shapes of $`P(\mathrm{log}PM,\mathrm{log}MM)`$ are likely signatures of sequence-dependent effects. However, any hypothesis is impossible to verify as the probe sequences are not released to the public. Nevertheless, there are some obvious suspects. First, the nontrivial susceptibilities $`k`$ and $`h`$ mentioned above depend on the areal density of probe, which is sequence-dependent by virtue of the varying efficiencies of the lithography process. Secondly, nucleic acids need to unstack the single-stranded probes in order to form each new duplex as they hybridize. Further, stacking energies are extremely sensitive to sequence details, which might result in large energy barriers. This would translate into kinetics constants that vary exponentially (à la Arrhenius) in these energies, and lead to important consequences as the hybridization reactions are not carried to full thermodynamic equilibrium.

Given a set of $`N`$ experiments, further insight can be obtained by following a pair $`\stackrel{}{P}_i=(\mathrm{log}PM_i,\mathrm{log}MM_i)`$ with $`i=1,\mathrm{},N`$ *across* the entire dataset (after subtracting B). Ideally, these points would fall on a curve parametrizable by the mRNA concentration. In reality, however, the observed patterns range from nearly one-dimensional to almost circular clouds. To classify probes, we computed the center of mass CM and inertia tensor $``$ of the set of points $`\{\stackrel{}{P}_i\}`$. The positive eigenvalues of $``$, $`I_1I_2`$ define the eccentricity $`e=\sqrt{I_1/I_2}`$ and largest excursion $`\lambda _1=\sqrt{I_1}`$. Pairs with high eccentricities are those carrying high S/N, whereas $`e1`$ characterizes a very noisy probe pair.

Fig. 2 illustrate the distribution of center of mass after different filtering for $`e`$ and $`\lambda _1`$. It turns out that Fig. 2a looks very similar to Fig. 1, which is not a priori evident. On the contrary, this similarity emphasizes that most probes behave in a very reproducible manner. For instance, probes lying below the $`PM=MM`$ diagonal at the high-intensity end do so in essentially all of the 86 experiments (leading to a CM that is also below the diagonal), instead of visiting different regions of the plot. Another striking result is that (i) selecting for $`e>3`$ eliminates most of the low-intensity probes (Fig. 2b), (ii) the remaining set contains two components: one consisting of the large $`\lambda _1`$ probes (Fig. 2c) lying mostly in the $`PM>MM`$ region; and the small $`\lambda _2`$ component forming an almost perfectly symmetric “tulip” structure (Fig. 2d), containing two forked branches plus the button described in Fig. 1.

Another troubling effect which deeply affects attempts at analysis is the very broad brightness distributions within probes belonging to the same gene. Fig. 4 shows that the $`PM`$ probe intensities span up to four decades. Possible reasons for such behavior are again sequence specific effect similar to those discussed in the context of the MM behavior.

The main practical challenge is reconstructing the target mRNA concentration from the probeset data. As we showed, the variability in the hybridization properties of the probes is larger than naively anticipated, therefore, it is unlikely that a single definitive procedure will be appropriate in all cases. On the contrary, it is desirable to have several analysis tools at hand for viewing the data from different angles. For instance, as a consequence of the strongly probe dependent susceptibilities $`\alpha ,k`$ and $`h`$, the differential $`PMMM`$ will not consistently be a good estimator of the true signal. Given the unclear information contained in the $`MM`$, one alternative we studied is not considering them at all. The mRNA expression level is then obtained from a robust geometric average of the PM-B values, after a careful estimation of B US . The use of geometric averages (rather that arithmetic) is dictated by the distributions in Fig. 4. Of course, using only PM probes neglects cross-hybridization effects that would be detectable by a working MM probe, and hence tends to be less sensitive at the low-intensity end. One the other hand, it allows to rescue probesets with a high number of misbehaving MMs.

A completely different approach, closer in spirit to the model-based method WONG , would be to extend the ellipsoid of inertia idea to the the full probeset. Concretely, one would take the matrix $`A_i^j=(\mathrm{log}PM_i^j,\mathrm{log}MM_i^j)`$ ($`j=1,\mathrm{},N_p`$ is the probe and $`i`$ the experiment index) and do a principal component analysis to identify the modes carrying the most signal. After singular value decomposition $`\widehat{A}=U\mathrm{\Lambda }V^T`$, where $`\widehat{A}_i^j=A_i^jm^j`$ and $`m^j=\frac{1}{N}_iA_i^j`$ is the center of mass, the signal $`s_i=_j(m_j+\widehat{A}_i^j)V_j^1`$ is given by the projection onto the largest direction of variation. A signal-to-noise measure for the entire probeset is given by $`\frac{S}{N}=\frac{\lambda _1}{\sqrt{_{j=2}^{N_p}\lambda _j^2}}`$. Preliminary testing of the method has lead to very promising results.

In conclusion, we showed that the hybridization of short length DNA sequences to single mismatched templates exhibits a far more diverse picture than what is usually assumed. These observations do not only point at interesting physics in the DNA hybridization process to short sequences with defects, attached to a glass surface; they also have strong consequences for designers of GeneChip analysis tools, especially when it comes to the level of noise rejection of different methods. We hope this will bolster interest in the physics of hybridization and mismatch characterization.

We’d like to thank E. van Nimwegen, E. Siggia, and S. Bekiranov for sharing interesting ideas. MM acknowledges support of the Meyer Foundation; FN is a Bristol-Myers Squibb Fellow in Basic Neurosciences and acknowledges support from the Swiss National Science Foundation.


format: text

text:
# Viscoelastic Analogy for the Acceleration and Collimation of Astrophysical Jets

## 1 Introduction

Viscoelastic behavior is characteristic of fluids such as polymer suspensions in which long molecular chains become entangled. The deformation and disruption of this entanglement that occurs when the fluid is subjected to shear, and the relaxation to a new entangled state, give such fluids a “memory” that ordinary Newtonian fluids do not have. Egg whites are one common household example. These fluids are in a sense midway between being elastic solids and Newtonian fluids, being more one or the other depending upon the ratio of the relaxation time of the entangled polymers to the deformation time scale. This ratio is known as the Weissenberg number $`We`$. Viscoelastic behavior can cause flows that are qualitatively quite different from the Newtonian case, such as the Weissenberg effect, in which a fluid climbs a spinning rod.

We propose that the effects of “mesoscale” magnetic fields in some instances of MHD turbulence can be modeled by additional viscous terms, and that the structure of these terms is similar to that of a viscoelastic fluid. (It has come to our attention that this has also been conjectured in a recent paper by Ogilvie Ogilvie (2001), who demonstrated that a viscoelastic model of MHD turbulence can dramatically affect the dynamics of eccentric accretion disks.) It has long been noted that the dynamics of vortex filaments of fully-developed turbulence bear a superficial likeness to a network of polymers; see Chorin (1988). There is also an obvious analogy between the tangled polymers of a polymer fluid and the tangled field lines in MHD turbulence in the case that most of the field energy is on intermediate to small scales. However, in contrast with polymers which are typically conserved, the magnetic energy on small scales is not conserved and may dissipate. We therefore restrict our attention to flows in which the rate of shear changes slowly with respect to the dissipation time scale. We assume that the large-scale (i.e. mean-field) magnetic field is rather weak, but the magnetic field on intermediate scales is advected by hydrodynamic turbulence Ott (1998). The particular case we consider is turbulence driven by shear, such as in a Balbus-Hawley magnetic shearing instability Balbus & Hawley (1998) in the case that the net magnetic flux is negligible.

One of the more startling effects seen in viscoelastic fluids is the reversal of the secondary flow around a rotating sphere at moderate rotation speeds. A Newtonian fluid will flow inwards along the rotation axis and be expelled in the equatorial plane, but in a viscoelastic fluid, the secondary flow is exactly the opposite. We propose this phenomenon as a model for the production of astrophysical jets powered in part by the angular momentum of the central object.

## 2 MHD “rheology”

In MHD, the magnetic field evolves according to the vector advection-diffusion equation

$$_t𝐁=𝐁𝐯𝐯𝐁𝐁(𝐯)+\eta ^2𝐁,$$
(1)
where $`𝐁`$ is the magnetic field, $`𝐯`$ is the fluid velocity and $`\eta `$ is the magnetic diffusivity, which is determined by the conductivity of the fluid. If the diffusivity $`\eta `$ is zero, the flux-freezing is perfect. The contribution of magnetic fields to the stress tensor $`\tau _{ij}`$ is given by

$$\tau _{ij}=\left(\frac{𝐁^2}{8\pi }\right)\delta _{ij}+\frac{1}{4\pi }B_iB_j.$$
(2)
Let us average spatially over some “subgrid” length scale, and concentrate on the curvature term. Under perfect flux-freezing, the tensor $`M_{ij}=B_iB_j`$ can be shown from eq. (1) to evolve according to

$$D_t(M_{ij})M_{ik}L_{kj}L_{ik}^TM_{kj}+2M_{ij}L_{kk}=0$$
(3)
where $`L_{ij}=_iv_j`$ and $`D_t`$ is the usual advective derivative. The first three terms in eq. (3) are known collectively in rheology as the upper-convected invariant derivative $`𝒟_t^{[u]}`$ of the tensor $`M_{ij}`$. Insofar as the viscous stress in a non-Newtonian fluid does not depend only on the velocity gradient $`L_{ij}`$, but on its history as well as higher-order derivatives of the velocity, the concept of the advective derivative $`_t+𝐯`$ must be extended to respect the material frame indifference of the stress tensor. The upper-convected invariant derivative $`𝒟_t^{[u]}`$ of the viscous stress tensor is a common material derivative in rheology, and it is used in many rheological models.

According to flux-freezing, as expressed in eq. (3), on- and off-diagonal components of the magnetic stress in a highly conducting plasma threaded with a small-scale magnetic field will grow in response to large-scale shear. A simple analogy is a square piece of rubber on which we have drawn a tangled, statistically homogeneous and isotropic scribble. If we pull the upper edge of this sheet to the right, say, and the lower edge of the sheet to the left, the area of the square will remain the same, but the scribble is now no longer isotropic. There is an overall tendency of the scribble to be composed of lines oriented diagonally from the lower left corner to the upper right corner.

With this picture in mind, suppose that at some time $`t=0`$ the individual components of the magnetic field are distributed as Gaussian random variables with equal variance $`\sigma `$ and zero mean. In a Cartesian reference frame, given a uniform shear $`\gamma `$ in the $`\widehat{x}`$ direction, with perfect flux-freezing, the mean field will evolve according to

$$M_{ij}=\sigma ^2\left(\begin{array}{ccc}1+(\gamma t)^2& \gamma t& 0\\ \gamma t& 1& 0\\ 0& 0& 1\end{array}\right).$$
(4)
The field lines are in reality constantly shifting and being arranged on small scales by small turbulent eddies, which we characterize by a relaxation time $`s`$. This time incorporates the effects of the turbulent stretching and twisting that generates the field, as well as the dissipative processes that are governed by the magnetic diffusivity $`\eta `$. In the case we consider, the magnetic field draws its energy from the shear, so it is natural to assume that $`1/s`$ is comparable to the rate of shear $`\gamma `$. We therefore expect that the product $`\gamma s`$, which is the Weissenberg number $`(We)`$, is of the order of magnitude of one or greater. (Ogilvie finds $`We1`$ based on a comparison with simulations.) If the turbulence is driven by other processes, such as convection, the effective Weissenberg number might be very different from $`1`$.

In the case of steady shear, we approximate the form of the tensor $`B_iB_j`$ by replacing $`t`$ in the above expression with the relaxation time $`s`$. The off-diagonal stress is proportional to $`\gamma s=(We)`$, with the constant of proportionality determined by the viscosity prescription. There is a normal stress difference in the direction of shear $`\tau _{xx}\tau _{yy}=(We)\tau _{xy}`$. This form of the stress tensor, modulo the magnetic pressure, is very similar to the second-order approximation to a viscoelastic Weissenberg fluid, with zero second normal stress difference $`\tau _{yy}\tau _{zz}`$. The non-Newtonian viscous stress in such a fluid may be written in terms of the Rivlin-Erikson tensors $`A^{[1]}`$ and $`A^{[2]}`$ as $`\tau =\mu A^{[1]}+\alpha _1A^{[2]}+\alpha _2(A^{[2]})^2`$ where $`\mu `$ is the usual viscosity and the second-order coefficients $`\alpha _1`$ and $`\alpha _2`$ satisfy $`\alpha _2=2\alpha _1`$. The first Rivlin-Erikson tensor $`A^{[1]}`$ is simply the symmetrized velocity gradient, $`A_{ij}^{[1]}=L_{ij}+L_{ji}`$, and subsequent tensors $`A^{[n]}`$ obey the recursion relation $`A_{ij}^{[n+1]}=D_tA^{[n]}+A_{ik}^{[n]}L_{kj}+A_{jk}^{[n]}L_{ki}`$. The $`n`$th Rivlin-Erikson tensor gives the $`n`$th time derivative of the length of an infinitesimal material line $`|d𝐱|`$ embedded in the fluid, so that $`d_t^{(n)}|d𝐱|^2=A_{ij}^{[n]}dx_idx_j`$. A second-order fluid is thus simply an approximation to a non-Newtonian fluid that is second-order in the velocity gradient, having contributions to the viscous stress that are proportional to the second spatial derivative of the velocity as well as the square of the first spatial derivative of the velocity, in addition to the normal Newtonian viscous stress which is proportional to the (symmetrized) first spatial derivative of the velocity. As has been noted by Böhme (1987), the term “second-order fluid” is thus somewhat unfortunate since it is the flow, and not the fluid, which is really second-order. If the material derivative of the rate-of-strain tensor itself changes on a time scale comparable or faster than $`s`$ then a more sophisticated model must be brought to bear.

In a cylindrical coordinate system with azimuthal shear, such as in the fluid surrounding a rotating sphere, the normal stress difference contributes to a global hoop-stress $`\tau _{\theta \theta }`$, which acts like a rubber-band, creating a radially-inward “non-Newtonian” specific force $`f_{\mathrm{nn}}=\tau _{\theta \theta }/(r\rho )`$. In the laboratory, the hoop-stress-induced forces cause a radially-inward flow close to the equatorial plane Giesekus (1963); Geisekus (1994). This has also been shown mathematically up to various orders in $`We`$ and $`Re`$ for a second-order fluid Giesekus (1963); Böhme (1987); Joseph (1990); Geisekus (1994) as well as to leading order for a viscoelastic Oldroyd model B’ fluid Thomas & Walters (1964) which has a spectrum of relaxation times. This inflow pushes and collimates the intervening material into bipolar axial outflows for moderate $`We`$ and Reynolds number $`Re`$.

## 3 Jet model

The astrophysical analog we propose is a spinning massive central object immersed in a viscous thick disk with a tangled magnetic field. The solutions in the rheology literature do not map directly onto astrophysical solutions because gravitational, compressibility, and radiative effects are important in the latter. Furthermore, the viscosity in the astrophysical case is probably not constant as we assume here. We therefore confine ourselves to general dimensional arguments. We will allow the central object to have an extended magnetosphere that acts like a rigid rotator, so that the effective radius of the object $`R_0`$ may in fact be much larger than the actual radius of the object, depending upon the strength of the magnetosphere of the object. We assume that the disk and central object are viscously coupled so that $`\mathrm{\Omega }_0=\mathrm{\Omega }(R_0)`$, and we define $`\stackrel{~}{\omega }`$ to be the angular speed of the central object $`\mathrm{\Omega }_0`$ normalized to the Keplerian frequency at the object’s effective surface $`\mathrm{\Omega }_K(R_0)`$. Note that we may have $`\stackrel{~}{\omega }>1`$. The central object rotates in the same direction as the disk. Also we define $`r=R/R_0`$. Define $`r_m`$ to be the dimensionless radius of the thick disk, and let the angular frequency of the disk $`\mathrm{\Omega }`$ scale as $`R^q`$. Our conclusions are not sensitive to this last assumption, but it provides a nice framework to discuss the viscoelastic effects. Note that if we take $`q>3/2`$ for $`\stackrel{~}{\omega }>1`$ and $`q<3/2`$ for $`\stackrel{~}{\omega }<1`$, the thick disk will match to a Keplerian disk at $`r^{q3/2}=\stackrel{~}{\omega }`$

The structure of thick disks is sensitive to the viscosity prescription, and the “correct” viscosity is unknown. We assume that $`R_0<H<r_mR_0`$ where $`H`$ is the thickness of the thick disk. As an upper and lower bound, we set

$$\nu =\alpha \beta R_0^2\left(\frac{\mathrm{\Omega }}{\mathrm{ln}R}\right)_0.$$
(5)
where $`\beta `$ is between $`1`$ and $`r_m^2`$ and $`\alpha 1`$ as usual. The traditional parameter $`\alpha `$ is thus redundant here. The effective Reynolds number as determined by $`R_0`$ and $`\mathrm{\Omega }_0`$ is then $`Re=1/(q\alpha \beta )`$, which is on the order of $`1`$ or less.

The viscous stress is $`\tau _{r\theta }=\rho \nu (R\mathrm{\Omega }^{})`$. The specific inward “non-Newtonian” radial force $`f_{nn}`$ due to the normal stress difference is then

$$f_{nn}=(We)\nu |\mathrm{\Omega }^{}|=(We)q^2\alpha \beta R_0^2\left(\frac{\mathrm{\Omega }^2}{R}\right).$$
(6)
It must not be forgotten that the viscous and hoop stresses occur in proportion to the magnetic pressure. From eqs. (2) and (4), and recalling that $`\tau _{r\theta }=\nu (R\mathrm{\Omega }^{})`$, we find

$$\frac{1}{\rho }P_{\mathrm{mag}}=\frac{1}{8\pi }\nu |R\mathrm{\Omega }^{}|\left(3We^1+We\right)=\frac{C_1}{8\pi }\nu |R\mathrm{\Omega }^{}|.$$
(7)
The specific centrifugal inertial force is $`f_{\mathrm{inert}}=\mathrm{\Omega }^2R`$ and the specific gravitational force is $`f_{\mathrm{grav}}=\mathrm{\Omega }_K^2R`$, so that the inward hoop stresses will dominate inertial forces for $`r^2<(We)q^2\alpha \beta `$, and they will dominate gravitational forces for $`r^{(2q1)}<(We)q^2\alpha \beta \stackrel{~}{\omega }^2`$. However, it is important to note that, in the incompressible case, or for a polytrope, the specific forces $`f_{\mathrm{grav}}`$, $`f_{\mathrm{inert}}`$, and $`f_{\mathrm{press}}`$ are all derivable from potentials, whereas $`f_{\mathrm{nn}}`$ is not. Note that the specific pressure force $`f_{\mathrm{press}}=\rho ^1P`$ may be written as the gradient of some function so long as there is a one-to-one relation between $`P`$ and $`\rho `$. We expect that there is typically no realistic combination of gravitational, centrifugal, or pressure forces that will balance the non-Newtonian hoop stresses for an arbitrary choice of $`v_\theta (R,z)`$, regardless of the value of $`\beta `$. For the choice of $`_zv_\theta =0`$, such as in a thin disk, the hoop stress $`f_{\mathrm{nn}}`$ depends only on $`R`$ and thus may be written $`f_{\mathrm{nn}}=\mathrm{\Phi }`$ for some $`\mathrm{\Phi }`$, but for a thick disk this will not typically be the case.

In contrast with azimuthal viscous forces which dissipate kinetic energy, $`\overline{f}_{\mathrm{nn}}\overline{v}`$ is positive for a radial inflow, so that the hoop stresses do positive work on the fluid. There are no hoop stresses along the rotation axis of the system, so the energy gained by the fluid is free to escape in this direction. Lastly, as in the laboratory case, it is hypothesized that these hoop stresses naturally confine and collimate the flow, if the thick disk is sufficiently larger than the object radius $`R_0`$.

The hoop stresses are powered in part by the mechanical accretion luminosity, and the specific work performed by them for a fluid element in the equatorial plane that moves from the outer radius $`R_m`$ to $`R_0`$ is

$$W=_{R_0}^{R_m}f_{\mathrm{nn}}𝑑R=\frac{1}{2}q\alpha \beta (We)\frac{GM}{R_0}(1r_m^{2q})$$
(8)
For $`\beta 1`$ this work can be much larger than the binding energy; some of this energy will be viscously dissipated, and some is work against the pressure gradient.

Extrapolating from the laboratory case for moderate $`We`$ and $`Re`$, the launch speed of the jet is expected to be of the order of

$$v_{\mathrm{jet}}=C_2\frac{\nu }{R_0}\left(\frac{\mathrm{\Omega }_0R_0^2}{\nu }\right)^2=C_2\frac{R_0\mathrm{\Omega }_0}{q\alpha \beta },$$
(9)
where $`C_2`$ is a constant of order unity. In escaping to infinity, the jet will lose its binding energy $`GM/R_0=(R_0\mathrm{\Omega }_0/\stackrel{~}{\omega })^2`$, but it will be propelled by the magnetic pressure $`P_{\mathrm{mag}}/\rho `$ given in eq. (7). If the thick disk is substantially thicker than the object effecive radius $`R_0`$ then the jet may also suffer some viscous drag, which we estimate by Poiseuille’s law. In traversing the thickness of the disk $`H`$ a jet of radius $`aR_0`$ (where we expect $`a1`$) will lose an amount of energy per unit mass of the order of $`\nu v_{\mathrm{jet}}H/(aR_0)^2=C_3HR_0\mathrm{\Omega }_0^2`$ where $`C_3`$ is of order unity. We thus very roughly estimate the jet speed at infinity as

$$\frac{1}{2}v_{\mathrm{}}^2R_0^2\mathrm{\Omega }_0^2\left(\frac{C_1q^2\alpha \beta }{8\pi }+\frac{C_2}{2(q\alpha \beta )^2}\frac{C_3H}{R_0}\frac{1}{\stackrel{~}{\omega }}\right).$$
(10)
The speed of the jet naturally scales as $`R_0\mathrm{\Omega }_0`$, and increases with the speed of rotation $`\stackrel{~}{\omega }`$, which must be close to breakup for $`\beta 1`$ for the jet to reach infinity. We caution however that further work is required to find to what degree the above admittedly very approximate relationship holds in the astrophysical case. This expression is intended only to show roughly the relative importance of the various effects we have considered. We have not considered any relativistic effects, and we have not addressed the importance of the accretion efficiency, nor any effects associated with the central object having an event horizon. Accretion onto a black hole will certainly complicate matters very much. As it stands, therefore, we believe our present arguments apply perhaps best to the case of jets from protostellar objects.

## 4 Discussion

We have proposed that the fluid dynamics of a highly-conductive turbulent medium with tangled magentic fields is similar in some respects to the viscoelastic fluid dynamics of a suspension of tangled polymers. We discuss the production of jet-like structures in viscoelastic fluids induced by a spinning sphere, and outline an analogous model for the production of astrophysical jets propelled by a spinning central object. The processes outlined above are thus hypothesized to propel jets by tangled magnetic fields, as has been previously suggested by Heinz & Begelman (2000).

We believe that this tentative model for the acceleration of jets has some advantages over previous models. Firstly, the jet process does not rely upon the insertion “by hand” of any large-scale magnetic field. The jet also appears to be self-collimating, and it is a solid-core jet, rather than a hollow-core jet. As pointed out by Li (2001), tangled magnetic fields may also stabilize the jet from the long-wave mode kink instability that plagues standard MHD jet models. Lastly, the jet mechanism we describe is powered in part by the angular momentum of the central object. We hope to improve upon this work with more detailed calculations in a future paper.

The majority of this work was completed when the author was a graduate student at the University of Texas at Austin. The author wishes to thank Craig Wheeler for helpful comments, John Scalo and Nairn Baliber for encouraging him over the years to write this up, and Robert Buchler for giving him the time and financial support necessary to do so.


text:
# ASCA Observation of the Nearest Gravitational Lensing Cluster Candidate – A 3408

## 1. Introduction

Clusters of galaxies are the largest gravitationally bound systems in the universe, which set clear constraints on the formation of the structure and composition of the universe. Furthermore, the gravitational mass of clusters is an important observational quantity to constrain the density parameter, $`\mathrm{\Omega }`$.

There are primarily two methods to estimate the total (gravitational) mass of clusters of galaxies. One method employs X-ray observations under the assumption of hydrostatic equilibrium of the hot gas. The other uses the configuration of gravitational lensing arcs seen in optical observations. However, the mass estimates from lensing observations have been found to be systematically about 2-times larger than the mass estimates from X-ray observations (Wu, Fung 1997). Several ideas have been proposed to solve this discrepancy, including the following: (1) uncertainty in selecting the cluster center (Allen 1998; Hashimotodani 1999); (2) effect of the cooling flow (Allen 1998); and (3) mass profiles different from the conventional $`\beta `$ model, e.g., the NFW density profile (Navarro et al. 1996).

In order to examine these possibilities, detailed X-ray observations of the central region of clusters are important. Such observations may be realized either with higher resolution instruments or equivalently by taking low-z gravitational lensing clusters as targets. This is our motivation to observe low-z gravitational lensing clusters.

Recently, several groups have found gravitational lensing clusters with low-z. For instance, Allen, Fabian, and Kneib (1996) discovered a z = 0.43 lensed arc in the massive cooling flow cluster PKS 0745–191 at z = 0.103. Blakeslee and Metzger (1999) discovered a lensed arc at z = 0.573 in the nearby cD cluster A 2124 at z = 0.066. However, the lensing properties of these low-z gravitational lensing clusters remain unclear. Especially, detailed X-ray observations of low-z gravitational lensing clusters have rarely been performed. A 3408 at z = 0.042 is the lowest-z gravitational lensing cluster candidate in which arclike objects are found. In this paper, we present the results of an X-ray observation of A 3408 with ASCA.

For our study, we adopt $`\mathrm{\Omega }_0=1`$, $`\lambda =0`$, and $`H_0=50h_{50}^1`$ km s-1 Mpc-1; 1 corresponds to 68 $`h_{50}^1`$ kpc at the distance of the cluster.

## 2. Discovery of an Arc-like Object in A 3408

Because A 3408 has a relatively low galactic latitude ($`b=17.57^{}`$), detailed optical observations have not been performed. According to the ACO catalog (Abell et al. 1989), the number of galaxies is 43 and the richness class is 0. Campusano and Hardy (1996) discovered an arclike object near the center of A 3408. Figure 1 displays the B \+ I band image of A 3408. An arc-like object is $`50^{\prime \prime }`$ away from the brightest cluster galaxy (BCG). The length of the arc-like object is $`11^{\prime \prime }`$, but its width is unresolved. The central elliptical galaxy has a redshift of 0.042, and the arc-like object has a redshift of 0.073. If the arc-like object is the lensed image of a background galaxy, A 3408 is the nearest gravitational lensing cluster that has been observed.

## 3. Lens Models of A 3408

Campusano et al. (1998) considered two lensing mass models to evaluate the projected mass within the arc radius, regarding the BCG as the center of the potential of A 3408.
The minimum mass case (The BCG model):
The mass distribution follows the light profile of the central brightest elliptical galaxy. The cluster component was not taken into account. The mass enclosed within $`r_{\mathrm{arc}}`$ (from the BCG to the arc) is $`M(r_{\mathrm{arc}})1.5\times 10^{12}M_{}`$, and the amplification of the source galaxy is $`\mathrm{\Delta }m=0.01`$.
The maximum mass case (The dark-halo model):
A typical massive dark-halo of cluster scale is assumed. They adopted a massive mass profile comparable to that found in A 2218. The mass enclosed within $`r_{\mathrm{arc}}`$ is $`M(r_{\mathrm{arc}})4.4\times 10^{13}M_{}`$, and the amplification of the source galaxy is $`\mathrm{\Delta }m=0.60`$.

Note that these models are based almost solely on a single arc image and that no X-ray information was employed. They could not determine which model is appropriate for the arc-like object.

## 4. ASCA Observations of A 3408

We observed A 3408 on 1999 March 13–14 with ASCA (Tanaka et al. 1994). With the exception of the ROSAT all-sky survey, no other X-ray observations had been carried out for this cluster. The data from the Gas Imaging Spectrometer (GIS; Ohashi et al. 1996; Makishima et al. 1996) were taken in the PH mode, while those of the Solid-state Imaging Spectrometer (SIS) were acquired in the 2 CCD Faint mode. Because the SIS spectrum was severely distorted, possibly by a dark current, and its field of view is smaller than that of the GIS, we concentrated on the GIS data in a following analysis. The data selection was done using the standard screening parameters. We excluded the GIS data with a cutoff rigidity lower than 6 GV and those with elevation angles below 5 from the night-earth rim and 25 from the day-earth rim. After screening, the total exposure time was $``$ 50 ks for the GIS.

Figure 2 shows the 0.7–10 keV GIS image of A 3408 (contour) overlaid on an optical image from the Digitalized Sky Survey (DSS). Extended emission from A 3408 is surely detected. The sky coordinates of the GIS image are adjusted to the point source, Seyfert 1 galaxy 1H0707–495, located at $`15^{}`$ south to the center of the GIS field of view. It is found that the X-ray surface brightness peak of A 3408 is not coincident with the position of the BCG. The X-ray peak is $`60^{\prime \prime }`$ away from the BCG (PGC 070830.4–491246). Considering a GIS pointing accuracy of $`24^{\prime \prime }`$ (Gotthelf et al. 2000), we conclude that the BCG is not located at the X-ray peak, which we usually assume to be the center of the potential.

## 5. Spectral Analysis

We estimated the background spectrum by the following method. The non-X-ray background was reproduced with a method introduced by Ishisaki (1996). The cosmic X-ray background was extracted from the Large Sky Survey (Ueda et al. 1999) observed during the ASCA PV phase.

We fitted the GIS spectra with a single-temperature thermal-emission model by Raymond and Smith (1977). We fixed the hydrogen column density to the galactic absorption of $`N_\mathrm{H}=6.0\times 10^{20}`$ cm-2 (Dickey, Lockman 1990) and the redshift to $`z=0.042`$. We first extracted spectra from the region $`r<250`$ kpc and $`250`$ kpc $`<r<800`$ kpc. The fit results are summarized in table 1.

The temperature of the region $`r<250`$ kpc is similar to that of the region $`250`$ kpc $`<r<800`$ kpc. We also fitted the spectrum of the region within the arc radius ($``$110 kpc). The temperature of the region $`r<r_{\mathrm{arc}}`$ is $`3.0_{0.3}^{+0.3}`$ keV. These results imply a nearly isothermal nature for the ICM of A 3408.

We next extracted spectra from the combined region $`r<800`$ kpc. Figure 3 shows the spectrum of the region $`r<800`$ kpc. The temperature of this region is $`2.9_{0.2}^{+0.2}`$ keV and the luminosity in the 2–10 keV energy band is $`L_\mathrm{X}=2.3_{0.2}^{+0.1}\times 10^{43}`$ $`h_{50}^2`$ erg s-1. The luminosity measured by the ROSAT all-sky survey in the 0.1–2.4 keV energy band is $`L_\mathrm{X}=5\times 10^{43}`$ $`h_{50}^2`$ erg s-1 (Ebeling et al. 1996), which is consistent with our result derived from the spectral fitting ($`L_\mathrm{X}=4.5_{0.1}^{+0.3}\times 10^{43}`$ $`h_{50}^2`$ erg s-1). The temperature of the ICM corresponds to a velocity dispersion of $`674_{24}^{+22}`$ km s-1, whereas we assume equal specific energy in the gas and galaxies.

## 6. Spatial Analysis

The X-ray surface brightness profile was extracted by accumulating photons within annular bins of $`0^{}.75`$ width. We considered the X-ray peak as the center. Figure 4 shows the X-ray surface brightness profile of A 3408. We fitted the X-ray surface brightness profile with an isothermal $`\beta `$ model. In the $`\beta `$ model, the gas density profile is described as $`\rho _{\mathrm{gas}}(r)=\rho _{\mathrm{gas0}}/[1+(r/r_\mathrm{c})^2]^{\frac{3\beta }{2}}`$, where $`\beta `$ is the slope parameter of the outer part of the distribution and $`r_\mathrm{c}`$ is the core radius of the cluster. Assuming this density profile, the X-ray surface brightness profile is described as

$$S(b)=S_0\frac{1}{(1+(\frac{b}{r_\mathrm{c}})^2)^{3\beta 1/2}}+\mathrm{background},$$
(1)
where $`S_0`$ is the central surface brightness and $`b`$ is the projected 2D distance from the cluster center. Because of the complicated PSF of the ASCA XRT, we used a ray-tracing simulation (Ikebe et al. 1996) to fit the radial profile. The fit results are shown in figure 4 and table 2. The profile is well fitted with a single $`\beta `$ model and doesn’t exhibit a significant excess X-ray emission from the central region, which is expected for a cooling flow. Note that the core radius and $`\beta `$ are similar to the typical value of nearby clusters determined by Einstein observations (Abramopoulos, Ku 1983).

We also fitted the same profile with an isothermal gas model in the NFW potential. The gas-density profile is analytically modeled as $`\rho _{\mathrm{gas}}(r)=\rho _{\mathrm{gas0}}/[1+(r/r_\mathrm{s})]^{B\frac{r_\mathrm{s}}{r}}`$ (Makino et al. 1998). Here, $`r_\mathrm{s}`$ is the scaling parameter of the cluster, and $`B`$ is the gradient of the gas density. The fit results are shown in figure 4 and table 3. There is little difference between two best-fit curves, and the $`\chi ^2`$ value of the NFW fit is similar to that of the $`\beta `$ fit. Instruments with a better spatial resolution and a larger effective area are needed to distinguish these models. In this paper, we adopt the $`\beta `$-model parameters to derive the mass distribution.

## 7. Mass Profile of A 3408

From the $`\beta `$-model parameters and the temperature of the ICM, we calculated the mass distribution of A 3408. Assuming hydrostatic equilibrium and neglecting the temperature gradient, the mass contained within a circle of radius $`b`$ on the sky plane is given by

$$M(b)=1.14\times 10^{14}\beta \stackrel{~}{m}(b)(\frac{kT}{\mathrm{keV}})(\frac{r_\mathrm{c}}{\mathrm{Mpc}})M_{},$$
(2)
where

$$\stackrel{~}{m}(b)=\frac{R_0^3}{1+R_0^2}_{b_0}^{R_0}x\sqrt{x^2b_0^2}\frac{3+x^2}{(1+x^2)^2}𝑑x.$$
(3)
In equation (3), $`b_0=b/r_c`$ and $`R_0=R/r_\mathrm{c}`$, where $`R`$ is the physical size of the cluster, for which we assumed 5 Mpc.

In order to estimate the uncertainty in our mass calculation, we employed the 90% confidence errors in the spectral and imaging parameters. Figure 5 shows the mass profile of A 3408 integrated along the line of sight. Each of the 27 lines corresponds to a combination of either the best-fit value or the 90% confidence lower or upper limit of the three parameters: $`\beta `$, $`r_\mathrm{c}`$ and $`kT`$. The projected mass within the arc radius, $`r_{\mathrm{arc}}`$, of 110 kpc is $`M(r_{\mathrm{arc}})=1.2_{0.4}^{+0.8}\times 10^{13}M_{}`$.

We also applied polytropic models (Sarazin 1988) in which some temperature gradient in the intracluster gas is assumed. If the hot gas is strictly adiabatic, its pressure and density have a simple relation of $`P\rho ^\gamma `$. In this model, $`\gamma =1`$ implies that the gas distribution is isothermal. The hydrostatic mass contained within a sphere of a radius $`r`$ is prescribed by

$$M(r)=\frac{\gamma kT_{\mathrm{gas}}(r)}{G\mu m_\mathrm{p}}r\frac{d\mathrm{ln}\rho _{\mathrm{gas}}}{d\mathrm{ln}r},$$
(4)
where

$$T_{\mathrm{gas}}(r)=T_{0,\gamma }(\frac{\rho _{\mathrm{gas}}(r)}{\rho _0})^{\gamma 1}=T_{0,\gamma }(1+\frac{r^2}{r_\mathrm{c}^2})^{3\beta (\gamma 1)/2}.$$
(5)
Markevitch et al. (1998) found universal temperature gradients in nearby rich clusters observed with ASCA, which are well described by the polytropic distribution with $`\gamma =`$1.2 – 1.3, implying temperature increases toward the center. Therefore, we introduced a polytropic gas distribution with $`\gamma =1.3`$. Although there is no observational evidence, we also tested the case with $`\gamma =0.9`$ for a temperature decrease toward the center. However, because the temperature gradient is significant in the region $`r>r_\mathrm{c}`$, the range of the estimated mass is almost the same as that in the isothermal case, $`M(r_{\mathrm{arc}})=1.2_{0.4}^{+0.9}\times 10^{13}M_{}`$ ($`\gamma =0.9`$) and $`M(r_{\mathrm{arc}})=1.4_{0.5}^{+0.6}\times 10^{13}M_{}`$ ($`\gamma =1.3`$).

## 8. Discussion

### 8.1. Comparison with Models by Campusano et al.

The existence of hot gas extending over a few Mpc in A 3408 rejects the BCG model by Campusano et al. (1998), in which only the BCG was taken into account. The projected mass of A 3408 within the arc radius derived based on our X-ray observation is $`1.2_{0.4}^{+0.8}\times 10^{13}M_{}`$. This is 18 – 45% of the mass obtained from the dark-halo (cluster) model by Campusano et al. (1998). However, it should not be taken straightforward, since the center of the potential that they assumed, i.e., BCG, differs from the X-ray center that we defined. Furthermore, they assumed that the cluster mass profile of A 3408 is similar to that of a distant cluster, A 2218, while our mass profile is based on the X-ray observation of A 3408, itself. On these two points, we consider that our mass estimation based on the X-ray observation has some advantage over the mass estimation by Campusano et al. (1998).

Our mass model is based on the following assumptions: an isothermal gas distribution, hydrostatic equilibrium, and a spherically symmetric gas distribution. However, as mentioned in section 7, the range of the estimated mass is little affected by possible temperature gradients with $`\gamma `$=1.3 or 0.9. Meanwhile, if A 3408 is a merging or post-merging cluster, the assumptions of hydrostatic equilibrium and a spherically symmetric gas distribution would be inappropriate. However, the relatively regular X-ray surface-brightness distribution of A 3408 suggests that the cluster is a moderately relaxed system. This is also supported by the fact that the morphological classification of A 3408 is Bautz-Morgan Type I-II (Abell et al. 1989). In order to confirm that A 3408 is really a relaxed system, a further investigation, such as concerning the velocity dispersion of member galaxies, is needed.

### 8.2. Potential Model Accounting for Both the X-ray Data and the Arc Image

For a given gravitational potential (cluster mass distribution) and a source position, the X-ray emission profile and the lensing arc image can be calculated. Conversely, we constrained the source position and the gravitational potential that reproduce the arc image and the X-ray data, as was done in Hattori et al. (1998) and in Hashimotodani (1999). We performed a similar analysis for A 3408. In addition to the potential of the cluster (we allow for deviations from spherical symmetry), the potential due to the BCG was considered, since, in several clusters of galaxies, bright galaxies near the arc are considered to play a significant role in lensing (Hattori et al. 1998, Hashimotodani 1999). The potential depth of the BCG was estimated by using the Faber–Jackson relation (Faber, Jackson 1976); it is $`300`$ km s-1 in terms of the velocity dispersion. Figure 6 shows the source position predicted by our model. The source position determined by Campusano et al. (1998) in their dark halo model is shown in the same figure as a reference. The magnification factor for the source galaxy is constrained to be $`\mathrm{\Delta }m=0.07\pm 0.03`$ mag. Thus, the source galaxy is only slightly magnified by the lens A 3408 in our potential model. The small magnification factor is also consistent with the slightly distorted image of the arc-like object.

Because of the small redshift of the source galaxy for the A 3408 lens, we have a good opportunity to study its intrinsic properties and the modification by the lens in detail. The characteristics of the source galaxy were discussed by Campusano et al. (1998) under the assumptions of their minimum mass case and the maximum mass case, for which they evaluated the source magnification factors to be $`\mathrm{\Delta }m=0.01`$ and 0.60, respectively. On the other hand, we provided $`\mathrm{\Delta }m=0.07`$, which implies that the absolute magnitude, $`M_B`$, of the source galaxy is $`18.7`$. According to the galaxy population statistics by Binggeli, Sandage, and Tammann (1988), about 75% of the field galaxies with the same absolute magnitude are spiral. Considering the small magnification factor, $`\mathrm{\Delta }m=0.07`$, of the lens, the source galaxy is not significantly modified in its shape, and should be regarded as an edge-on spiral originally. In this case, except for the small probability to realize an edge-on configuration, the low rotation velocity of the source galaxy measured by Campusano et al. (1998) might imply a somewhat extraordinary nature of the source galaxy. This is because even the smallest Sc galaxies have twice the observed rotation velocity of the source galaxy, as mentioned by Campusano et al. (1998). Although further study is required concerning this point, this kind of examination is generally important for a nearby lensing object.

## 9. Summary

We observed the nearest gravitational lensing cluster candidate, A 3408, with ASCA. Our observation determined the ICM temperature and X-ray surface brightness profile of A 3408 for the first time. We also showed that the gas distribution of A 3408 is isothermal at the central region of A 3408, even within the arc radius.

We compared our mass model based on the X-ray observation with that of Campusano et al. (1998) based on the optical image. The mass derived by us is 18 – 45% of that by Campusano et al. (1998). The primary cause of this discrepancy is in their assumption that the center of the cluster potential coincides with the position of the BCG, while we reveal the BCG position to be 60′′ outside of the X-ray center.

As shown in this paper, X-ray observations are essential to study gravitational lens phenomena and to determine mass profiles in clusters of galaxies. Hashimotodani (1999) found that there is a significant number of samples in which the BCG position deviates from the X-ray peak of the cluster. He suggested that this can be one of the primary causes for the discrepancy in the mass estimates from the X-ray data and from the lensing images.

One advantage of observing low-z gravitational lensing clusters is that it enables us to study X-ray profiles even with an instrument having only moderate spatial resolution. Further detailed studies of the central part of A 3408, such as a temperature gradient and gas distribution within the arc radius, will be enabled with better spatial resolution and larger effective-area instruments, like XMM-Newton. Another advantage is that it becomes possible to examine a lensed image of the source galaxy, if the source redshift is small, as in our case. In fact, we showed that the weak shear of the arc-like object is consistent with the small magnification factor which we derived.

We are grateful to the ASCA team for operating the spacecraft and supporting the data analysis. H.K. is supported by JSPS Research Fellowship for Young Scientists. The computation code for the lens image and cluster potential used in the last subsection was provided by M. Hattori and M. Murata, Tohoku University.

## References

Abell, G. O., Corwin, H. G., Jr. & Olowin, R. P. 1989, ApJS, 70, 1

Abramopoulos, F., & Ku, W. H. -M. 1983, ApJ, 271, 446

Allen, S. W. 1998, MNRAS, 296, 392

Allen, S. W., Fabian, A. C., & Kneib, J. P. 1996, MNRAS, 279, 615

Binggeli, B., Sandage, A., & Tammann, G. A. 1988, ARA&A, 26, 509

Blakeslee, J. P., & Metzger, M. R. 1999, ApJ, 513, 592

Campusano, L. E., & Hardy, E. 1996, in IAU Symp, 173, Astrophysical applications of gravitational lensing, ed. C. S. Kochanek & J. N. Hewitt (Dordrecht: Kluwer), 125

Campusano, L. E., Kneib, J. -P., & Hardy, E. 1998, ApJ, 496, L79

Dickey, J. M., & Lockman, F. J. 1990, ARA&A, 28, 215

Ebeling, H., Voges, W., B$`\ddot{o}`$hringer, H., Edge, A. C., Huchra, J. P., & Briel, U. G. 1996, MNRAS, 281, 799

Faber, S. M., & Jackson, R. E. 1976, ApJ, 204, 668

Gotthelf, E. V., Ueda, Y., Fujimoto, R., Kii, T., & Yamaoka, K. 2000, ApJ, 543, 417

Hashimotodani, K. 1999, Ph. D. Thesis, The University of Osaka

Hattori, M., Matuzawa, H., Morikawa, K., Kneib, J. -P., Yamashita, K., Watanabe, K., B$`\ddot{o}`$hringer, H., & Tsuru, G. T. 1998, ApJ, 503, 593

Ikebe, Y. 1996, Ph. D. Thesis, The University of Tokyo (RIKEN IPCR CR-87)

Ishisaki, Y. 1996, Ph. D. Thesis, The University of Tokyo (ISAS-RN 613)

Makino, N., Sasaki, S., & Suto, Y. 1998, ApJ, 497, 555

Makishima, K., Tashiro, M., Ebisawa, K., Ezawa, H., Fukazawa, Y., Gunji, S., Hirayama, M., Idesawa, E. et al. 1996, PASJ, 48, 171

Markevitch, M., Forman, W. R., Sarazin, C. L., & Vikhlinin, A. 1998, ApJ, 503, 77

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563

Ohashi, T., Ebisawa, K., Fukazawa, Y., Hiyoshi, K., Horii, M., Ikebe, Y., Ikeda, H., Inoue, H. et al. 1996, PASJ, 48, 157

Raymond, J. C., & Smith, B. W. 1977, ApJS, 35, 419

Sarazin, C. L. 1988,“X-ray emissions from clusters of galaxies”, (Cambridge: Cambridge University press)

Tanaka, Y., Inoue, H., & Holt, S. S. 1994, PASJ, 46, L37

Ueda, Y., Takahashi, T., Inoue, H., Tsuru, T., Sakano, M., Ishisaki, Y., Ogasaka, Y., Makishima, K. et al. 1999, ApJ, 518, 656

Wu, X. -P., & Fung, L. -Z. 1997, ApJ, 483, 62

Figure Captions

Fig. 1. B \+ I band image of an arc-like object in A 3408 taken with the 0.9 m CTIO telescope (Campusano et al. 1998). The image size is $`2^{}\times 2^{}`$. An arc-like object is seen 50′′ away from the BCG.

Fig. 2. GIS image of A 3408 (contour) overlaid on an optical image from DSS (Digitalized Sky Survey). The GIS energy range is 0.7–10 keV. The GIS image was smoothed by a $`\sigma =0_.^{^{}}5`$ Gaussian filter. The spacing of the contour is linear with a step size of $`4.2\times 10^4`$ count s-1 arcmin-2. The lowest contour is at $`8.6\times 10^4`$ count s-1 arcmin-2. The circle represents the GIS pointing accuracy.

Fig. 3. GIS spectrum of the region $`r<800`$ kpc.

Fig. 4. X-ray surface brightness profile of A 3408. The solid line is the best-fit curve of the $`\beta `$-model fitting and the dotted line is that of the NFW model fitting.

Fig. 5. Projected mass distribution of A 3408.

Fig. 6. Magnifications of a background source at $`z=0.073`$ produced by our lens model of A 3408. The image size is $`100^{\prime \prime }\times 100^{\prime \prime }`$, and the center (0,0) corresponds to the BCG position. The contours indicate magnifications of $`\mathrm{\Delta }m=`$0.04, 0.05, 0.07, 0.11, and 0.19 from the outer to the inner. S1 and S2 correspond to the position of the source galaxy derived by Campusano et al. (1998) and by our model, respectively. The cross represents the X-ray peak.


text:
# Evidence for Slow Velocity Relaxation in Front Propagation in Rayleigh-Bénard Convection

## I Introduction

Although the propagation of a front into an unstable state plays an important role in various physical situations ranging from the pearling instability  to dielectric breakdown , detailed experimental tests of the explicit theoretical predications, especially those for the velocity of so-called “pulled” fronts are scarce. One of the reasons lies in the difficulty in preparing the system in the unstable state.

If the initial front profile is steep enough the propagating front converges to a unique shape and velocity. Theoretically, one distinguishes two regimes for front propagation into unstable states: the so-called “pushed” regime, where the front is driven by the nonlinearities and the so-called “pulled” regime where the asymptotic velocity of the propagating front, $`v_{as}`$, equals the spreading speed $`v^{}`$ of linear perturbations around the unstable state: $`v_{as}=v^{}`$. Pushed fronts are by definition those for which the asymptotic speed $`v_{as}`$ is larger than $`v^{}`$: $`v_{as}>v^{}`$. It is thus as if a “pulled” front is literally “pulled” by the leading edge whose dynamics is driven by linear instability of the unstable state ; the nonlinearities merely cause saturation behind the front. We focus here on the experimental tests of the dynamics of such pulled fronts; since $`v^{}`$ is determined by the equations linearized about the unstable state, the front velocity of pulled fronts can often be calculated explicitly, even for relatively complicated situations.

There have been two experiments aimed at testing the predictions for the speed of pulled fronts. Almost 20 years ago Ahlers and Cannell  studied the propagation of a vortex-front into the laminar state in rotating Taylor-Couette flow. The measured velocities were about 40% smaller than expected from the theoretical predictions. A few years later, however, Fineberg and Steinberg (FS)  published data which appeared to confirm the expected velocity in a Rayleigh-Bénard convection experiment to within about 1%. The issue then seemed to be settled when it was also shown that the discrepancy observed by Ahlers could be traced back to slow transients .

The theoretical developments of the last few years give every reason to reconsider the old experiments by FS: It has been shown  that the convergence of the velocity of pulled fronts is always very slow, in fact with leading and subleading universal terms of $`𝒪(1/t)`$ and $`𝒪(1/t^{3/2})`$ with prefactors which follow from the linearized equation. This slow relaxation implies that it will in general be very difficult to measure the asymptotic front speed to within a percent or so in any realistic experiment. Hence, from this new perspective, the proper question is not why in the Taylor-Couette experiment the measured velocity was too low, but why apparently in the Rayleigh-Bénard experiment of FS the asymptotic front speed was measured.

The main purpose of this paper is to address this issue, and to reanalyze the experiments in the light of the present theory. We will conclude that the data of FS actually do show signs of the predicted power law convergence of the front velocity to an extrapolated asymptotic value which is about 15% larger than their transient value. This of course implies that there is then a discrepancy of order 15% between the value of $`v_{as}=v^{}`$ as extrapolated from their data, and the one claimed in the original experiments. We will argue that the most likely reconciliation of the two results is that the value of the correlation length $`\xi _0`$ in the experimental cell of FS is somewhat larger than the theoretical value used by FS to interpret their data.

Of course, only new experiments can settle whether the interpretation we propose is the correct one. We do consider new experiments along the lines of FS in fact very desirable, not so much as they might settle the numerical value of the velocity, but more because they hold the promise of being the first accurate experimental test of the universal power law relaxation of pulled fronts.

In section II we will first summarize the relevant theoretical predictions for the velocity of pulled fronts. Then we will discuss the experiments of FS in the light of these results in section III, where we will also reanalyze their data. Finally, in section IV we turn to a brief discussion of the wavenumber of the pattern selected by the front. Here, the results of FS were not quite consistent with the predictions for the asymptotic wavenumber from the Swift-Hohenberg equation. As we shall discuss, the wavelength of the pattern is affected by various effects which are not easily controlled, but the most likely interpretation of the data of FS is that they did not observe the asymptotic wavelength behind the front, but the local wavenumber in the leading edge of the front. Indeed it is in general difficult to test the theory by studying the asymptotic pattern wavelength and the convergence to the asymptotic value experimentally.

## II Summary of theoretical predictions

### A Asymptotic speed and power law convergence

Just above the onset of a transition to stationary finite wavelength patterns, for small dimensionless control parameters $`ϵ`$ the slow dynamics on length scales larger than the wavelength of the pattern can be described by the Ginzburg-Landau amplitude equation

$$\tau _0_tA=ϵA+\xi _0^2_x^2Ag|A|^2A.$$
(1)
The time scale $`\tau _0`$ and length scale $`\xi _0`$ as well as the nonlinear saturation parameter $`g`$ depend on the particular system under study.

The asymptotic spreading speed $`v^{}`$ of linear perturbations around the unstable state is in general obtained from the linear dispersion relation $`\omega (k)`$ of a Fourier mode $`e^{i\omega t+ikx}`$ through

$$\frac{\mathrm{Im}\omega }{\mathrm{Im}k}|_k^{}v^{}=0,\frac{\mathrm{Im}\omega }{\mathrm{Re}k}|_k^{}=0,\frac{\mathrm{Im}\omega (k^{})}{\mathrm{Im}k}=v^{}.$$
(2)
This yields for the Ginzburg-Landau equation

$$v^{}=2ϵ^{1/2}\xi _0\tau _0^1.$$
(3)
and

$`\mu ^{}`$ $``$ $`\text{Im}k^{}=\sqrt{ϵ}/\xi _0,`$ (4)
$`D`$ $``$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{^2\mathrm{Im}\omega }{(\mathrm{Im}k)^2}}|_k^{}={\displaystyle \frac{\xi _0^2}{\tau _0}}.`$ (5)
In the experiments with which we will compare, the scaled velocity

$$\stackrel{~}{v}=v\frac{\tau _0}{\xi _0\sqrt{ϵ}}$$
(6)
is often used. According to (3), for a pulled front the asymptotic value $`\stackrel{~}{v}_{as}=2`$.

The above results were known in the eighties, at the time when the experiments were done. The crucial insight of the last few years is the finding that the convergence or relaxation towards the asymptotic velocity $`v^{}`$ of pulled fronts is always extremely slow: the general expression for the time dependent velocity $`v(t)`$ emerging from steep initial conditions (i.e., decaying faster than $`e^{\mu ^{}x}`$) is given by

$$v(t)=v^{}\frac{3}{2\mu ^{}t}+\frac{3}{2\mu ^2t^{3/2}}\sqrt{\frac{\pi }{D}}+𝒪(t^2).$$
(7)
For the case of the Ginzburg-Landau amplitude equation, we then have

$$\stackrel{~}{v}(t)=2\frac{3}{2ϵt/\tau _0}+\frac{3\sqrt{\pi }}{2(ϵt/\tau _0)^{3/2}}+\mathrm{}$$
(8)
It is important to stress that the above expression for $`v(t)`$ is exact but asymptotic — this is illustrated by the fact that at time $`ϵt/\tau _0=\pi `$ the subdominant $`t^{3/2}`$ term is equal to the first correction term of order $`t^1`$ in absolute value, but of opposite sign. Thus, although for any sufficiently long time $`t`$, the above expression will always become accurate, at any finite time, however, the expression might only yield a good estimate. In fact, in practice one usually has to go to dimensionless times $`ϵt/\tau _0`$ of order 10 or larger for the asymptotic expression to become accurate, while for dimensionless times $`\epsilon t/\tau _0`$ in the range 3 to 10 the first correction term yields a reasonable order-of-magnitude estimate . As we shall see below, in the Rayleigh-Bénard experiments  the maximum dimensionless time $`ϵt/\tau _0`$ that can be probed is about 3 to 4. In comparing with experiments and in making order of magnitude estimates, we will therefore only use the first correction term. Finally it is important to realize that after how long a time these expressions become really accurate, depends also on the initial conditions.

### B Dependence on initial conditions

In order to illustrate how accurate these expressions are in practice, we numerically integrate the real Ginzburg-Landau equation starting with an exponentially decaying initial condition with steepness $`\mu `$ : $`A(x,t=0)=ae^{\mu x}`$. The result is shown in Fig. 1 for various values of $`\mu `$. We see that initially the velocity falls off quickly and then approaches the asymptotic velocity from below, and that the asymptotic expression (8) in practice does yield a reasonable estimate of the time dependent velocity for values of the scaled time of order 3 and larger.

We also note that the theoretical analysis shows that for initial profiles falling off exponentially with steepness $`\mu <\mu ^{}`$, the asymptotic velocity lies above $`v^{}`$ and is given by $`v_{as}(\mu )=\mu +\frac{1}{\mu }`$. The dotted line in Fig. 1 shows an example of a case with $`\mu `$ slightly less than $`\mu ^{}`$, for which $`\stackrel{~}{v}_{as}=2.05`$. In this case, the time-dependent scaled velocity is approximately equal to 2 at times of order 3 to 4.

In the experiments, the typical protocol was to simultaneously increase the heat flux from an initial state at $`ϵ_i<0`$ (with a typical value of $`ϵ_i=0.015`$) to a supercritical value $`ϵ_f`$, and switching on the end heater. A rapid convergence to the asymptotic velocity due to special initial conditions with $`\lambda \lambda ^{}`$ would have required setting $`ϵ_iϵ_f`$ for all $`ϵ_f`$ and would probably also have required switching on the end heater before the value of $`ϵ`$ was changed, so that effectively at the end a convection pattern was prepared whose amplitude decayed exponentially into the bulk. Thus, it is theoretically possible to select special initial conditions so as to get a scaled velocity around 2 at some finite time, but the finetuning necessary to do so is so sensitive that we consider it very unlikely that experimental observations done over a range of values of $`ϵ`$ are due to initial condition effects. As stressed before, only new experiments can completely rule out this possibility, however.

## III Reexamination of the Rayleigh-Bénard fronts

In the long quasi-one-dimensional Rayleigh-Bénard cell of FS  the front propagation was initiated by simultaneously increasing the heat flux to a supercritical value $`ϵ>0`$ ($`ϵ=\frac{\mathrm{\Delta }T}{\mathrm{\Delta }T_c}1`$) and switching on a heater at the ends of the long cell. A vortex front is induced near the this heated end-wall and the propagation of this front into the unstable conductive state is then studied. Thin fins were attached to the long sides of the cell to avoid both induction of long rolls and pinning. Both because of the fact that the initial perturbation was caused by heating at the end and because the state before bringing the temperature difference beyond its critical value was unrelated to the final value of $`ϵ`$, there is every reason to believe that the experimental protocol did not create any special initial conditions that can not be considered sufficiently steep or localized.

As FS point out, since the front velocity $`v^{}`$ grows as $`\sqrt{ϵ}`$ while the growth rate in the bulk grows as $`ϵ`$, fronts can only be observed and in fact dominate the dynamics for small enough $`ϵ`$. In practice the pattern could be distinguished from the bulk noise up to a time $`t_{bg}=n\tau _0ϵ^1`$ where the numerical factor $`n`$ is of order 3–4. This determines some upper limit $`ϵ_0`$ on $`ϵ`$ for which the front can advance of the order of $`(1/m)`$th of the cell length $`l`$ before bulk growth takes over:

$$ϵ_0=\left(\frac{\stackrel{~}{v}\xi _0mn}{l}\right)^2.$$
(9)
Let us now estimate, using our analytical estimate (8), the relative importance of the correction term at the latest times of order $`t_{bg}=n\tau _0ϵ^1`$ at which measurements can be taken. Substitution gives

$$\stackrel{~}{v}(t_{bg})2\left(1\frac{3}{4n}\right).$$
(10)
Thus, for the latest time accessible in the experiment, we obtain with the empirical experimental value $`n=`$3–4 a velocity which is of order $`20\pm 5\%`$ below the asymptotic value. Although the asymptotic formula may not be accurate yet at such early times, the numerical results of Fig. 1 lead one to expect corrections of the same order of magnitude: these times also correspond to scaled times of order $`n`$ in the numerical simulation plots of Fig. 1, and as we have seen, over this time range, the velocity is also suppressed by about 15% relative to the asymptotic value. Hence, by any reasonable estimate, the slow relaxation can not be negligible in the experiments!

As stated in the introduction, FS quoted a measured velocity $`\stackrel{~}{v}=2.01\pm 0.02`$ for their experiments, contrary to our expectation that their scaled velocity should have been in the neighborhood of 1.7. According to the theoretical results, one should expect to get to within 1% of the asymptotic velocity only around a dimensionless time of order 100!

In our view, the most plausible explanation for the origin of this discrepancy is that the value $`\xi _0`$ was actually larger in the experiments than the value used in the analysis. The value of $`\tau _0`$ was experimentally confirmed to be very close to the theoretical value . For $`\xi _0`$, however, the theoretical value was used without independent experimental check . Because of the special design of the cell with side-fins to create one-dimensional patterns, a different value of $`\xi _0`$ might not be unexpected. In fact, an indication that $`\xi _0`$ in the experiments was larger than the theoretical value used in the analysis, comes from the observation of FS that the value of the wavelength at onset was 13% larger than the theoretical value. This might indicate that all lengths in the experiments are a factor 1.13 larger than the theoretical values, and this is precisely the factor needed to reconcile the front data with the theoretical expectations! One should keep in mind, though, that $`\xi _0`$ is determined by the the curvature of $`ϵ`$ versus $`k`$ tongue around the critical wavenumber, and that it is not guaranteed that both are changed by the same factor; only independent measurements can fully settle this issue.

We now show that a reanalysis of the data of FS actually gives quite convincing evidence for slow convergence effects in the experimental fronts.

FS measured the velocity by comparing the front with itself at various time intervals, by appropriately shifting the traces back. This yielded a set of points in the $`\mathrm{\Delta }x`$, $`\mathrm{\Delta }t`$ plane which appeared to lie on a straight line. However a possible relaxation of the velocity was masked since points from early and later times will approximately fall in the same place in the plane, and because the front shape also has an asymptotic $`1/t`$ relaxation. We therefore have tried to reanalyze the raw data of figure 2 from the FS paper; this space-time plot of a propagating from is reproduced in the top panel of our Fig. 2. We define the position of the front as the point where the interpolation of the maxima of the profile equals some fraction of its maximum in the bulk (we chose 0.4). Our data for $`\stackrel{~}{v}(t)`$ obtained this way are shown in the lower panel of Fig. 2. Whereas the local velocity initially slightly decreases, an increase for dimensionless times larger than 0.5 is evident. We stress that this qualitative behavior is independent of the choice of the parameters. In order to compare quantitatively to the predictions for the relaxation, we have used the value for $`\tau _0`$ given by FS but increased the value of $`\xi _0`$ by 13% on the basis of the argument given above. Clearly, with this choice, the data are certainly consistent with the analytical as well as numerical estimates of the velocity relaxation — in fact, in a way the data are the first experimental indications for the universal power law relaxation of pulled fronts .

## IV Relaxation of wavelength

FS also studied the problem of the selection of the wavelength $`\lambda `$. As we discussed before, the actual value of the wavelength $`\lambda _c`$ of the patterns is at criticality about 13% off from the theoretical value; however, we are not interested here in the absolute value, but in the relative variation of $`\lambda _c/\lambda `$.

The difficulty of comparing theory and experiment on the variation of the wavelength is that the only theorically sharply defined quantity is the wavelength sufficiently far behind the front, $`\lambda _{as}`$, and that one has to go beyond the lowest order Ginzburg-Landau treatment to be able to study the pattern wavelength left behind. E.g., if we use a Swift-Hohenberg equation for a system with critical wavenumber $`k_c`$ and bare correlation length $`\xi _0`$,

$$_tu=\frac{(\xi _0k_c)^2}{4}\left(1+\frac{1}{k_c^2}\frac{^2}{x^2}\right)^2u+ϵuu^3,$$
(11)
then a node counting argument  yields for the asymptotic wavelength $`\lambda _{as}`$ far behind the front :

$`{\displaystyle \frac{\lambda _c}{\lambda _{as}}}`$ $`=`$ $`{\displaystyle \frac{3\left(3+\sqrt{1+24ϵ/(k_c^2\xi _0^2)}\right)^{3/2}}{8\left(2+\sqrt{1+24ϵ/(k_c^2\xi _0^2)}\right)}}`$ (12)
$``$ $`1+ϵ/(2k_c^2\xi _0^2)(ϵ1).`$ (13)
In the Rayleigh-Bénard experiments, $`k_c2.75/d`$, where $`d`$ is the cell height; the theoretical value is $`\xi _0=0.385d`$, so if our conjecture that the value is some 15% larger is correct, we get $`\xi _00.44d`$. This then gives

$$\frac{\lambda _c}{\lambda _{as}}1+0.34ϵ.$$
(14)
As we stressed already above $`\lambda _{as}`$ is the wavelength far behind the front; for a propagating pulled front, there is another important quantity which one can calculate analytically, the local wavelength $`\lambda ^{}`$ measured in the leading edge of the front. For the Swift-Hohenberg equation, one gets for this quantity

$`{\displaystyle \frac{\lambda _c}{\lambda ^{}}}`$ $`=`$ $`\sqrt{1+{\displaystyle \frac{1}{4}}\left(\sqrt{1+24ϵ/(k_c^2\xi _0^2)}1\right)}`$ (15)
$``$ $`1+{\displaystyle \frac{3ϵ}{2k_c^2\xi _0^2}}(ϵ1)`$ (16)
$``$ $`1+ϵ,`$ (17)
where in the last line we have used the experimental values.

Let us now discuss the experimental findings in the light of these results. In Fig. 3 we show both the experimental values of $`\lambda _c/\lambda _{exp}`$ as determined experimentally by FS, and results for the Swift-Hohenberg equation with the value (14) relevant for the experiments.
(i) For small $`ϵ`$, the experimental values are roughly linear in $`ϵ`$, but when a fit is made over the whole range of $`ϵ`$ values studied, a square root behavior, as proposed by FS, would probably be better.
(ii) The experimental values for the wavelength ratio deviate about a factor of 3 from the theoretically expected value for the ratio far behind the front, $`\lambda _c/\lambda _{as}`$.
(iii) Just like the front speed converges very slowly to the asymptotic value, so does the local value of the wavelength behind the front . The relaxation of the velocity is plotted with a dashed-dotted line in Fig. 3 for times about $`6\tau _0/\sqrt{ϵ}`$, which is the time it takes for a front to propagate to close to the center of a system about as large as the experimental cell. Clearly, the wavelength ratio due to slow relaxation lies below the asymptotic value, and hence further away from the experimental data for small $`ϵ`$. Also, numerical results for the wavelength of the fourth “roll” measured at the same time (indicated by the dotted line, lie below the asymptotic curve.
(iv) FS have measured the wavelength for very small values of $`ϵ`$, down to $`\mathrm{4\hspace{0.33em}10}^4`$. However, already for values as small as 0.01, the coherence length $`\xi =\xi _0/\sqrt{ϵ}`$ is about $`2\lambda _c`$ in the experiments. The total front width is several times this number, and for even smaller values of $`ϵ`$ the front width is even smaller. Since the total experimental cell was about $`12\lambda _c`$ long in the experiments, it is clear that over much of the small $`ϵ`$ range, one would not expect to see well-developed fronts in the experiments. In other words, it is very unlikely that over the experimental range of $`ϵ`$-values, one has a chance to measure $`\lambda _{as}`$ of a well-developed pattern behind a front.
(v) The grey line in the plot shows our analytical result for the wavelength ratio in the leading edge of the front. Clearly, this line follows the data for small $`ϵ`$ quite well. In view also of point (iv) above that it will be hard to obtain well-developed fronts for small $`ϵ`$, we propose as a tentative explanation of the data that in the small $`ϵ`$ range, one actually measures the emergent roll pattern associated with the leading edge of a front. Of course, only new experiments can decide on the validity of this suggestion.
(vi) We mention that the variation of the wavelength ratio with $`ϵ`$ depends also on the third order derivative term in the expansion of the dispersion relation around the critical wavenumber. This term is not modeled correctly in the Swift-Hohenberg equation, but may have to be taken into account in a full comparison of theory and experiment.
(vii) We finally mention that in the experiments there was an up-down asymmetry in the rolls. We have investigated whether this could be a source of the discrepancy between the asymptotic wavelength ratio and the observed one, by studying a Swift-Hohenberg equation with a symmetry-breaking quadratic term. However, with this term, the wavelength ratio apears to decrease away from the experimental values.

## V Conclusion

It was recently discovered that quite generally pulled fronts relax very slowly to their asymptotic velocity. Comparison of the experimental data for the velocity with numerical simulations and analytical estimates give, in our view, evidence that these experiments provide clear signs of the presence of such slow relaxation effects, although the time scales that can be probed experimentally are too short to test the general power law relaxation. Theoretically the only other viable option to reconcile theory with the interpretation originally proposed  is that somehow special initial conditions created an initial convection profile with precisely the right spatial decay into the bulk. As we discussed in section III, in our view this is an unlikely interpretation, but only new experiments can settle this issue completely.

While measurements of the wavelength of the pattern generated by a front are even more difficult to interpret than those of the velocity, our analysis indicates that in the small-$`ϵ`$ regimes a well-developed front does not fit into the experimental cell, and that as a result one probes the local wavelength in the leading edge of the front rather than the well-developed asymptotic wavelength behind it. The analytical estimates are consistent with this suggestion.

We hope that this work will trigger new experimental activity to investigate these issues — experiments along these lines hold the promise of being the first ones to see the universal power law relaxation of pulled fronts.

## VI Acknowledgement

We are grateful to Jay Fineberg and Victor Steinberg for correspondence about their work and about the issue discussed in this paper. WvS would also like to thank Guenter Ahlers for urging him to analyze to what extent slow convergence plays a role in the Rayleigh-Bénard experiments. JK is grateful to the ‘Instituut-Lorentz’ for hospitality.


text:
# Evidence of local magnetic order in hcp iron from Raman mode splitting

## Abstract

Experimental measurements of Raman spectra for hcp iron at high pressure show two modes over a considerable pressure range in contrast to the prediction of one doubly degenerate mode for the hcp lattice. We use density functional theory to investigate the influence of magnetic order on the Raman active modes of hcp iron. We find an antiferromagnetic state that lifts the degeneracy of the transverse optical mode, and yields stable antiferromagnetic moments up to approximately 60 GPa (55 Bohr3). The resulting frequencies of the two transverse optical modes are in good agreement with the experimental Raman shifts, lending support to the existence of local antiferromagnetic order in hcp iron.

The presence of magnetism strongly influences material properties and phase relations in the $`3d`$ transition metals. This is particularily true for iron where the competition between magnetic and non-magnetic contributions to the internal energy, differences in the vibrational and magnetic entropy, and differences in volumes all contribute to phase stability moroni . Ferromagnetism stabilizes the body centered cubic structure over the close packed phases at ambient condition bagno , and an understanding of the spin-density wave ground state in its high temperature polymorph, cubic close packed (fcc), has important implications on its physical properties uhl . Magnetism in the high pressure polymorph, hexagonally close packed (hcp), is less well characterized and has long been thought to be absent, but its proposed presence steinle  may play an important role in our understanding of material properties at pressure, with applications to the study of planetary interiors stixrude  and impact phenomena.

The hcp lattice has two atoms per unit cell giving rise to optical phonon modes through Brillouin zone folding. Advances in experimental techniques have made it possible to measure the Raman active transverse optical ($`TO`$) mode in iron through optical spectroscopy merkel ; olijnyk . In contrast to the fundamental prediction that this mode ($`E_{2g}`$) be doubly degenerate, one set of experiments at high pressure show two peaks in the Raman spectrum up to pressures of 40 GPa merkel . While amplitude and sharpness of the two peaks differ significantly, this observation nonetheless suggests that the symmetry of hcp iron is lower than the atomic arrangement; spin-phonon interactions provide a symmetry-breaking mechanism. Based on first principles theory the lowest energy state found to date is an anti-ferromagnetic (afm) structure with orthorhombic symmetry and four atoms in the unit cell (afmII) steinle . The afmII structure corresponds to a spin wave with wave vector $`q=(0,1/2,0)`$, at the $`M`$ point on the Brillouin zone boundary (Fig. 1). This structure has been predicted to be stable for hcp iron up to pressures of almost 60 GPa steinle . The magnetic structure results in two $`TO`$ zone center modes ($`A_{2g}`$) both of which are Raman active.

Here we investigate the influence of afmII ordering on the zone center $`TO`$ phonon mode over the compression range where finite moments arepredicted ($``$60 Bohr3/atom). We calculate the zone center Raman frequencies by the frozen phonon approximation: energy changes are evaluated in response to small displacements along the phonon eigenvector, the second order term yielding the frequency. We base our analysis on total energies obtained with the full-potential linearized-augmented plane-wave method (LAPW) singh  using a generalized gradient approximation pbe  to the exchange correlation potential. We treat $`3s`$, $`3p`$, $`3d`$, $`4s`$, and $`4p`$ states as valence electrons for all volumes and use $`R_{MT}`$=2.0 Bohr for the muffin tin radii, $`R_{MT}K_{max}=9.0`$, a 12$`\times `$6$`\times `$12 special k-point mesh, and a temperature broadening of 5 mRy. With this set of computational parameters we have previously established convergence of relative energies to within $`0.1`$ mRy/atom and magnetic moments to better than $`0.05`$ $`\mu _B`$/atom steinle . The $`TO`$ modes are characterized by displacements of the close packed planes with respect to one another. The eigenvectors for the the two $`A_{2g}`$ modes in the afmII structure are along the orthorhombic $`a`$\- and $`b`$-axes (Fig. 1). For comparison we also compute the $`TO`$ mode frequency for the non-magnetic case. To efficiently evaluate the energetics of the system we fixed the axial ratio in our calculations to $`c/a=1.6`$, close to the experimentally jephcoat ; mao90  and theoretically steinle  estimated equilibrium value.

Inspection of the afmII structure reveals that fundamentally different spin interactions are involved in the two $`TO`$ modes (Fig. 1). For displacements along the orthorhombic $`a`$-axis ($`TO_a`$) atoms approach nearest neighbors with unlike spin. For displacements along the orthorhombic $`b`$-axis ($`TO_b`$) atoms alternately move towards a nearest neighbor with like spin and towards a pair of nearest neighbors of opposite spin. The resulting energy - displacement relations reflect the antiferromagnetic spin interactions (Fig. 2): for unlike spins approaching the energy is reduced with respect to the non-magnetic case, and for same spin nearest neighbor interactions (along positive $`b`$) there is an additional repulsion. The calculated equilibrium magnetic moments are consistent with these findings (Fig. 2): They increase for displacements along $`a`$, as they do along negative $`b`$; for positive $`b`$ they decrease considerably as the like spin nearest neighbors approach.

We find that the $`TO_a`$ mode frequency agrees well with that of the lower frequency, higher amplitude peak found in the Raman experiments, and that the $`TO_b`$ mode frequency corresponds to the experimentally observed satellite peak at higher frequency (Fig. 3). The magnitude of the predicted $`TO`$ mode splitting decreases as the afm moment is reduced by compression, in excellent agreement with the observation in the Raman experiments (Fig. 3). The magnitude of splitting is related to the directionality of the spin wave considered here (with wave vector $`M`$), but other spin waves with wavevectors along high symmetry directions in the base of the hexagonal Brillouin zone will generally result in splitting of the $`TO`$ mode. The systematic offset of the calculated frequencies by approximately 20 cm-1 ($`<`$ 10%) is typical for a comparison of computed and measured phonon frequencies znolijnyk ; althoff .

Spin-phonon interactions have been found to have a strong effect on Raman scattering in a number of systems, including cupric oxide chen95  and the copper-ruthenium oxide RuSr2GdCu2O8 fainstein . The general character of the effect of spin-phonon interaction on the Raman spectra in these materials is consistent with the observations for hcp iron: broad, low amplitude, satellite peaks appear as the sample is cooled below the Curie or Neel temperature, upon further cooling such peaks generally sharpen and gain amplitude as the degree of magnetic ordering increases. Observation of an increase in amplitude of the satellite peak in hcp iron upon cooling would lend support to the predicted magnetic structure.

Other experimental investigations of possible magnetic states in hcp iron have been inconclusive. Magnetism in hcp iron exists at over-expanded volumes in epitaxially grown multilayers on a ruthenium substrate maurer  although the nature of the spin ordering is still controversial lager ; knab . Nuclear X-ray absorption experiments rueff  show significant loss of moment across the phase transition from bcc to hcp iron, but can not unambiguously be interpreted to show no moments in the high pressure polymorph: the change in absorption spectra is due to changes in the density of states as well as to spin related satellites. Diamond anvil cell in-situ Mössbauer measurements on hcp iron have not shown any ordered magnetic ground state in its stability field nicol ; cort ; taylor1 ; taylor2  but it was recognized that magnetism cannot unambigously be ruled out based on the data nicol ; taylor1 . The Mössbauer data suggest that correlation times of spin fluctuations must be relatively short, or that competing contributions to the effective field cancel each other as is the case in CaRuO3 perovskite gibb ; felner . The in-situ observation of super-paramagnetism in hcp-Fe taylor1  may support the notion that the effective field is approximately cancelled, and suggests, following the procedure of Felner et al. felner , a Mossbauer experiment on an appropriately doped hcp-Fe sample (e.g. with Ru). The inherent frustration of the triangular lattice might lead one to suspect more complex spin arrangments than those considered here, such as incomensurate spin waves, as in the case of fcc iron tsunoda ; uhl , non-colinear structures, spin glass, or a combination of these cohen01 .

Acoustic modes can also be influenced by spin-phonon interactions as is evidenced by anomalous phonon dispersion in fcc iron near the zone center zaretsky . Indeed we find that the afmII ordering also influences the compressional properties of hcp iron steinle : the bulk modulus at pressure is reduced, bringing it in closer agreement with experimental data jephcoat ; mao90  than non spin-polarized calculations. Better agreement of aggregate elastic properties between afmII and experimental data lends further, independent, support to the presence of magnetism in hcp iron.

The possible existence of magnetic correlations in hcp iron may be important in our understanding of the recent observation of superconductivity in hcp iron under pressure shimizu . Long thought to be antithetical, superconductivity and magnetism have recently been observed simultanously in a number of systems saxena ; pfleiderer ; aoki . Of particular interest is the observation for ZrZn2 where magnetism and superconductivity appear to be directly coupled pfleiderer  as is evidenced by the loss of superconductivity and magnetism at the same pressure, leading to the speculation that hcp iron might behave similarly saxena01 .

In conclusion, we report first-principles results on the zone center $`TO`$ phonon frequency of hcp iron for non-magnetic and afm structures at ambient and high pressure from a frozen phonon calculation. We find good agreement with experimental results for the absolute frequencies, and spin-phonon interactions provide a quantitative explanation for the split of the $`TO`$ mode in hcp iron which is observed up to a pressure of 40 GPa. In combination with a considerably better agreement of elastic properties between experiment and the orthorhombic afm structure, as compared to non-magnetic calculations, this suggests that magnetism may play an important role in the physical behavior of hcp iron at high pressure.

###### Acknowledgements.

We greatly appreciate helpful discussion with A. Goncharov, R. Hemley, I. Mazin, S. Merkel, H. Olijnyk, and V. Struzhkin, and thank H. Krakauer and D. Singh for use of their LAPW code. This work was supported by the National Science Foundation under grants EAR-9980553 (LS) and EAR-998002 (REC), and by DOE ASCI/ASAP subcontract B341492 to Caltech DOE W-7405-ENG-48 (REC). Computations were performed on the Cray SV1 at the Geophysical Laboratory, supported by NSF grant EAR-9975753 and by the W. M. Keck Foundation.