text:
# Probing complex RNA structures by mechanical force

RNA secondary structures of increasing complexity are probed combining single molecule stretching experiments and stochastic unfolding/refolding simulations. We find that force-induced unfolding pathways cannot usually be interpretated by solely invoking successive openings of native helices. Indeed, typical force-extension responses of complex RNA molecules are largely shaped by stretching-induced, long-lived intermediates including non-native helices. This is first shown for a set of generic structural motifs found in larger RNA structures, and then for Escherichia coli’s 1540-base long 16S ribosomal RNA, which exhibits a surprisingly well-structured and reproducible unfolding pathway under mechanical stretching. Using out-of-equilibrium stochastic simulations, we demonstrate that these experimental results reflect the slow relaxation of RNA structural rearrangements. Hence, micromanipulations of single RNA molecules probe both their native structures and long-lived intermediates, so-called “kinetic traps”, thereby capturing –at the single molecular level– the hallmark of RNA folding/unfolding dynamics.

Keywords: RNA folding/unfolding; Single molecule experiments; Stochastic simulations; Non-native helices and kinetic traps; 16S ribosomal RNA.

Introduction

Recent developments of micromechanical experiments on single biomolecules have provided structural insights into alternative structures of DNAchatenay ; bustamante ; leger1 ; allemand  and mechanical properties of proteinsprot1 ; prot2 ; prot4 . In principle, such techniques could also provide new tools to probe RNA structures which remain by and large refractory to many crystallization schemes. However, this prospect requires one to relate mechanically induced unfolding pathways to RNA structural features. Although it could be done successfully for small RNA structures by solely invoking successive openings of native helicesliphardt1 , probing more complex RNA structures by mechanical force is expected to involve non-native structural rearrangements of the initial secondary structure upon stretchingmontanari ; bundschuh ; lubensky . Local rearrangements, such as the formation of simple stem-loops, occur quite fast ($`<`$ 1 ms) under low pulling force (or in the absence of force) and the number of possible hairpins (with small loop) is proportional to the length of stretched region of the RNA molecule. Thus, alternative hairpins, not present on the initial structure, should inevitably form under partial stretching of long RNA molecules (e.g., $`>`$ 1000 nucleotides). Conversely, more global rearrangements, which involve the coordonated removal and formation of different sets of helices, might occur much more slowly (e.g., after a few minutes)williamson1 ; woodson . Hence, under typical pulling rates (i.e., full extension within a few seconds), most stretching experiments likely occur under out-of-equilibrium conditions and should exhibit unfolding/refolding hysteresis curves.

To study the full potential and limitations of these micromechanical techniques so as to probe complex RNA structures, we have combined single RNA molecule stretching experiments and out-of-equilibrium stochastic simulations. Three small artificial structures, M1, M2 and M3 (Fig 1), representing prototypes for the main structural modules of larger RNA secondary structures, were first designed and studied in details. The mechanical response of E. coli’s 1540-base long 16S ribosomal RNA was then studied using the same experimental setup and a somewhat simplified numerical approach. The generic structural motifs M1, M2 and M3, correspond to three different arrangements of two 15 base pair long helices consisting almost exclusively of either GC or AU base pairs, Fig 1. M1 corresponds to two adjacent stem-loops with respect to the external single strand joining the molecule ends. By contrast, M2 and M3 present the same nested organisation with either the strong (GC) helix or the weaker (AU) helix connected to the external single strand.

The 5’ and 3’ ends of either these small RNA motifs or E. coli’s 16S rRNA were hybridized to two pUC19 dsDNA extensions labelled, respectively, with biotin and digoxygenin (see Materials and Methods). The force-extension experiments were then done by grafting the ends of these extended molecular constructs between the antidigoxygenin coated glass surface of a capillary and a micrometer size silica bead coated with strepavidin. The capillary was moved by a piezo-electric stage (50 to 300 nm/s) and the resulting force exerted by the molecule on to the bead was measured with an optical tweezer.

Results

Single molecule stretching experiments of small RNA motifs

When structural motifs M1, M2 or M3 are inserted in the molecular construct, one or two force drops occur on the force-extension curve, Fig 2. A histogram of the rupture force and a set of unfolding and refolding force-extension curves are shown for each motif on Fig 2A and Fig 2B. For each set, variations between force-extension curves correspond to stochastic fluctuations between either successive stretchings on the same molecule or different experiments on equivalent molecules.

A comparison of the different rupture force histograms and the corresponding unfolding curves (black on Fig 2B) shows that M1 and M2 present very similar unfolding responses with two sequential drops or inflexion regions around 11 $`\pm `$ 3 pN and 22 $`\pm `$ 3 pN, whereas M3 presents a single and larger force drop at about 22 $`\pm `$ 3 pN. These values are in very good agreement with ref.gaub  although these latter experiments concern the opening of DNA hairpins. For M1 and M2, these results can be simply attributed to the first opening of the weak (AU) helix followed by the stronger (GC) helix at a higher applied force. Indeed, the applied tension being uniformly distributed along the external single strand joining the molecules ends, the weaker (AU) helix is expected to break first on M2, while it should certainly do so by construction on M1. Besides, by calibrating the stiffness of the optical trap, both force drops on these curves can be converted into a distance released by the molecule, taking into account the angular inclination of the setup (30 to 40). This corresponds to the expected 20 nm in both cases. Substracting the net free-energy contribution stored in the stretched single strandcocco , we find in term of pairing energy, around 1.7 kT/bp for AU and 3 kT/bp for GC, in good agreement with known parametersturner . By contrast for M3, the strong (GC) helix shields the weaker (AU) stem from the applied force and no significant unzipping is observed until the whole molecule suddenly unfolds at the critical force to break GC stacking base pairs.

The refolding curves for M1 (blue) and M2 (red) show most often a small hysteresis below the force drop associated with the strong (GC) stem’s opening. For M2, a second small refolding hysteresis occurs also usually below the force drop associated with the weaker (AU) stem’s opening. By contrast, a much stronger hysteresis is systematically observed for M3 (magenta), even at the lowest loading rate achieved, 3 pN/s. Moreover, in this case, the refolding event around 10 pN does not usually fold back onto the initial stretching curve. This suggests that the stretching of M3 involves long-lived intermediate structures including non-native helices (see Fig 3 and next section on stochastic simulations of small RNA motifs). Still, all three molecules eventually fold back in their initial native structure after a few seconds, as shown by the reproducibility of force extension curves in successive pulls on the same molecule.

Stochastic unfolding/refolding simulations of small RNA motifs

We have performed stochastic simulations of these out-of-equilibrium unfolding/refolding experiments for the short M1, M2 and M3 structural motifs. The heart of the numerical method, following the approach detailed inisambert , consists in simulating the stochastic unfolding and refolding of helices not only present on the initial RNA structure but also for all other helices which can possibly pair on the RNA sequence of interest (see also RNA Kinefold server at http://kinefold.u-strasbg.fr ). Common pseudoknots (i.e., helices interior to loops) are also allowed following the structural modeling approach proposed inisambert . In addition, the region of the RNA structure under direct mechanical tension (corresponding to the “on-net” backbone inisambert ) is modelled as an inextensible wormlike chain with a 1.5 nm persistence lengthbustamante  and 0.7 nm/base contour lengthcontour . Stretching is induced by a slowly varying rigid constraint on the end-to-end distance of the RNA-dsDNA-tweezer construct (rate $`\pm `$ 300 nm/s in $`\pm `$ 2 nm steps). It is also important to take into account the acquisition rate (300 Hz) and to model the statics and dynamics of the optical tweezer trap, although time scale separation allows to consider that the micromechanical setup responds to a slow time average of the fast RNA dynamics, which corresponds to stochastic closing and opening of single helicesisambert . To avoid overfitting with non essential parameters, we have simply modelled the trapped bead and the two dsDNA extensions of the construct as an ideal spring with a slow viscous relaxation time (1 ms) and an effective stiffness (typically 0.1 pN/nm) fitted on the individual force extension curves.

The simulated force-extension responses for the M1, M2 and M3 motifs (Fig 2C) are in good agreement with the unfolding and refolding experimental results (Fig 2B). In particular, they allow for the identification of likely intermediate structures involved in the refolding hysteresis, which primarily correspond to the formation of two non-native helices originating from each strand of the strong (GC) helix, Fig 3. The transition from these alternative helices back to the strong initial (GC) stem is facilitated under high external force, hence the small hysteresis for M1 and M2. By contrast for M3, the transition can only occur at a lower force after the weak (AU) stem has refolded and is, therefore, slower, as observed experimentally. Note, the small experimental differences in the hysteresis responses of M1 and M2 are well reproduced on their simulated force-extension curves, suggesting that elementary unfolding/refolding events are reliably captured by these stochastic simulations devised to probe RNA molecular dynamics on second to minute time scales. It should be emphasized that such long time scale simulations could not be achieved for other molecules, such as proteins, for which elementary unfolding/refolding transitions are not so easily defined and also much more frequent. For instance, the best molecular dynamic simulations of proteins are currently limited to around 100 nsfersht .

The force-induced unfolding of these three generic structural motifs M1, M2 and M3, reveals the potential and limitations of single molecule experiments to probe the main folded features of more complex RNA structures. The comparison between M2 and M3’s force-extension responses illustrates that the order of helix stability along a single secondary structure branch can be readily identified, while the bifurcation arrangements of helices or the presence of multibranched loops are not so easily distinguished from single branches with increasing helix stability (as in M1 versus M2). Moreover, the formation of non-native rearrangements under stretching likely affects the force-extension responses of most RNA structures (as for M3).

In this context, combining experimental and numerical approaches to study RNA mechanical unfolding pathways seems promising insofar as transient structural rearrangements (under stretching) are difficult to probe with traditional chemical or enzymatic techniques.

On the other hand, for long RNA molecules (e.g., $`>`$1000 bases), it has been arguedmontanari ; bundschuh ; lubensky  that such structural rearrangements under stretching should ultimately smooth out the observed characteristics completely by continuous adjustments to the applied constraint, assuming that quasi-equilibrium stretching is achieved.

To investigate this issue and test whether large structures of biologically relevant RNA molecules are also amenable to convergent studies in both single molecule experiments and stochastic simulations, we decided to study the mechanical unfolding of Escherichia coli’s 1540-base long 16S ribosomal RNA.

Single molecule stretching experiments of E. coli 16S rRNA

The force-induced stretching of E. coli 16S ribosomal RNA was studied using a similar molecular construct and micromechanical setup as for the stretching of the small structural modules M1, M2 and M3 (see Materials and Methods). No ribosomal proteins which associate to 16S rRNA to form the 30S subunitramakrishnan  of the ribosomenoller1  were included for these stretching experiments. As the piezo stage is displaced, the force begins to rise due to the elastic response of the DNA handles. The results on Figs 4-5 show a well-structured and reproducible unfolding pathway under mechanical stretching, in about 50% of the tested constructs for which more than two unfolding/refolding rounds could be performed before molecular breakage. In these cases (total 44 stretching curves), a $``$ 1 $`\mu `$m-long quasi plateau is observed around 11-15 pN, with force fluctuation amplitude of about 20%. This signal is the signature of 16S rRNA unfolding by mechanical force. Other stretching curves exhibit somewhat more erratic behaviors, presumably due to non-specific interactions of the construct with the glass surface of the capillary (data not shown) Extension beyond the unfolding plateau corresponds to the combined elastic response of the dsDNA handles and the opened ssRNA molecule. Most refolding curves exhibit strong hysteresis depending on stage velocity (50-300 nm/s).

We have quantitatively evaluated the statistical reproducibility of unfolding curves between successive stretchings of the same 16S molecule (Fig 4) and between independent unfolding curves from different 16S molecules (Fig 5) (i.e., different sample preparations and different micromechanical experiments). The analysis is based on the calculation of Spearman nonparametric correlation coefficient $`r_s`$ recipes : $`r_s=\mathrm{\Sigma }_i^N(R_i\overline{R})(S_i\overline{S})/(\mathrm{\Sigma }_i^N(R_i\overline{R})^2\mathrm{\Sigma }_i^N(S_i\overline{S})^2)^{1/2}`$, where $`R_i`$, $`\overline{R}`$ and $`S_i`$, $`\overline{S}`$ are chosen as the rank-ordered deviations and averages from the median line fits of the unfolding plateaux (Fig 6). Such relative Spearman correlation coefficient is much more sensitive to the specific unfolding signals, as compared to the absolute Spearman correlation of the actual unfolding curves which yields higher, yet less discriminating correlation coefficients (see Fig 6 for details).

In the context of comparing 16S unfolding curves, we found that this relative Spearman correlation coefficients correspond to good correlations above 50% and excellent ones above 70%, while $`|r_s|<`$15% reflects little or no correlation between unfolding pathways irrespective of the overall inclination of their unfolding plateaux. Stochastic reproducibility between successive stretchings of the same molecule is remarkably high ($`r_s`$=75$`\pm `$4.6%, Fig 4) and still quite good between unfolding curves of different 16S RNA ($`r_s`$=53$`\pm `$9.1%, Fig 5) despites inherent variations between different sample preparations and different micromechanical experiments. For instance, correlation between independent green and red unfolding curves on Fig 6 is: $`r_s`$=60.5%, while the correlation distributions of all independent curves of Fig 5 (colors) with their average unfolding response (black) is even higher: $`r_s`$=70$`\pm `$4.7%. See Table 1 for further correlation data and quantitative comparison with stochastic simulations.

Before discussing the stochastic unfolding/refolding simulations of 16S rRNA, we want to emphasize that force fluctuations from the plateau median line cannot be attributed to dehybridization of the DNA handles for the following reasons: i) we never recorded such signals on simple pUC19 dimer without RNA insert (see Fig 1 inset). ii) experiments done by other groups with the same infrared laser power have preserved nucleic acids’ integritybustamante ; svoboda . iii) recalling that we pull on opposite DNA strands, the average force magnitude at which the plateau appears is too low to originate from DNA denaturationchatenay ; bustamante ; leger1 . Moreover, no torque is applied on the molecules using our optical tweezerallemand . In addition, experimental force fluctuations cannot correspond to the unzipping of a long structureless double stranded RNA molecules; an analysis based on G+C contents as in ref.bock  does not account for the experimental signal, nor does a thermal equilibrium energy calculations.

Stochastic unfolding/refolding simulations of 16S ribosomal RNA

The unrestricted stochastic simulations discussed above to model the mechanical unfolding and refolding of small structural motifs are numerically unpractical in the case of much larger RNA structures like those of ribosomal RNAs. Hence, we have made the following three additional assumptions to study the force-induced stretching of E. coli’s 1540-base long 16S ribosomal RNA: i) the initial structure before mechanical stretching is assumed known from independent sources; ii) Unfolding and (re)folding dynamics is restricted to the formed helices under direct mechanical tension and to all potential helices that would be under direct tension once formed. Hence, large scale structural rearrangements can only originate and propagate from helices directly coupled to the applied mechanical tension, as expected under strong stretching conditions; iii) For each intermediate structure along the unfolding pathway, the actual base pair extent of each helix under direct tension is not globally optimized to best fit the end-to-end molecular extension imposed by the mechanical setup (this would become exponentially difficult in the number of such helices). Instead, a local heuristics extending the most stable base pair stacks and shrinking the weakest helix ends is used iteratively to minimize free energy. This approach, which yields a linear optimization in the number of helices under direct tension, is usually very good as long as there are few mutually incompatible helices competing for the same bases, a typical situation under strong mechanical stretching. Overall, we found that these restricted stochastic simulations give virtually identical results for the small M1, M2 and M3 motifs (results not shown).

Adopting this heuristic numerical approach for the bare E. coli

16S rRNA, we simulated the force-induced unfolding pathway starting either from the known native secondary structure inside the ribosomeramakrishnan ; noller1  or from a low free-energy structure predicted by mfold, referred at, hereafter, as the “control structure” (http://bioinfo.math.rpi.edu/$``$mfold/). The comparison between the force-extension responses of these two structures was primarily intended to probe the stochastic simulation’s sensitivity to the initial structure. In both cases, the role of helices not initially formed on the starting structure was studied, by allowing a variable number of helices to form and break during different stochastic simulations. The results in Fig 7 and Table 1 show that a reasonable agreement exists between the experimental measurements (black and brown curves) and the simulated force-extension responses starting from the known native structure, when all helices longer than 3bp and more stable than 10kT (i.e., 6 kcal/mol) are included a priori in the simulations (red curves). This demonstrates that some of these 18,000 different non-native helices play a significant structural role along the unfolding pathway. More quantitatively, cross-correlations amonsgt 4 independent simulated stretching curves (i.e., about 3 weeks of CPU on a 1.2GHz PC) reveal a good stochastic reproducibility in these simulations starting from the native structure and including about 18,000 possible helices: $`r_s`$=57$`\pm `$ 3.8%. This is comparable to observed variations between experimental unfolding curves (see Table 1). Then comparing these individual simulations with the experiment average curve (black curve on Fig 5 from 4 independent experiments), we obtain a significative correlation coefficient: $`r_s`$=47$`\pm `$8.8%, while correlating the experiment average (black) directly to the simulation average clearly reflects common features between the experimental response and the simulated unfolding pathway curve starting from the native structure: $`r_s`$=55%. In addition, restricting simulations to the 6,500 possible helices longer than 3bp and stronger than 15kT (magenta curves) or including even fewer helices (blue and green curves), produces marked differences on the simulated stretching curves (Table 1 shows lower averages and larger standard deviations for the simulation stochastic reproducibility and for the correlations with experimental response). By contrast, equivalent stochastic simulations starting from the control structure (Fig 7) present clearly distinct results from experimental observations, i.e., $`\left|r_s\right|<`$15% (even for a large number of possible helices included in the simulations). The fact that the stochastic reproducibility of these control simulations happens to decrease with the number of possible helices taken into account (Table 1) reflects the concomitant decrease of specific unfolding signal relative to the median line fit of the plateau (Fig 7). A decreasing signal over noise ratio naturally leads to a lower reproducibility of the simulated curves. The same trend is also visible between blue and magenta curves for simulations starting from the native structure.

Fig 8 compares more closely a simulated force-extension response of the known native structure (red) and an experimental stretching curve (black). Again, both simple visual comparison and calculation of their correlation coefficient as above (here $`r_s`$=61%) strongly suggest that the experimentally probed structure shares, indeed, more structural features with the actual native structure than with the control structure (Fig 7), in spite of the absence of ribosomal proteins in these single molecule stretching experiments of E. coli’s 16S rRNA. Analysing the unfolding pathway during the simulated force-extension response reveals that the main predicted unfolding events (corresponding to abrupt force drops on the red curve) are either related to the cooperative opening of several native helices (as in the unfolding of M3) or to the simultaneous rearrangements of mainly non-native helices leading to a stepwise increase of the predicted extension of the molecule along the direction of pulling (violet curve). This is illustrated with 12 successive snapshots of intermediate structures along the stretching-induced unfolding pathway. In particular, the 3’ major (III) and 3’ minor (IV) domains are shown to break and partially rearrange at the start of the stretching plateau (intermediates 2 to 4) while the 5’ domain (I), partially unfolded between intermediates 4 and 5, remains then largely intact until most other native and non-native helices have been opened under stretching (intermediate 11). Finally, the central domain (II) exhibits a more distributed unfolding fate which extends from intermediates 1 to 10. Hence, mechanical breaking of the native structure does not occur through successive openings of entire native domains. Instead, native helices contribute to a more complex (yet largely reproducible) sequence of force drops, reflecting also the rearrangements of non-native helices. For instance, this is the case for the recorded signal between intermediates 8 and 10 which is largely caused by successive rearrangements of weak non-native helices between 10kT and 15kT (compare magenta and red curves for the native structure on Fig 7 and experimental and simulated curves on Fig 8). In retrospect and more generally, these results underline the possible pitfalls in attempting to assign specific structural features of large RNA molecules by studying the mechanical unfolding of their independently folded domains separately.

Discussion

We have measured the force range to unfold RNA secondary structures by mechanical stretching experiments. It extends from 10 pN for AU rich to 25 pN for GC rich regions in agreement with intermediate values reported for intermediate G+C contents. We also showed that non-native rearrangements have a large influence on force-extension measurements of complex RNA structures, as in the case of E. coli’s 16S rRNA presented here.

Interestingly, this force-induced unfolding process of the bare 16S rRNA’s domains seems to mirror, only in reverse order, the predominant 5’ to 3’ polarity of the in vitro assembly of 16S rRNA into 30S ribosomal subunitsnoller . From a more general perspective, the high reproducibility of the mechanical unfolding curves shown here (e.g., Figs 4-5) sharply contrasts with the multiple folding and misfolding pathways usually experienced by RNA molecules of this size during thermal renaturation. This reflects the fact that unfolding/refolding pathways under mechanical constraint solely explore a restricted number of possible intermediate structures, as compared to unconstrained denaturation/renaturation folding experiments. In other words, single molecule unfolding and refolding experiments under mechanical control probe particular, well-defined pathways due to the slowly varying external constraint applied onto the ends of the RNA molecule. In addition, we found that the overall unfolding curves did not critically depend on the rate of pulling used (typically 300 nm/s); for instance, imposing

an extension rate twice as fast or twice as slow did not significantly modify the force-extension curves (data not shown). In retrospect, this restricted set of unfolding pathways and their relative insensitivity to the precise values of external parameters also explain why we could simulate these force-induced unfolding pathways starting from a given secondary structure, while predicting such 1540-nucleotide initial structure a priori is still beyond the current limitations of secondary structure prediction algorithms.

Despite clear similarities, the agreement between simulated and experimental force-extension responses in Fig 8 is uneven. In fact, variations between predicted (red) and measured (black) curves might reflect real differences between the probed structure and the actual native secondary structure inside the ribosomeramakrishnan ; noller1  used here as the intial structure in the simulations. In particular, deviations at the beginning of the stretching plateau might originate from alternative base pair (re)folding of the 3’ major domain (III) due to the absence of essential ribosomal proteins (e.g., as s7nowotny ) and Mg<sup>2+</sup> ions. Moreover, the relatively short time scales (few seconds) of these stretching-induced unfolding/refolding experiments might not be sufficiently long to let 16S rRNA find its lowest free energy structure between successive pulls.

These results illustrate what should be expected, in general, when RNA secondary structures are probed by mechanical force. Strong helices resist until their breaking exposes weaker regions, which are unable to withstand such high forces. This leads to the unfolding of a significant domain with a concomitant force drop. A fraction of the unpaired bases then typically reform different helices, which compensate, in part, for the sudden relaxation of the mechanical tension. Yet, force-extension responses are not completely smoothed out, as initially suspectedmontanari ; bundschuh ; lubensky , by these local rearrangements. Instead, they reveal the slow dynamics of large scale cooperative changes in complex RNA structures. Tertiary interactions, likely marginal here due to the absence of Mg<sup>2+</sup> ions, are expected to strengthen the unfolding cooperativity between interacting domains and, concomitantly, increase the reformation of non-native helices upon stretching.

Local rearrangements of RNA molecules, similar to those reported here, likely occur in vivo as well, in particular, during translation when large domains of messenger RNAs become unfolded upstream of the ribosome. In fact, the influence of long-lived intermediate structures is likely ubiquitous to the RNA folding problem itself, as slow structural rearrangements are known to occur in the context of both in vitro and in vivo RNA folding processeswilliamson1 ; woodson . New experimental tools are needed to better understand the strategies of RNA molecules in circumventing such kinetic traps (for instance through specific interactions with ions or proteinstreiber , through RNA chaperoneswoodson  or co-transcriptional encoded folding pathwaysisambert ).

By exploring RNA structure energy landscapesfontana , micromanipulations combined with appropriate stochastic simulations can help address such questions, reflecting both structural and metastability features of single RNA molecules.

Added note:

Onoa et alonoa  have recently reported experimental results on the mechanical unfolding of the L-21 derivative of Tetrahymena thermophila ribozyme, a 390-nucleotide catalytic RNA. By contrast with the present study which strictly focuses on the RNA secondary structure level (no Mg<sup>2+</sup> added), Onoa et al primarily investigate the tertiary fold of this selfsplicing ribozyme in the presence of Mg<sup>2+</sup>. A variety of hysteresis responses to the applied force is presented for various parts of the molecule or in the presence of specific antisense oligos. A direct correlation between Mg<sup>2+</sup>-dependent unfolding events and the opening of specific native helices is proposed.

Materials and Methods

Sample displacement: Sample displacement is driven and monitored by a nanometer resolution piezoelectric stage with capacitive position sensor (P530-3, Physik Instrument). The piezoelectric stage position is controlled and monitored by a 0-10V voltage.

Optical tweezer: The optical tweezer consists of a Nd:Yag infra-red laser beam (TOPAZ, SpectraPhysics) focussed inside the capillary by a 1.3 N.A. x100 objective (Zeiss). The laser is always set at full power (2.5 W) and the stiffness of the trap is controlled by the amplitude of an acoustic wave generated by an acousto-optic modulator (A-A) placed right after the laser head. The experiments described here were performed with a 50% attenuation of the laser intensity, which sets the optical trap stiffness around 7$`\times `$10<sup>-5</sup> N/m. The bead displacement from the laser beam focus point is measured as follows: after passing through the sample, the bead diffused light is collected by a 0.6 N.A. x40 objective (Zeiss). The objective back focal plane is imaged by a lense of 40 mm focal length onto a two-quadrant photodetector (S5980, Hammamatsu). The whole experiment setup is mounted on an invar table so as to minimize thermal position drift. The photodiode electric currents $`I_A`$ and $`I_B`$ are converted into voltage and amplified by a home made amplifier. The voltage difference $`V_AV_B`$ which is proportionnal to the distance of the bead away from the trap center is further amplified and filtered at 300 Hz by low noise amplifier (SR-50, Stanford Research Instrument). The total light intensity that is collected by the x40 objective measured by the voltage sum $`V_A+V_B`$ is also amplified.

Data acquisition: The monitoring voltage coming out of the piezoelectric driver, the voltage difference $`V_AV_B`$ and the voltage sum $`V_A+V_B`$ coming out of the low noise amplifier are each directed into a separate channel of an acquisition board (ATMIO-16X, National Instrument). The driving voltage of the piezoelectric stage is generated by the same board. The acquisition rate is 300 Hz which sets the duration of the stretching/relaxing experiment around 10-20 seconds.

Calibration: The fourier power spectrum of a free bead inside the trap follows a lorentzian law as expected for brownian fluctuations. Fitting this curve with two parameters provides both the trap stiffness and the voltage/distance conversion factor. In the case of a pulling experiment, these two parameters are used to convert the ratio $`\left(V_AV_B\right)/\left(V_A+V_B\right)`$ directly into piconewtons. The maximum force that can be measured with our setup is 60 pN. The bead position resolution inside the trap is $`\pm `$ 5 nm which sets the force resolution at $`\pm `$ 0.4 pN. The bead is captured at about 500-1000 nm from the capillary interior surface. The pUC19 dimer contour length is 1742 nm (0.33 nm/pb). The geometry imposes to displace the piezo stage by 48.8-56 nm and 985-1108 nm to completely unfolded the small RNA motifs and 16S rRNA, respectively.

Molecule synthesis and functionalization: RNA molecules were synthesized by in vitro ”run off” transcription of EcorV linearized DNA plasmids. These plasmids were constructed by inserting DNA oligomers (IBA GmBh) starting with a T7 promotor region inside the BamHI-PstI region of pUC19. The RNA sequence was flanked by 12 nucleotides at both extremities to allow for the ligation with the double-stranded DNA arm extensions. In the case of 16S rRNA, the gene was isolated by PCR from pKK3535 plasmid (courtesy of K. Lieberman and H.F. Noller). It was cut by BstEII-BclI and then inserted in pUC19 together with oligomers carrying a T7 promotor with the DNA arm extensions and the complementary ends of the 16S sequence. The reconstructed plasmids were produced in competent DH5 alpha bacteria and were extracted and purified using Jetstar purification kit. They were further sequenced. Due to the small length of RNA molecules, there were extended with digoxygenin or biotin labelled dsDNA at, respectively, the 3’ and 5’ ends to enable grafting between the capilary glass surface and the silica bead. In practice, DNA oligomers (Fig 9) were first ligated to Pst I restricted digoxygenin labelled pUC19 to yield a 12-nucleotide 3’ extension complementary to the RNA 3’ end. Then, the DNA/RNA hybridization and ligation protocol was the following: RNA was heated to 90C for 5 minutes then quenched on ice. It was incubated with the former prepared pUC19 (molar ratio 100/1) at 70C for 20 minutes and then slowly cooled ($``$0.6C/min) to 16C. At this temperature, T4 DNA ligase and buffer were added and the ligation reaction was carried over 4 hours at 16C. The band corresponding to the pUC19 molecular weight on a 0.7% agarose gel was purified using Qiaquick (Qiagen). The whole procedure was repeated with the RNA 5’ end using a BamH I restricted biotin labelled pUC19 DNA (Fig 9). The band corresponding to a pUC19 dimer molecular weight on a 0.7% agarose gel was purified using Qiaquick (Qiagen). In the case of the 16S RNA, the functionalization protocol was slightly modified. The oligomers were first hybridized with the RNA 3’ end following the heat-cooling protocol described above. The excess oligomers were washed away by 2 consecutive centrifugations at 4000 $`\times g`$ and 16C using GS-200 microspin column. The Pst I restricted digoxygenin labelled pUC19 DNA, T4 DNA ligase and buffer were added and the ligation reaction was carried over 4 hours at 16C. The same procedure was repeated with the 5’ end and the molecule was purified on an 0.7% agarose gel by cutting the band corresponding to a pUC dimer molecular weight. Prior to the experiment, the molecules are incubated with the streptavidin coated beads (Bangs Laboratories) for 30 minutes. The solution is introduced in the rectangular capillary by a peristatic pump which allows buffer circulation. All experiments were performed at room temperature and in Tris 10 mM pH 7 NaCl 250 mM buffer.

Acknowledgements

We thank K. Lieberman and H.F. Noller for kindly providing us with the pKK3535 plasmid, D. Evers and R. Giegerich for the use of their “RNAMovies” software, and L. Bourdieu, C. Ehresmann, S. Lodmell, T. Pan, M. Poirier and E. Westhof for discussions and suggestions. This work was supported in part by an ACI “Jeunes Chercheurs” grant from Ministère de la Recherche (France), an NOI grant from the CNRS, and by the “Physique et Chimie du Vivant” program of the CNRS.


format: text

text:
# Non conservation of the magnetization current across magnetic hetero-structures

## Abstract

We show that when the magnetizations on opposite sides of a junction are not collinear, the magnetization current is not conserved as the junction is crossed. Thus the usual treatment of this problem needs to be modified. We argue that this is due to an implicit assumption of an external torque that is required to clamp the magnetization in place. The physical consequence of this is explored.

There is much recent interest in spin polarized transport in heterostructures such as tunnel junctions and metallic multilayers. A general approach to these types of problems is to solve for the steady state solutions. To incorporate the transmission and reflection of the injecting current at the interface, it is necessary to solve the boundary value problem. Central to treating the transmission and the reflection is how the charge current $`J`$ and magnetization current $`J_𝐌`$ go from one side of the heterostructure to the next. For a simple heterostructure with no magnetic impurities both types of currents are usually assumed to be conserved as they go across the barrier. An equation that is usually employed in the treatment of this phenomena is of the formValet

$$J_𝐌=\underset{s}{}(\alpha _s\mu _{Ls}\beta _s\mu _{Rs})$$
(1)
where $`\alpha `$, $`\beta `$ are coefficients, $`\mu _{L(R)s}`$ are the chemical potentials of the electrons of spin channel $`s`$ on the left (right) of the junction. In this equation, it is assumed that $`J_𝐌`$ is unique because it is conserved. Recently, it was discovered that when the magnetizations on opposite sides of a metallic barrier structure are noncollinear, a current can induce torques on the magnetizationsCornell ; Tsoi ; Slon1 ; Berger . This motivated us to examine the boundary condition across a heterostructure more closely. We consider the strong barrier limit corresponding to tunnel junctions and the weak barrier limit with zero interface resistance, corresponding to metallic multilayers. We find that when the magnetizations on opposite sides of the barrier are noncollinear, the magnetization current is not conserved as it goes across the barrier. Thus eq. (1) cannot always be used. An examination of the physics suggests that implicit in the assumption of the noncollinear magnetizations is a torque to maintain the orientation of the magnetization. It is this torque which changes the magnetization current. We now discuss our results in detail.

We first remind the reader the meaning of the direction of magnetization. The general quantum statistical properties of a system is described by the density matrix $`\rho `$. The general state of magnetization is therefore described by the $`2\times 2`$ polarization matrix $`|M_0|\rho `$ where $`|M_0|`$ is the magnitude of the magnetization. As is well knownMerzbacher  there is always a direction $`𝐌`$ so that $`\rho =(1+\sigma 𝐌)`$ where $`\sigma `$ are the Pauli matrices $`(\sigma _x,\sigma _y,\sigma _z)`$. The density matrix is diagonal along $`𝐌`$. We call this the direction of magnetization. Along this direction, the off-diagonal components of the density matrix are zero and there is no phase coherence between the spin up and the spin down states.

We next evaluate the magnetization current. For simplicity we take a one dimensional geometry where the currents are uniform in the $`x`$-$`z`$ plane and move in the $`y`$ direction. The evaluation of the magnetization current is simplest when we pick a coordinate system along the direction of magnetization. The magnetization current then becomes the difference between the spin up and the spin down current. We first consider the strong barrier limit of spin polarized tunnelling. Our starting point is the well accepted tunnelling Hamiltonian $`H_T=T_{ss^{}}c_{Lks}^+c_{Rqs^{}}+c.c.`$ where $`s`$, $`s^{}`$ are with respect to the magnetization directions on the left (labelled by the subscript $`L`$) and on the right (labelled by the subscript $`R`$). $`k`$, $`q`$ are other labels for the electronic states. We are interested in the case when the directions of magnetization on opposite sides of the junction are not the same. The spins states $`|s`$ along one magnetization direction are related to those ($`|s^{}`$) along the other direction by a rotation matrix $`𝐑_{s,s^{}}`$. From a single particle quantum mechanical calculation we obtain an estimate of the tunnelling matrix element

$$T_{ss^{}}=𝐑_{ss^{}}\mathrm{exp}[_0^d𝑑y(\mathrm{}^2(U\mu _{Ls}+\mu _{Rs^{}}Vy)/2m)^{0.5}].$$
(2)
where $`U`$ is the tunnelling barrier, $`V`$ is the external voltage, $`m`$ is the electronic mass. In linear response, the term $`\mu _{Ls}+\mu _{Rs^{}}Vy`$ is absent.

Stimulated by the experimental success in superconductivity tunnelling there has been much studies done with this type of tunnelling Hamiltonian. We shall follow the work by Cohen, Falicov and PhillipsCFP , as is explained in detail in MahanMahan . The magnetization current on the left can be evaluated from the time rate of change of the left magnetization due to the tunnelling Hamiltonian: Since $`J_s=N_s/t`$, we get $`J_M^L=(N_+^LN_{}^L)/t=M_L/t`$ (we have used a unit so that the Bohr magneton is 1) with a similar expression for $`J_M^R`$. The time derivative of the particle number of spin $`s`$ can be evaluated as its commutator with the tunnelling Hamiltonian. The details of this for each spin component is essentially the same as in previous calculations and thus will not be repeated here. We find thatcf

$$J_M^L=4\pi e\underset{ss^{}}{}s|Tss^{}|^2(n_L(\mu _{Ls})n_R(\mu _{Rs^{}})),$$
(3)
and

$$J_M^R=4\pi e\underset{ss^{}}{}s^{}|Tss^{}|^2(n_L(\mu _{Ls})n_R(\mu _{Rs^{}})).$$
(4)
where $`\mu _{Ls}`$ ($`\mu _{Rs^{}}`$) is the chemical potential for the spin $`s`$ ($`s^{}`$) electrons on the left (right) side. $`n(\mu )`$ is the density of particles at energy $`\mu `$.

When the directions of magnetization on both sides of the junction are collinear, $`T`$ is diagonal and we find that $`J_M^L=4\pi e_ss|Tss|^2(n_L(\mu _{Ls})n_R(\mu _{Rs}))=J_M^R.`$ The magnetization current is conserved. When the magnetization directions are not collinear, $`T`$ is no longer diagonal. In the linear response regime, $`T_{s,s^{}}=T_{s^{},s}`$ is a symmetric matrix. After renaming the variables of summation, one has $`J_M^L=4\pi e_{ss^{}}s^{}|Tss^{}|^2(n_L(\mu _{Ls^{}})n_R(\mu _{Rs}))`$. $`J_M^LJ_M^R`$ because $`\mu _{Ls^{}}\mu _{Ls}`$ for $`ss^{}`$. The latter comes about because there is a splitting of the spin up and spin down chemical potential caused by the spin accumulation effect6 ; 7  and the interplay between the charge and magnetization dipole layers1 .

$`J_M^L`$ and $`J_M^R`$ correspond to magnetizations along different directions. In the above we have shown that the magnitude of these two currents are not the same. A more appropriate condition would be requiring the magnetization currents for magnetizations along the same direction be equal: $`J_M^L=J_M^R`$. This requirement obviously cannot be satisfied because for example, on the right hand side it would require a magnetization current with a magnetization component that is perpendicular to the magnetization on the same side. From what we just derive, this is zero. The reason for the nonconservation of the magnetization current is easy to understand. Consider a wave function $`|sL`$ of an electron of spin $`s`$ on the left. With respect to the magnetization direction on the right, this wavefunction can be written as $`_s^{}𝐑_{ss^{}}|s^{}R`$. The corresponding density matrix for a state like this is given by

$$\rho _0=\left(\begin{array}{cc}R_{++}^2& R_+^2\\ R_+^2& R_{}^2\end{array}\right)$$

with non-zero off-diagonal elements. But by our definition of the magnetization direction, we have assumed that after this electron has moved from the left to the right, the phase coherence between the up and down states are lost and the off-diagonal components become zero. This happens because we have implicitly assumed that there is some external torque that clamps the magnetization in certain directions. This torque can come from the crystalline anisotropy energy of the system, for example. It is this implicit torque that causes the magnetization current not to be conserved. In addition to this torque $`t_0=(J_M^RJ_M^L)/a\mu _B`$ ($`\mu _B`$ is the Bohr magneton, $`a`$ is atomic spacing) there may be other contributions which we shall not discuss in this paper. The reasoning in this section also applies to metallic multilayers. We turn our attention to that case next.

For metallic multilayers, eq. (3) and (4) are no longer valid. For simplicity we consider the other limit where the interface resistance is zero. For that case, the current is controlled by the voltage drop along the metal. We assume that the length of heterostructure is less than the spin diffusion length and apply a parallel circuit model to illustrate the essential physics. In that case, we get $`J_s^{L(R)}=\sigma _s^{L(R)}E^{L(R)}`$ where $`\sigma _s^{L(R)}`$ is the conductance for spin $`s`$ on the left (right). $`E`$ is the electric field that drives the current. The spin current on opposite sides of the heterostructure are related by $`J_s^{}^L=_s|𝐑_{s,s^{}}|^2J_s^R`$, where, again, $`𝐑`$ is the rotation matrix. There are two ways that one can see this. The current on the right for spin index $`s^{}`$ is just the number of electrons with $`z`$ angular momentum $`s^{}`$ moving past per unit time per unit area. These electrons come from the left side. The factor $`R^2`$ takes into account the probability that the $`z`$ angular momentum is $`s^{}`$. Formally, in terms of electron creation and destruction operator the current operator is just $`J_s^{}^R=vc_s^{}^+c_s^{}`$ where $`v`$ is the velocity. Now the creation and destruction operator can be expressed with the z axis along $`M_L`$ as $`c_s^{}=𝐑_{s^{},s}c_s`$, $`c_s^{}^+=𝐑_{s^{},s}c_s^+.`$ Because $`𝐑`$ is real, it occurs in the expressions for both the creation and destruction operator. By our assumption of the magnetization direction, we have set such terms as $`<c_s^+c_t>`$ with $`st`$ to be zero. Hence the final expressions quoted above is obtained.

It is straightforward to show that $`E^{L(R)}=\sigma ^{L(R)}E/(\sigma ^L+\sigma ^R)`$ where the total conductance $`\sigma ^{L(R)}=_s\sigma _s^{L(R)}`$, $`E=E^L+E^R`$ measures the total voltage drop. Just as explained in the previous paragraph, because the magnetization points in different directions on opposite sides of the heterostructure, the magnetization currents are not equal to each other. The magnitudes of the two currents are equal to each other, however.

In the presence of a finite interface resistance, in linear response the current is linearly proportional to the chemical potential. We thus expect the physics to be governed by equations of the form

$$J_M^L=\underset{s}{}r_{1s}^L\delta \mu _s^Lr_{2s}^L\delta \mu _s^R$$
(5)
$$J_M^R=\underset{s}{}r_{1s}^R\delta \mu _s^Lr_{2s}^R\delta \mu _s^R$$
(6)
where the $`r`$’s are the effective interface resistances. We next explore the physical implication of our result.

In typical experimental arrangements, only the directions of the magnetizations far away from the interface are specified. Because of the torque, these directions will be rotated as the interface is approached. A domain wall is created. We estimate this rotation here for a simple case. We assume that far from the interface the magnetzations are in the $`x`$-$`z`$ plane along the easy axis of anisotropy at angles $`\theta _L`$, $`\theta _R`$ so that $`\theta _R\theta _L=\mathrm{\Phi }`$. At the interface, the left (right) magnetization has rotated by an angle $`\delta \theta _R`$ ($`\delta \theta _L`$). From minimizing the energy, the equation governing the angle $`\delta \theta `$ at a distance $`y`$ away from the interface is given by $`J_id^2\delta \theta /dy^2K_i\delta \theta =0`$ where $`i=L,R`$; $`J_i`$ ,$`K_i`$ are the exchange and anisotropy constants. The effective field $`H_{eff}`$ due to the exchange, $`J_id^2\delta \theta /dy^2`$, is balanced out by that due to the anisotropy, $`K_i\delta \theta `$. The solution of this equation is $`\delta \theta (y)=\delta \theta _iexp(|y|/d_i)`$ where the magnetic length $`d_i=\sqrt{J_i/K_i}`$. At the interface, the effective field from the exchange becomes $`Jd\delta \theta /dy`$ while that from the anisotropy remain the same. The net torque $`t=E/\theta =M_0\delta \theta _i(K_iJ_i/d_i)=M_0\delta \theta _i[K_i(J_iK_i)^{0.5}]`$. This creates the difference in the magnetization current discussed above. We thus get $`t=t_0`$. $`\delta \theta _i=\mathrm{\Delta }J_M/M_0a[K_i(J_iK_i)^{0.5}]`$. Thus we can relate $`\delta \theta _i`$ to $`\mathrm{\Delta }J_M`$.

Experimental measurements are carried out under constant charge current conditions. It is necessary to relate $`\mathrm{\Delta }J_M`$ to the current $`J`$. To achieve this goal it is necessary to find solutions for the currents for each side of the junction and match them across the junction. In the treatment of spin polarized transport, we should use eq. (3) and (4) instead of eq. (1). These equations will be complemented by equations which relate the currents to the chemical potentials inside the material on either side of the junction. An example of such a equation is the modified Landau Gilbert equation with a source term equal to the divergence of the magnetization current:EB

$$\frac{𝐌}{t}\gamma 𝐌\times 𝐇+\widehat{𝐉}_\mathrm{M}=\frac{\delta 𝐌}{\tau }$$
(7)
where $`\gamma `$ is the gyromagnetic ratio, $`\widehat{𝐉}_\mathrm{M}`$ is the spin current (tensor), $`\tau `$ is the spin relaxation time, and $`𝐇`$ is the effective field describing the precession of the magnetic moments. The magnetization and charge currents are given by Fick’s law as sums of terms proportional to the gradients of the charge ($`\delta n`$) and magnetization ($`\delta 𝐌`$) densities and the local internal electric field. They are

$`𝐉`$ $`=`$ $`\sigma V\beta D\delta 𝐌𝐩_0\alpha _1D\delta n`$ (8)
$`\widehat{𝐉}_\mathrm{M}`$ $`=`$ $`\beta \sigma V𝐩_0D\delta 𝐌\alpha _2\beta D\delta n𝐩_0`$ (9)
where $`\alpha _1`$ and $`\alpha _2`$ are two phenomelogical parameters, $`\beta `$ measures the asymmetry of the spin up and spin down current in the metal, $`\sigma `$ is the conductivity of metal, and $`D`$ is the diffusion coefficient. $`𝐩_0`$ is the unit vector along the direction of magnetization. $`V`$ is the sum of the electic potential $`V_e`$ describing the external electric field and the local electric (screening) potential $`W`$ due to other electric charges determined self-consistetly

$$W(𝐫)=d^3𝐫^{}U(𝐫𝐫^{})\delta n(𝐫^{}).$$

From these equations and the charge conservation

$$𝐉=\frac{\delta n}{t}$$
(10)
we obtain for each side of the heterostructure the charge and magnetization dipole layers at the interface given mathematically by

$$\begin{array}{cc}\delta n^R=\underset{i=1}{\overset{4}{}}\delta n_{i0}^Re^{(y\frac{d}{2})/l_i},\hfill & \\ \delta 𝐌^R=\underset{i=1}{\overset{4}{}}\delta 𝐌_{i0}^Re^{(y\frac{d}{2})/l_i}\hfill & \end{array}$$
(11)
where the four length scales correspond to the screening length and the magnetic lengths for the longitudinal and left and the right rotating transverse magnetizations. $`d`$ is the thickness of the barrier. The superscript $`R`$ denotes the right hand side. A similar expression can be written for quantities on the left hand side. Because of the coupling between the charge and spin degrees of freedom, the length scales $`l_i`$ are renormalized from their bare value. The coefficients $`\delta n_{i0}`$ are found to depend on $`\delta 𝐌_{i0}`$, while $`\delta 𝐌_{i0}`$ are determined by matching the boundary condition at the interface. We illustrate this approach with a calculation of the steady state solutions in the low field linear response limit. After this boundary matching is carried out, we found for metallic multilayers in linear response a large change in the longitudinal component of the magnetization that is of the order of the voltage times the density of states at the Fermi surface (in units with $`\mu _B=1`$). The transverse magnetization change is much smaller. As a result, there is no additional torque generated along the ferromagnet. All the torque comes at the interface.

When the material parameters on both sides are the same, we found

$$\begin{array}{cc}e\mathrm{\Delta }J_M=\beta J\mathrm{cos}(\mathrm{\Phi }/2)\mathrm{sin}\mathrm{\Phi }\hfill & \\ \times \frac{A_1+2\sqrt{1\beta ^2}\xi \left[1+A_2\mathrm{sin}^2(\mathrm{\Phi }/2)/(1\beta ^2)\right]+\xi ^24\mathrm{sin}^2\mathrm{\Phi }/2}{1\beta ^2+2\sqrt{1\beta ^2}\xi +\xi ^2\mathrm{sin}^2\mathrm{\Phi }}\hfill & \end{array}$$
(12)
where $`\xi `$ is the ratio of the conductance of the interface to that of the metal. $`A_1`$, $`A_2`$ are some dimensionless constants of the order of unity.

In summary, we argue that the magnetization current is not conserved at the interface of a heterostructure when the magnetizations on opposite sides of the structure are not collinear. This is due to the implicit assumption of a torque which creates a domain wall at the interface. This is illustrated with a steady state solution of the of the problem in the limit of linear response. In that limit, all the torque comes at the interface, no additional torque is introduced in the metal.

Z L was supported in part by China NNSF.


text:
# Flow of He II due to an Oscillating Grid in the Low Temperature Limit

## Abstract

The macroscopic flow properties of pure He II are probed in the limit of zero temperature using an oscillating grid. With increasing oscillation amplitude the initially pure superflow changes abruptly: the resonant frequency decreases and the response becomes strongly nonlinear, attributable to the growth of a boundary layer of quantized vortices that increases the effective mass of the grid. On further increase of oscillation amplitude, the flow undergoes a transition to turbulence.

preprint: APS/123-QED

Although the exotic flow properties of He II have been a subject of intensive investigation ever since the discovery of superfluidity, accumulating a vast amount of experimental data and theoretical knowledge Tough ; Newton , many important features still remain to be explained. In the limit of low velocity, He II flow is very well described in terms of Landau’s two-fluid model. On increasing the flow velocity beyond a certain threshold, however, quantized vortices appear. Their presence couples the originally independent normal and superfluid velocity fields via a mutual friction. Above $``$1.2 K, where He II contains an appreciable proportion of normal fluid, numerous investigators observed that on exceeding a suitably defined Reynolds number, He II flow acquires an increasingly classical character: (i) the surface of a rotating bucket forms a nearly classical meniscus Osborne:50 ; (ii) flow past a microsphere displays both laminar and turbulent drag Wilfried , as well as (iii) a drag crisis vSci ; (iv) the energy spectrum of turbulent He II involves an inertial range with a classical K41 Kolmogorov roll-off exponent of $`5/3`$ Tabeling ; (v) the decay of both grid-generated and counterflow turbulence displays a classical character Skr ; Skr1 ; Vinen ; QTReview . Although these features are seen over a wide range of normal fluid/superfluid density ratios, it is impossible to exclude the possibility that this classical-like behaviour is caused by the presence of the viscous normal fluid. Following oscillating experiments with a tiny sphere Wilfried  and a very thin vibrating wire Wire , yielding interesting results attributable to single vortices, there is a clear call to study the macroscopic flow properties of He II in the $`T0`$ limit where normal fluid is (almost) absent and pure superflow can be investigated systematically.

Our tool for performing such investigations is the thin electroformed nickel grid grid  of $`70`$% transparency and mesh size of $`127`$ $`\mu `$m shown schematically in Fig. 1. It constitutes a circular membrane $`2R=8`$ cm in diameter, tightly stretched on a circular mild steel carrier, mounted horizontally, equidistant from two gold-plated plane copper electrodes, each drilled with 170 holes of 2 mm ID to connect the grid space to the rest of the experimental chamber. The electrodes are separated by $`d=2`$ mm, immersed in a 1.5 liter volume of isotopically pure 4He in an experimental cell attached to the mixing chamber of the dilution refrigerator. The apparatus can be operated at all pressures up to the solidification pressure. A static potential, typically $`V_0=500`$ volts, is applied to the grid. A driving potential $`V_1=V_{10}\mathrm{cos}\omega t`$ ($`V_{10}`$ $`V_0`$) applied to the upper electrode produces a driving force on the grid of form $`f_d=2\epsilon _0\epsilon _r\pi R^2VV_1/d^2`$, where $`\epsilon _0`$ and $`\epsilon _r`$ denote respectively the relative permittivities of free space and the liquid 4He. The grid can thus be considered as an oscillating membrane Morrell  under uniform tension. In this study we approximate its motion as one-dimensional and assume that the oscillation amplitude is uniform across its area. Oscillations of amplitude $`\mathrm{\Delta }D`$ induce a signal of amplitude $`V_2=V\mathrm{\Delta }D/(2d)`$ on the lower electrode, which can be monitored with a lock-in amplifier (which we use to investigate the low-drive linear response of the grid) or directly with a memory oscilloscope, allowing visualization of transient processes resulting from the strongly nonlinear response at higher drives. Allowing for a reduction in the induced voltage $`V_2`$ by a factor of $`(1+C_c/C)^10.065`$, where $`C_c700`$ pF is the capacitance of the connecting cable and $`C47`$ pF is the capacitance between the grid and the lower copper electrode, the response amplitude $`|V_2|`$ provides a direct measure of the amplitude and spatially averaged peak velocity $`|v_g|=|\omega \mathrm{\Delta }D|`$ of the oscillating grid.

At the lowest oscillation amplitudes, the grid displays linear behaviour in that its frequency response to the drive is a Lorentzian of narrow width caused predominantly by nuisance damping, with a quality $`Q`$ factor typically exceeding 5000. The resonant frequency $`f_{\mathrm{res}}`$ can be altered temporarily by increasing the pressure, which we attribute to quantized vortices generated e.g. by the jet from the filling capillary: they probably become pinned to inhomogeneities on the grid surface, and between the grid and surrounding electrodes. Moving the grid violently at high drive amplitude, however, is found to stabilize $`f_{\mathrm{res}}`$ typically to within $`\pm 0.1`$ Hz, presumably by shaking off most of this pinned vorticity.

Fig. 2 shows our central experimental observation, which for the “cleaned” grid has been repeated for several pressures spanning 0.3 $`p`$ 24.8 bar without appreciable change. As the drive increases and the grid amplitude reaches the first threshold (typically 10–20 mVpp, corresponding to a mean grid velocity of $`0.3<v_g^{(1)}<0.6`$ cm/s), the oscillation amplitude at resonance continues to rise in proportion to the drive (Fig. 3), but the resonant frequency suddenly starts decreasing and the resonance curves acquire highly nonlinear features. If one follows e.g. the resonance curve for the 0.05 Vpp drive down from 1091 Hz (Fig. 2) by slowly sweeping the drive frequency (in practice, digitally, in steps of 0.001 Hz), the system displays the usual nearly Lorentzian stationary response while the amplitude remains below the first threshold. On further decreasing the drive frequency, the response amplitude increases above this first threshold until, shortly below 1086 Hz the amplitude suddenly collapses down to a lower stable branch. On increasing the drive frequency again from this point the system stays on the lower branch until about 1087.15 Hz, where a transition to the stable upper branch occurs. These hysteretic loops are robust to temperature increase, at least up to our maximum of 130 mK.

Within the hysteretic parameter range, beat phenomena are occasionally seen between two apparently stable amplitude levels, with envelopes of period $``$1 s (see, eg., Fig. 2, drive 0.05 V$`pp`$ around 1086.5 Hz); once established, they are stable on the scale of hours. Small changes of frequency do not kill them, but modify the upper and lower amplitude levels between which beating occurs. Only with further change of the driving frequency does the beating disappear. On restoring the original frequency the response is not completely reproducible; within some small frequency range beating might not appear at all; and sometimes it reappeared only later, or not at all. Despite considerable effort, we did not succeed in establishing any fully repeatable pattern or well-defined conditions for the appearance of beats beating .

The fall in frequency with increasing drive (Fig. 2) reaches typically 2 Hz for all investigated pressures (see Fig. 4), ceasing at a second threshold amplitude (typically $`200`$ mVpp corresponding to a mean grid velocity of about $`v_g^{(2)}6`$ cm/s). Observation of the two well-defined resonant frequencies, at all investigated pressures shifted by $`2`$ Hz, is a remarkable feature of the superflow that, to our knowledge, has not previously been reported. With further increase in drive, the oscillation amplitude at resonance initially remains almost constant, while the widths of the resonance curves increase rapidly (see Fig. 2). Only for drive levels exceeding by about an order of magnitude that for the second threshold does the amplitude at resonance grow again; this time approximately in proportion to the square-root of the drive, as shown in Fig. 3. In this high-drive regime the linewidth increases rapidly while the resonant frequency decreases gradually, qualitatively in the manner expected for increasing damping. Pronounced beatings were never observed in this high-drive regime. If the oscillation amplitude is measured as a function of drive level at fixed frequency, one observes clear hysteresis within the frequency range containing the two stable branches of the grid response (Fig. 3). Outside this range, the drive dependences are nonlinear, but without any hysteresis.

Measurements of the grid response at low temperature in vacuum (smooth curve in Fig. 2) demonstrate that the presence of the He II has two effects. First, it down-shifts the low-drive resonant frequency by about 30 Hz from its vacuum value of $`f_0=(1117.2\pm 0.05)`$ Hz, depending on the pressure $`p`$ in the cell (see Fig. 4). This classical effect arises from hydrodynamic enhancement of the mass $`M`$ of the grid by $`\mathrm{\Delta }M_H=\beta \varrho _{He}(p)M/\varrho _G`$, grid ; Brooks . The measured $`f_0`$ and low drive resonant frequencies yield $`\beta =3.01\pm 0.05`$. Secondly, it is indeed the He II that is responsible for the observed nonlinearities vacuum .

We believe that the behaviour of the oscillating grid ought to be understandable in terms of quantized vortices generated in He II by its motion. As one possibility, we may speculate that, on exceeding the first threshold, a “boundary layer” of vortex loops (perhaps of horseshoe form kus ; Samuels ) grows on the grid, thus increasing its effective mass. Correspondingly, the resonant frequency shifts down and the resonance curves acquire strongly nonlinear features. Provided the vortex loops remain pinned and that they do not reconnect to create free vorticity, the mass enhancement would be non-dissipative, as observed (see below). We note that small mass enhancements were observed earlier for a sphere vibrating in He II above 1 K hotsphere , but without associated critical thresholds.

It is interesting to characterize the first threshold by a superfluid Reynolds number $`Re_s=U_{\mathrm{max}}G/\kappa `$, where $`U_{\mathrm{max}}`$ stands for a peak critical flow velocity through the grid window vmax , $`G`$ is its linear size and $`\kappa `$ denotes the circulation quantum. Observed values of $`Re_s10`$ compare well with the critical $`Re_s=UD/\kappa 20`$ ($`U`$ is the transport velocity and $`D`$ the diameter of the pipe) found as a temperature independent threshold in pipe flow of He II at much higher temperature when the flow of the normal component was inhibited by superleaks at both ends of the pipe ToughPipe . The generation of quantized vortex loops on the surface of the grid by macroscopic superflow around it probably involves growth from remanent vortex lines Awschalom , or from a “plasma” of half vortex rings kus . This idea is supported by the following observation. While measuring the response dependence on increasing the drive amplitude at the resonant frequency, the system encounters a nucleation problem when passing the first threshold: on some occasions the resonance response stopped growing with increasing the drive level, and jumped on the response/drive curve (Fig. 3) only later. With decreasing drives this feature disappears, and the response remains proportional to the drive level.

The frequency down-shift $`\mathrm{\Delta }f(p)=f_{\mathrm{res}}f_{\mathrm{sh}}`$ 2 Hz between the two critical thresholds might be considered in terms of a boundary layer of thickness $`\lambda `$ that enhances the hydrodynamic effective mass of the grid by $`\mathrm{\Delta }M_g^\lambda =A\lambda \varrho _{He}(p)`$, where $`A`$ denotes the surface area of the grid. Requiring that the downward shift of the resonance frequency corresponds to those observed experimentally leads to

$$\lambda =\frac{M+\mathrm{\Delta }M_g^\lambda }{A\varrho _{He}}\left(\frac{f_{\mathrm{𝐫𝐞𝐬}}^2}{f_{sh}^2}1\right)=\sqrt{\frac{2\nu _{\mathrm{𝐞𝐟𝐟}}}{\omega }}=0.53\pm 0.05\mu m$$

which would correspond to an effective kinematic viscosity $`\nu _{\mathrm{𝐞𝐟𝐟}}10^5`$ cm2/s. Extending the analogy to classical fluids, we may estimate the expected resonant linewidth due to the drag of such a hypothetical viscous fluid. A straightforward calculation LL  based on estimating the flow velocity gradient in the direction inside the fluid by $`v_g/\lambda `$ leads to a linewidth $``$1 Hz, consistent with resonances for (e.g. Fig. 2, 0.2 and 0.5 Vpp data) near the second critical threshold. However, this picture seems unable to account for the re-entrant character of the nonlinear resonances because the resultant viscous broadening would tend to smear out the effect in question.

We infer that, beyond the second critical threshold (Fig. 2), the vortex loops start to reconnect. Free vorticity is then shed into the liquid, corresponding to dissipation of the vibrational energy and leading to the observed increase in linewidth. Under these flow conditions, He II behaves in close analogy with a classical Navier-Stokes fluid: the square-root behaviour (Fig. 3) of the resonant response as a function of drive amplitude in this range is typical of classical turbulent drag scaling, and it is therefore likely that this threshold marks the onset of turbulence. After leaving the grid, the quantized vorticity probably “evaporates” 3HeBarenghi  as recently proposed in connection with turbulence in superfluid 3He-B Fisher .

The particular challenge posed by these results is to develop a quantitative description of vortex dynamics in the macroscopic flow around the moving grid, showing how it can give rise to the observed amplitude-dependent mass enhancements and re-entrant resonance curves, without dissipation. But the problem also carries wider significance, and will repay careful study, because the dynamical processes in question are probably fundamental to the generation of quantum turbulence.

We are grateful to N.J. Fullwood for making the electron micrograph in Fig. 1, and to D.I. Bradley, A. M. Guénault, R.P. Haley, M. Krusius, D.G. Luchinsky, G.R. Pickett, P. Skyba and W.F. Vinen for fruitful discussions. The research was supported by the Engineering and Physical Sciences Research Council (UK) and the Czech Grant Agency under grant GAČR 202/02/0251.


text:
# 1 Original image from the Advanced Camera for Surveys (ACS) on board the Hubble Space Telescope. The image shows an edge-on spiral and a dwarf galaxy with extended tails of starlight. The image is 70" × 72" in size.

Galaxy Disruption in a Halo of Dark Matter

a

Duncan A. Forbes1, Michael A. Beasley1, Kenji Bekki2, Jean P. Brodie3, Jay Strader3

1Centre for Astrophysics & Supercomputing, Swinburne University of Technology, PO Box 218, Hawthorn, VIC 3122, Australia
2School of Physics, University of New South Wales, Sydney NSW 2052, Australia
3Lick Observatory, University of California, Santa Cruz, CA 95064, USA

Abstract
The relics of disrupted satellite galaxies around the Milky Way and Andromeda have been found, but direct evidence of a satellite galaxy in the early stages of being disrupted has remained elusive. We have discovered a dwarf satellite galaxy in the process of being torn apart by gravitational tidal forces as it merges with a larger galaxy’s dark matter halo. Our results illustrate the morphological transformation of dwarf galaxies by tidal interaction and the continued build-up of galaxy halos.

A long standing problem with the theory of a hierarchical universe consisting of cold dark matter halos (1,2) is that it over-predicts the number of small dwarf satellite galaxies (3). Some satellites have been accreted by their host galaxy, and the stars from these disrupted satellites have been found in the Milky Way (4,5,6) and Andromeda galaxies (7). However, the accretion time for most satellites is sufficiently long that we would expect many of them to survive to the present day (8). As a satellite orbits within a halo it will be subject to tidal interaction effects which cause its outer stars to be stripped away forming extended tails (9, 10). These tails trace out the orbit of the dwarf, becoming a relic stream within the larger galaxy’s halo. Dwarfs containing a disk will have their disk disrupted. A tidal encounter may also create a bar which drives gas to the galaxy center, inducing a burst of star formation. This raises the central surface brightness and makes the galaxy more compact.

As part of its early release observations, the new Advanced Camera for Surveys (ACS) installed on the Hubble Space Telescope obtained a deep multi-color image of UGC 10214 (the Tadpole galaxy) in April 2002 (11,12). The ACS image also contained an edge-on spiral galaxy $``$1.5$`^{^{}}`$ northeast of UGC 10214. The spiral is barely visible on the Digital Sky Survey and is uncatalogued in the NASA Extragalactic Database. We have determined a position of $`\alpha `$ = 16:06:11, $`\delta `$=+55:26:57 (J2000). To the north of the spiral, at a projected separation of 6.4" and close to its minor axis, lies a small galaxy with extended tails of starlight (Fig. 1). There is also an indication that the spiral galaxy disk is warped (Fig. 2), a feature that may be caused by an interaction with a satellite galaxy. Although suggestive, confirmation of a physical association requires that the dwarf galaxy and the large spiral have a similar redshift.

During an observing run in March 2003 on the Keck telescope, we obtained a short exposure image in which the spiral, the dwarf and its tails were all visible. A spectrograph slit was aligned across the nuclei of both galaxies. Spectra with a 1 hour exposure time were taken and reduced using standard procedures (13). We obtained a spectrum of the dwarf with a signal-to-noise ratio of $``$5 and measured a heliocentric recession velocity of 43,445 $`\pm `$ 225 km/s (or redshift z = 0.145). For the spiral we measured 43,433 $`\pm `$ 99 km/s. The two galaxies have statistically indistinguishable velocities confirming their physical association. Assuming a $`\mathrm{\Lambda }`$CDM Universe (with $`\mathrm{\Omega }_\mathrm{\Lambda }`$ = 0.7 and Hubble constant Ho = 75 km/s/Mpc), this redshift implies that the two galaxies have a projected separation of 16 kpc (1" $``$ 2.5 kpc).

The ACS images consist of deep exposures in the filters F475W (g$`^{^{}}`$), F606W (broad V) and F814W (I). They have been aligned to within 1/3 of a pixel (each ACS pixel is $``$0.05"). Using zero points on the AB system (11), and a K-correction for redshift, we derive rest-frame magnitudes and surface brightnesses (after correcting for a (1+z)4 dimming effect). No correction for Galactic extinction is applied as it is small (AV $`<`$ 0.03). With an aperture specially-designed to match the shape of the spiral galaxy, we estimate its total Johnson V band magnitude to be 18.0. This corresponds to an absolute magnitude MV of –21.0 (correction to a face-on magnitude could make the galaxy more luminous by about one magnitude). Assuming an intrinsic disk flattening of 0.1, we estimate that the spiral is of type Sb and is inclined at 70 $`\pm `$ 5o to the line-of-sight.

The dwarf galaxy consists of an elongated main body, which has a diameter of $``$1.5" (3.8 kpc), from which the tails extend. The northern tail is better defined. After about 7.8" (20 kpc) it bends and extends an additional 7.5” (19 kpc). This tail has an average V band surface brightness of 26 mag/sq. arcsec. The southern tail is visible for 6.4” (16 kpc) where it becomes confused with the spiral galaxy starlight. There is a hint that it appears on the other side of the spiral (Fig. 2) with a surface brightness of 26.5 mag/sq. arcsec. For the main body of the dwarf, we estimate an absolute magnitude MV = –16.0. The total luminosity of the tails is more difficult to quantify, but we estimate that there is at least as much starlight in the tails as in the main body of the dwarf. Thus the original galaxy may have been twice as luminous as it is now, i.e. MV $``$ –16.8. From the local galaxy scaling relation, this magnitude would correspond to a mean metallicity of \[Fe/H\] $``$ –0.9 (14).

In the dwarf’s central regions its color remains fairly constant with radius at B–V = 0.34 and V–I = 0.5. These relatively blue colors suggest the presence of young stars. We compared the dwarf colors to a stellar population model (15) to determine its age. The model indicates a mean metallicity \[Fe/H\] = –0.7 with a range of 0 to –1, and a luminosity-weighted age of 2 $`\pm `$ 1 Gyr. Thus, the model suggests that the dwarf had a burst of star formation a few billion years ago. The tails are 0.1-0.2 magnitudes bluer than the main body of the dwarf.

We have fitted the dwarf galaxy isophotes with the IRAF program ellipse. The resulting surface brightness, ellipticity and position angle profiles of the model fit (Fig. 3) show that the dwarf has a position angle twist from near 0o at the center to about –20o in its outer parts. The ellipticity also varies over this radial range from nearly circular to elliptical at radii where the tails begin to dominate the starlight. A Sersic profile (16) has been fit to the galaxy surface brightness within the main body (further out, the galaxy has excess light compared to a Sersic profile). The fit for each filter gives similar results. For the F606W filter we measure a Sersic n value of 1.28 $`\pm `$ 0.08. This value would suggest a central black hole mass of $``$5 x 106 solar masses from the local galaxy correlation (17). We measure an effective radius (Re) of 0.37 $`\pm `$ 0.03" (0.93 +/- 0.08 kpc) and the surface brightness in the V band at this radius ($`\mu _e`$) of 22.7 $`\pm `$ 0.1 mag/sq. arcsec. The central surface brightness ($`\mu _o`$) is estimated to be 20.8 $`\pm `$ 0.2 mag/sq. arcsec. These structural and photometric properties for the dwarf main body resemble those of dwarf elliptical (dE) galaxies (18, 19).

Fig. 4 compares the photometric properties of the dwarf main body with nearby dwarf galaxies. It lies in a region of the diagram occupied by dIrr and dE galaxies. As the starburst in the dwarf fades and it continues to lose mass from tidal disruption it will move towards the region of the diagram occupied by dSph galaxies (20). The figure also shows the location of NGC 205, which has similar photometric properties to our ACS dwarf galaxy. This satellite of M31 is classified as a dE/dSph with a bright nucleus. It also reveals evidence for a young starburst and twisted outer isophotes (21, 22). These features suggest that NGC 205 has experienced gas inflow which induced a nuclear starburst and has undergone tidal shredding of its outer regions. The location of the Sagittarius dSph, which is currently undergoing its final disruption and accretion by our Galaxy (4), is also shown.

The relatively high surface brightness of the tails suggests that the progenitor galaxy contained a disk (10). This is further supported by our estimates of a young starburst, which requires gas and hence a disk. Thus the progenitor may have been a dwarf irregular (dIrr) galaxy. Simulations of galaxy interactions have shown that dIrr galaxies can be tidally disrupted in the halos of large galaxies and transformed to resemble dE and dSph galaxies (10). This process also explains why dE and dSph galaxies are preferentially found in the extended halo of giant galaxies, while dIrr galaxies are generally located in isolated regions (23).

Whether or not the dwarf fully merges with the host spiral depends on the dynamical friction time-scale. This can be estimated, based on the observed luminosity of the dwarf and spiral, and assuming a mass-to-light ratio of 10 for both galaxies. The main uncertainty in estimating this time-scale is the true separation between the galaxies at the start of the dwarf’s orbit. If we assume a separation of 50 kpc (corresponding to the extent of the northern tail) then the time-scale to fully merge is twice the current age of the Universe. As the dwarf loses mass by tidal disruption the dynamical friction time-scale becomes even longer, slowing the orbital decay further.

Simulations of dwarf satellites suggest they have highly elongated orbits (24) and typically produce tails of low surface brightness, e.g. 28-30 mag/sq. arcsec (9,10). So the satellites spend much of their orbit at a large distance from the host galaxy and generally produce tails that are extremely difficult to detect. Here we are witnessing a dwarf that is currently close to its host galaxy and has relatively high surface brightness tails from its disk stars.

Until now, clear observational evidence for a dwarf satellite actually in the process of being tidally disrupted within the halo of a larger galaxy was lacking. Our results indicate that spiral galaxy halos are still being built hierarchically as recently as 2 Gyrs ago (the look-back time for z = 0.145), providing further evidence for the diversity of the stars in galaxy halos (25). Our results also provide observational support for the suggestion that dIrr galaxies can be morphologically transformed into dE and dSph galaxies (20).

References and notes
1. White, S., Rees, M. Mon. Not. R. Astron. Soc. 183, 341-358 (1978).
2. White, S., Frenk, C. Astrophys. J. 379, 52-79 (1991).
3. Kauffmann, G., White, S., Guiderdoni, B. Mon. Not. R. Astron. Soc. 264, 201-218 (1993).
4. Ibata, R., Gilmore, G., Irwin, M. Nature. 370, 194-198 (1994).
5. Lynden-Bell, D., Lynden-Bell, R., Mon. Not. R. Astron. Soc. 275, 429-442 (1995).
6. Helmi, A., Nature. 412, 25-26 (2001).
7. Ibata, R., Irwin, M., Lewis, G., Ferguson, A., Tanvir, N. Nature. 412, 49-52 (2001).
8. Colpi, M., Mayer, L., Governato, F. Astrophys. J.. 525, 720-733(1999).
9. Bekki, K., Couch, W., Drinkwater, M., Gregg, M. Astrophys. J.. 557, 39-42 (2001).
10. Mayer, L. et al. Astrophys. J. 559, 754-784 (2001).
11. Tran, H. et al. Astrophys. J. Letters 585, 750-755 (2003).
12. de Grijs, R., Lee, J., Clemencia, M., Fritze-Alvensleben, U., Anders, P. New Astron. 8, 155-171 (2003).
13. Strader, J., Brodie, J., Forbes, D., Beasley, M., Huchra, J. Astron. J. 125, 1291-1297 (2003).
14. Brodie, J., Huchra, J. Astrophys. J.. 379, 157-167 (1991).
15. Girardi, L. et al. Astron & Astrophys. 391, 195-212 (2002).
16. Sersic, J. Obs. Astron. Univ. Nac. Cordoba. (1968).
17. Graham, A., Erwin, P., Caon, N. , Trujillo, I. Astro-ph/0206248. (2003).
18. Binggeli, B., Jerjen, H. Astron & Astrophys. 333, 17-26 (1998).
19. Graham, A., Guzman, R., Astron. J. 125, 2936-2950 (2003).
20. Grebel, E., Gallagher, J., Harbeck, D. Astron. J. 125, 1926-1939 (2003).
21. Choi, P., Guhathakurta, P., Johnston, K. Astron. J. 124, 310-331 (2002).
22. Hodge, P. Astrophys. J.. 182, 671-696 (1973).
23. Grebel, E. IAU Symposium 192 The stellar content of local group galaxies, ed. P. Whitelock & R. Cannon, San Francisco ASP, (1999).
24. Ghigna, S., Governato, F., Lake, G., Quinn, T., Stadel, J. Mon. Not. R. Astron. Soc. 300, 146-162 (1998).
25. Harris, W., Harris, G. Astron. J. 122, 3065-3069 (2001).
26. This paper is based on observations made with the Hubble Space Telescope which is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract No. NAS5-26555, and the W. M. Keck Observatory, which is operated jointly by the California Institute of Technology and the University of California. The authors wish to thank the ACS team, R. de Grijs for supplying the combined ACS images, and A. Russell for her comments on the text.