text:
# Ferromagnetic materials in the zinc-blende structure

## I Introduction

For magnetoelectronic devices, highly spin polarized ferromagnets have a large potential in applications prinz . Half metallic ferromagnets (HMF) are the most desirable materials in this respect. These materials have one spin channel that is metallic, while the other spin channel has a gap at the Fermi level. Hence, 100 $`\%`$ spin polarization of the conducting electrons is expected from these materials. Based on band-structure calculations, de Groot et al. degroot  first predicted half metallic behavior of C1b type Heusler alloys of NiMnSb and PtMnSb. Since then, there have been predictions of many half metallic materials. A practical aspect of these materials is the recent breakthrough utilizing the spin of electrons in existing semiconductor electronic devices. To manipulate spins of electrons inside the semiconductors, one first needs to inject spins inside the semiconductors. Spin injection can be done either from ferromagnetic metals or from diluted magnetic semiconductors (DMS) doped with magnetic ions. Ohno et al. holeinject  successfully injected spin polarized holes from a p-type III-V DMS to a semiconductor. From theoretical calculations of magnetotransport, Egues egues  proposed a spin filter for spin-dependent optoelectronics applications. Later on, Fiederling et al. fiederling  successfully injected spin polarized electrons from an n-type II-VI semiconductor. However, according to Schmidt et al. schmidt , the spin injection from a metal to a semiconductor is practically impossible due to the huge conductivity mismatch between a metal and a semiconductor. The use of a DMS is instead rather advantageous. Mn doped GaAs is a prototype of this class of systems but the highest available Curie temperature is only 110 K ohno  due to the limited solubility of Mn in GaAs. The use of other half metallic systems e.g. CrO2 and La0.7Sr0.3MnO3 is also restricted due to lower Curie temperature. Very recently, there have been reports on room temperature ferromagnetism in Mn doped GaP compounds mngap . Hence, the quest for new half metallic materials with room temperature ferromagnetism continues.

Recently, there have been efforts to grow binary compounds in the zinc-blende structure using molecular beam epitaxy. One of the advantages of these materials is that the zinc-blende structure is maintained when grown on a zinc-blende semiconductor substrate. Akinaga et al. akinaga  grew CrAs on GaAs substrates and found CrAs to be a ferromagnet having a Curie temperature higher than 400 K. The experimental results were complemented by ab-initio electronic structure calculations showing that the magnetic moment was 3 $`\mu _B`$. Also CrSb was grown in the zinc-blende structure crsb . These binary compounds do not exist naturally in the zinc-blende structure but the non equilibrium growth by MBE has facilitated the creation of these compounds with novel magnetic properties.

Sanvito and Hill sanvito  examined the ground state properties of MnAs in a NiAs and a zinc-blende structures by first principles pseudopotential calculations. Shirai shirai  performed electronic structure calculations for the zinc-blende compounds of transition metal (TM) and As. The early TM-As compounds were shown to be ferromagnets and the later compounds were antiferromagnets. Also, electronic structure calculations have been performed for different crystal structures of MnP, MnAs and MnSb by Continenza et al. conti . Galanakis galan  examined the half metallicity of the bulk and surface of zinc-blende CrAs in various lattice constants by electronic structure calculations. A full potential screened Korringa-Kohn-Rostoker method was used. Galanakis concluded that a Cr terminated CrAs (001) surface remains half metallic both for the experimental lattice constants of GaAs and InAs. On the contrary, an As terminated surface loses its half metallicity due to the presence of the surface dangling bond states.

The motivation of this work is to design new materials in a crystal structure that is suitable for spin electronics applications and provide a guideline for experimental work that attempt to grow these materials. The paper is organized as follows : In section ’Computational details’, we describe in detail the computational part. In the results section, a subsection ’Electronic structure and magnetism’ describes the electronic structure and magnetism of the binary compounds considered followed by the discussions about the interatomic exchange interactions and critical temperatures in some of the systems studied. Finally, structural stabilities of MnAs and CrSb compounds in various structural phases are presented.

## II Computational details

A zinc-blende structure has two fcc sublattices displaced from each other along the body diagonal by a vector ($`\frac{1}{4}\frac{1}{4}\frac{1}{4}`$). For calculations in the ferromagnetic phase, the unit cell used in the calculations is the same as in the fcc lattice and has 1 TM at the (000) and the group V atom at the ($`\frac{1}{4}\frac{1}{4}\frac{1}{4}`$) position. Calculations have been performed by an ab-initio plane wave pseudopotential code (VASP) vasp . Pseudopotentials generated using the Projector augmented wave (PAW) paw  method were used throughout the calculations. These pseudopotentials are suitable for transition metals (TM) and the calculations yield accurate results tmgaas . A kinetic energy cut-off of 600 eV was used for the plane waves included in the basis set. Perdew-Wang GGA (generalized gradient approximation) exchange -correlation pw  was used. According to the calculations of Continenza et al. conti , GGA is essential for obtaining accurate equilibrium volume and magnetic moments. A 20x20x20 k-points grid was used in the Monkhorst Pack scheme mp  which yielded 770 k points in the irreducible Brillouin zone and subsequently used for the Brillouin zone integration within the tetrahedron method. Local properties such as local density of states and local magnetic moments were calculated by projecting the wave functions into spherical harmonics as described in ref. eichler . The values of the experimental lattice constants used in these calculations were taken from Ashcroft and Mermin ash . Calculations were done for the TM-V compounds for the lattice constants of the III-V semiconductors having the same Group V element. In summary, we considered lattice constants of GaAs, GaP, GaSb, InAs, InP and InSb compounds which covers a wide range (5.45 Å for GaP to 6.48 Å for InSb).

The pseudopotential calculations have been complemented footnote  with calculations based on non-collinear tight-binding linear muffin-tin orbital method (TB-LMTO-ASA)andersen  including combined correction terms in the one-electron Hamiltonianskriver . The generalized gradient approximation (GGA) were used for the exchange-correlation potentialpbe . Equal Wigner-Seitz sphere radii were used for all atoms, as well as for the empty spheres representing the two types of interstitial sites in the zinc-blende structure. A dense 27x27x27 k-point grid was used in the Brillouin zone integration and to speed up convergence each energy eigenvalue was associated with a Gaussian of width 85 meV. Spin wave excitations have been studied using the frozen-magnon approachrosengaard ; halilov ; pajda . Then the total energy $`E(𝐪,\theta )`$ is calculated for spiral magnetic structures where each magnetic moment is defined by the Euler angles

$$\theta _i=const;\varphi _i=\text{q}\text{R}_i,$$
(1)
where $`\theta _i`$ and $`\varphi _i`$ are the polar and azimuthal angles, $`\text{R}_i`$ the position of ion $`i`$ and q is the spiral propagation vector. Spin spirals possess a generalized translational symmetrysandratskii , which allow for calculations in the chemical unit cell without need for a supercell. It can be shown, that the true spin wave excitation energies are related to the spiral energies through the relation

$$\omega (\text{q})=\frac{4}{M}\frac{\mathrm{\Delta }E(\text{q},\theta )}{sin^2\theta },$$
(2)
where $`\mathrm{\Delta }E(\text{q},\theta )=E(\text{q},\theta )E(0,\theta )`$. The last equation is valid for one atom per cell. Thus the formula is correct if one consider the magnetic response from the transition metal only and neglect the small induced magnetic moment on the group V element. In our calculations, we use an angle $`\theta =20^{}`$ and the spiral energies were obtained by employing the magnetic force theorem. We tested other angles as well as self-consistent calculations and found that the differences in total energy is within a few percent.

In order to obtain thermodynamic properties, the spin wave excitation energies were mapped to an effective Heisenberg Hamiltonian of classical spins

$$H=\underset{ij}{}J_{ij}\text{e}_i\text{e}_j,$$
(3)
where $`J_{ij}`$ are exchange interactions between the transition metal atoms. The critical temperatures $`T_c`$ were estimated both from mean field theoryash  (MFA) and Monte Carlo simulations (MC). For the MC simulations, we used the standard single flip Metropolis algorithm and the critical temperature were obtained from the ’cumulant crossing method’ Binder . In this method, the reduced fourth order cumulant of the order parameter (magnetization) $`U_L`$, defined as

$$U_L=1\frac{M^4}{3M^2^2},$$
(4)
is calculated for different sizes of the simulation box. The curves of $`U_L`$ have a common intersection point at a fixed point $`U^{}`$, which corresponds to the critical point. Hence, a value of $`T_c`$ is obtained from the intersection point of $`U_L`$ for different lattice sizes.

## III Results and discussions

### III.1 Electronic structure and magnetism

The zinc-blende lattice introduces a tetrahedral crystal field under which the transition metal atomic $`d`$-states are split into triply degenerate $`t_2`$ and doubly degenerate $`e`$ states. $`e`$ states lie lower in energy than $`t_2`$ states in a tetrahedral crystal field in contrast to the octahedral crystal field. The $`t_2`$ states hybridize with the anion $`p`$-states. As a result, bonding and antibonding states appear. Also, depending on the exchange splitting, the bonding and antibonding states may lie in different positions for different spin channels. $`e`$ states are of non-bonding character and do not take part in the hybridization.

In Fig. 1(a-f), we present the spin resolved DOSs of V, Cr and Mn compounds having the lattice constants of GaAs, InAs, GaSb, InSb, GaP and InP. In figs. (a), (b), (c), (d) and (f), all the compounds show half metallic behavior since there is no state at the Fermi level for the spin down DOS. The only exception is the case depicted in fig. 1(e) which doesn’t exhibit half metallic ferromagnetism as the Fermi level cuts both the spin bands. For the lower lattice constant of GaP, the $`pd`$ overlap and hybridization of the spin down states broadens the band to cut the Fermi level. Also, as the width of the gap in the spin down DOS is quite large, the Fermi level stays within the gap with the variation in the lattice constant. It is worthy to be mentioned here that these compounds can be regarded to be the derivatives of diluted magnetic semiconductors, e.g. MnxGa1-xAs for a Mn doped GaAs system having $`x`$=1. Several first principles calculations sanprb01 ; schultess ; zhao ; pashaprl  confirmed that the total spin magnetic moment of Mn doped GaAs for various concentrations of Mn is 4 $`\mu _B`$/cell with a half metallic nature. It was also shown tmgaas  that Cr doped GaAs is a half metallic ferromagnet with a total spin magnetic moment of 3 $`\mu _B`$/cell. We find here a similar trend for the extreme limit at a high TM concentration. A systematic variation in the DOSs of V, Cr and Mn is seen. As the number of valence electrons increases by going from V to Mn, the main features of the DOSs shift to lower energies relative to the Fermi level. This can be seen from figs. 1 (b) and 1 (f). The exchange splitting is maximum for Mn and here one finds a maximum separation between the spin up and spin down bands.

A close inspection of Fig. 1(a) reveals that VAs behaves like a magnetic semiconductor for the lattice constant of GaAs, i.e. the Fermi level passes through a gap for both spin up and spin down channels. We varied the lattice constant of VAs around the GaAs lattice constant and found that the small gap is sensitive to the lattice spacing. Only in a very small region around the GaAs lattice constant, do we find that the gap persists. This is an interesting result since it is one of quite few ferromagnetic semiconductors found in a stoichiometric composition. In addition, it is to be noted that unlike CrAs, MnAs doesn’t show half metallic behavior for the GaAs lattice constant. It is instead half metallic for the InAs lattice constant (see inset of Fig. 6(b)).

In Table I, we present the total magnetic moment/cell for all the compounds considered. The half metallic ferromagnetic moments are marked in bold. Half metallicity may be monitored from the integer moments of the unit cell. It is seen that from all the combinations, VAs and CrAs are HMFs for both GaAs and InAs lattice constants. VSb and CrSb are also half metallic for GaSb and InSb lattice constants. VP, CrP and MnP are only HMFs at the InP lattice constant. This shows that the lattice constants and the resulting electronic structures play an important role in determining half metallicity, a finding that is also evident from Fig. 1. The moments listed in Table I are consistent with the calculated moments of Cr and Mn impurities in a GaAs host, if one neglects the effect of other defects e.g. As antisites pashaprl . To be more specific we note that Mn as an impurity in GaAs carries a magnetic moment of 4 $`\mu _B`$/Mn atom in all theoretical first principles calculations based on local density approximation (LDA) that neglect other lattice defects such as As antisites. If one considers As antisite defects, or Mn interstitial atoms that couple antiferromagnetically to the Ga substituted Mn atoms, the effective moment becomes $``$2.4 $`\mu _B`$/Mn atom in agreement with experimental data pashaprl . For the currently investigated materials only a limited set of experimental data have been published. However, CrAs has been synthesized in the zinc-blende structure and the magnetic moment was measured to be 3 $`\mu _B`$/Cr atomakinaga  in agreement with the present theory and previous calculationsakinaga . Since Mn has one extra electron in the spin up band, it is natural that the moment of MnAs is 4$`\mu _B`$/Mn atom, which is in agreement with the moment of dilute Mn in GaAs (if one neglects defects).

The explanation of the half metallicity and the integer moments is similar to that of a III-V diluted magnetic semiconductor. In a III-V DMS, a cation is removed by creating a vacancy containing 3 holes in the valence band. When a transition metal occupies this cation site, it donates 3 electrons to the anion dangling bonds to saturate the bonding. The remaining unpaired d electrons of the TMs are responsible for the magnetic moments. Due to the large exchange splitting according to Hund’s rule, the spin up and spin down d bands are separated from each other having different degree of hybridization with the valence bands. Hence, V, Cr and Mn have magnetic moments of 2, 3 and 4 $`\mu _B`$, respectively from the unpaired spin up d electrons. However, the widths of the $`d`$-bands are fairly large indicating a substantial amount of hybridization with the anion-$`p`$ orbitals.

TABLE I. Total magnetic moment/cell (in $`\mu _B`$). Each entry corresponds to the binary alloy consisting of a transition metal from the leftmost column and the group V element in the subsequent columns. The bulk semiconductors whose lattice constants are considered for each calculation are shown in suffices.

|  | AsGaAs | SbGaSb | PGaP | AsInAs | SbInSb | PInP |
| --- | --- | --- | --- | --- | --- | --- |
| Ti | 0.00 | 0.00 | 0.01 | 0.29 | 0.00 | 0.04 |
| V | 2.00 | 2.00 | 1.90 | 2.00 | 2.00 | 2.00 |
| Cr | 3.00 | 3.00 | 2.74 | 3.00 | 3.00 | 3.00 |
| Mn | 3.53 | 3.89 | 3.19 | 4.00 | 4.02 | 4.00 |
| Fe | 2.36 | 2.69 | 1.06 | 3.20 | 3.26 | 3.15 |
| Co | 0.02 | 0.23 | 0.00 | 1.15 | 1.39 | 0.94 |
| Ni | 0.00 | 0.02 | 0.01 | 0.04 | 0.02 | 0.11 |

We also investigated compounds of Ti, Fe, Co and Ni, and found that they do not show half metallic ferromagnetism (see Table I). Fe compounds have fairly big moments in the ferromagnetic phase except for FeP in the GaP lattice constant. CoSb in the InSb lattice constant has the highest moment among the Co based compounds. Ti and Ni based compounds do not have prominent ferromagnetism. As only V, Cr and Mn based compounds show interesting and useful magnetic behavior, we will now concentrate on those compounds only.

TABLE II. Total energy differences in eV/Mn atom ($`\mathrm{\Delta }E=E_{AFM}E_{FM}`$) where $`FM`$ and $`AFM`$ are the ferromagnetic and antiferromagnetic arrangements between Mn atoms respectively. The symbol ’-’ is for the compositions where calculations were not performed.

|  | GaAs | InP | InAs | GaSb | InSb |
| --- | --- | --- | --- | --- | --- |
| VAs | 0.17 | - | 0.19 | - | - |
| CrAs | 0.27 | - | 0.32 | - | - |
| VP | - | 0.2 | - | - | - |
| CrP | - | 0.31 | - | - | - |
| MnP | - | 0.14 | - | - | - |
| VSb | - | - | - | 0.18 | 0.12 |
| CrSb | - | - | - | 0.31 | 0.27 |

### III.2 Interatomic exchange interactions

In this section, we compare the total energies of the HMFs found from our calculations, in the ferromagnetic (FM) and antiferromagnetic (AFM) phases. For briefness, calculations were only made for selected V, Cr and Mn compounds (see Table II). In our calculations, we considered an 8 atom (TM)4X4 cell in a layer geometry in the (001) direction with 2 TM (TMI) atoms in the plane at z=0.0 and 2 TM atoms at z=a/2 (TMII), where a is the lattice parameter. Here TM and X indicate transition metal (V, Cr and Mn) and the group V element (As, P and Sb) respectively. TMI and TMII atoms can couple ferromagnetically or antiferromagnetically.

The energy differences are given in Table II and are also plotted in Fig. 2. The energy differences, which reflects the interatomic exchange coupling strengths are fairly large and are found to always stabilize ferromagnetism. From Fig. 2 and Table II, it is clear that the strongest exchange interaction is found for the Cr based compound, where the energy difference is quite appreciable, indicating a high Curie temperature. An inspection of the DOSs for the FM and the AFM structures reveal that for the AFM structure, the DOS at the Fermi energy is non zero for both spin up and spin down bands whereas in the FM alignment, only spin up states exist at the Fermi level. This stabilizes the FM alignment, since the eigenvalue-sum in general lowers the total energy when electron states are removed from the Fermi energy.

In this section, we also present results of the calculations of critical temperatures using spin wave spectra and statistical mechanics. We have calculated $`T_c`$ for VAs, CrAs and MnAs for lattice constants of GaAs and InAs by extracting exchange interactions from calculated spin wave spectra, which subsequently were used in a classical Heisenberg Hamiltonian.

In Fig. 3, the calculated spin wave spectra for the systems are shown. MnAs in the lattice constant of GaAs was found to be antiferromagnetic and is therefore not shown in the figure. The other systems were found to be ferromagnetic and half metallic with total magnetic moments of 2.0, 3.0 and 4.0 $`\mu _B`$, respectively for VAs, CrAs and MnAs. In Table.III, calculated values of $`T_c`$ are shown. In the case of CrAs, a very high $`T_c`$ around 1000 K is estimated which could be very promising from an application point of view. $`T_c^{MC}`$ is smaller than $`T_c^{MFA}`$ due to the fact that fluctuations are absent in the latter. Normally, MC simulation gives $`T_c`$ values closer to experimental values than those obtained from mean field theoryrosengaard .

TABLE III. Calculated critical temperature of VAs, CrAs and MnAs from the MFA and MC simulations. Values at the left correspond to the lattice constant of GaAs and the right correspond to InAs.

As is clear from Table III, the critical temperatures of VAs and CrAs are in general larger for the substrate with the smallest lattice constant (GaAs). Hence it is clear that the critical temperatures of VAs and CrAs follow the same trend as that of rare-earth systems, where in general a reduced lattice constant results in enhanced critical temperatures (as for instance measured by an applied external pressure rare ). In rare-earths this phenomenon is easily understood from the RKKY interaction, which is known to be responsible for the magnetic ordering in these elements. By analogy we can say that the presently studied materials (VAs and CrAs) have rare-earth like characteristics, in the sense that the exchange interaction between localized moments is mediated by more diffuse valence states.

In fig. 4, we show the exchange interactions ($`J_{ij}`$) in CrAs as a function of Cr-Cr distance. It is clear that the first nearest neighbor exchange interaction is the dominant one and is ferromagnetic in nature. $`J_{ij}`$ decreases rapidly with increasing distance between Cr atoms. Further neighbor exchange interactions are antiferromagnetic in nature and are very weak. Similar kind of observation was reported by Sandratskii and Bruno sand  in the case of Mn doped GaAs.

## IV Band structure and calculation of reflectance

In Fig. 5(a) and (d), we show the band structures of MnAs in the InAs lattice constant for two symmetry directions viz. $`\mathrm{\Gamma }`$-L and $`\mathrm{\Gamma }`$-X directions respectively. In (b) and (c), we show the same for bulk InAs calculated for the experimental lattice parameter. Our calculation does not produce a gap at the $`\mathrm{\Gamma }`$ point, so the top of the valence band and the bottom of the conduction band merge with each other. This is a well known limitation of using LDA or GGA within density functional theory. For MnAs, the spin up bands cross the Fermi level for the two symmetry directions shown the figure and also for other directions (omitted in the figure). There is gap in the spin down channel. As described above, the $`e`$ bands are dispersionless as they do not mix with others. The antibonding $`t_2`$ bands show large dispersion as they are hybridized strongly with the As $`p`$ bands.

We now give a rough estimate of the reflectance and transmittance from a ferromagnetic metal/semiconductor interface from the band matching of the individual species. We use a simple free electron model, that was used also by Kilian and Victora kilian2000 . For an electron traveling through a one-dimensional potential step, the reflectance is given by

$$R=\frac{(k_{fm}k_{sc})^2}{(k_{fm}+k_{sc})^2},$$
(5)
where $`k_{fm}`$ and $`k_{sc}`$ are the wave vectors in the ferromagnetic metal and semiconductor respectively. As the spin up bands are the only states available at the Fermi energy, we consider the following for the spin up electrons only. For InAs, as the conduction band minimum is at the $`\mathrm{\Gamma }`$ point, the lowest energy levels available to accept spins from the metal are at the $`\mathrm{\Gamma }`$ point. We calculate $`R[100]`$, $`R[110]`$ and $`R[111]`$ for  ($`\mathrm{\Gamma }`$-$`X`$), ($`\mathrm{\Gamma }`$-$`K`$) and ($`\mathrm{\Gamma }`$-$`L`$) directions, respectively. It is known from experiments that for e.g. an interface of Au/n-InAs, the Fermi level is pinned 0.1 eV above the conduction band minimum. We use this information in our model assuming that the properties of the metal do not alter strongly the band alignments. Also, the spins within ($`k_BT0.025eV`$) of the MnAs Fermi level are the only ones available for transport properties. From this information we can calculate $`k_{sc}`$ for the three different directions. $`k_{fm}`$ for the three different directions are available from the band structures of MnAs directly. Then we calculated the reflectance from equation 1. The calculated values are $`R[100]`$=0.73, $`R[110]`$ = 0.89 and $`R[111]`$ = 0.71. This means that the transmissions in these directions will be rather small. However, one should bear in mind that in reality the interface is not ideal. There are interface states and different scattering mechanisms which can further decrease or increase the transmission through the interface. A full transport calculation taking into account a proper model of the interface could give a more realistic estimate, and our model should be viewed as an estimate only.

## V Structural stability

In this section, we present the relative structural stabilities of these compounds. For briefness, we choose to illustrate our conclusion only for 2 out of the 9 HMFs shown in Table I, i.e. CrSb and MnAs. The conclusions reached for these two compounds should be valid for all the HMFs listed in Table I. We compare the total energies of two crystal structures, hexagonal NiAs (that is the low temperature structure) and the zinc-blende (ZB) structure.

From Fig. 6, it is evident that for both MnAs and CrSb, the NiAs structure is more stable than the ZB structure. Similar results were obtained for MnAs by Sanvito and Hill sanvito . This is in accordance with naturally occurring compounds of CrSb and MnAs in the NiAs structure. However the total energy curves for the NiAs structure crosses the total energy curves of ZB structures for certain volumes. This occurs for a volume of 46.3 $`\AA ^3`$/f.u. in MnAs and 59.6 $`\AA ^3`$/f.u. in CrSb. Incidentally, these volumes are quite close to that of GaAs in case of ZB-MnAs and that of GaSb in case of ZB-CrSb compounds. For larger volumes, the ZB structure is more stable than the NiAs structure. Also, these compounds become half metallic ferromagnets according to the calculations. Therefore, ZB heterostructures of suitable combinations should be possible to use to obtain half metallic materials with a large critical temperature. Hopefully this fact can be a useful guideline for experimental work in trying to grow heterostructures for spin-electronics applications. Finally, the variation of magnetic moment of MnAs with the variation in lattice constant is shown in the inset of Fig. 6(b) where it is clear that it becomes half metallic only near the vicinity of a lattice constant of 6 Å .

The results presented above are in excellent agreement with the existing ones, wherever available. However, there are certain limitations of our approach. One of them is the underestimation of band gaps in semiconductors within the use of LDA or GGA. A more appropriate way is to use GW approximation gw  for these materials which corrects the band gaps extremely well. To our knowledge, no theoretical calculation exists for diluted magnetic semiconductors using this more complicated theory. The other issue is the effects of native defects which are always present in the growth processes. One should also take into account these effects to model a more realistic situation. But we believe that the results presented in this paper are quite significant as far as trends are concerned and also a guideline for future experiments.

## VI Conclusion

In conclusion, we have identified a number of possible half metallic ferromagnets in the zinc-blende structure. Ferromagnetism has been shown to be more stable than an antiferromagnetic structure and exchange couplings are calculated to be sufficiently strong to produce high Curie temperatures, especially for Cr based compounds. We also give a parameter range in substrate lattice constants where MnAs and CrSb are expected to grow and to have half metallic properties. Finally, we predict that VAs grown at suitable conditions should be a ferromagnetic semiconductor with a large critical temperature. We hope to motivate experimental work in trying to grow such materials for future applications. From a theoretical point of view, the future scope of this work is to investigate the interface properties of these binary compounds and the zinc-blende semiconductors e.g. CrAs and GaAs. A detailed investigation including structural relaxations of these interfaces is in progress.

###### Acknowledgement 1

This work was done during the stay of one of the authors (B.S.) in Max-lab, Lund, Sweden. Support from the Swedish Science Foundation (VR), Swedish Foundation for Strategic Research (SSF) and the Göran Gustafsson foundation are acknowledged. Also, we acknowledge support from National Supercomputer Center (NSC).


format: text

text:
# Topological constraints on spiral wave dynamics in spherical geometries with inhomogeneous excitability

## I Introduction

Geometry and inhomogeneity influence pattern formation in chemical and biological systems. Murray (1989); Kapral and Showalter (1996) One example where these two factors play a crucial role is in the experimental observations of distinctive spiral wave dynamics on the surfaces of spherical beads, which are excitable inhomogeneous chemical media Maselko and Showalter (1989, 1990). A biological example is the origin of abnormal cardiac rhythms in the human heart which depend on the anatomical substrate. The heart possesses a complex nonplanar geometry with multiple chambers, with holes corresponding to valves and blood vessels. Some serious arrhythmias are associated with circulating spiral waves similar to those observed in chemical media Winfree (2001). Since an abnormal anatomical substrate is a common finding in patients with some types of cardiac arrhythmias, and interventions that modify the anatomy are an accepted form of therapy Stein et al. (2002), theoretical analyses of the interplay between geometry of the substrate and dynamics may help in the therapy of cardiac arrhythmias.

In this paper we study spiral wave dynamics on (punctured) spheres with spatially inhomogeneous excitability. We show for punctured spheres that the sum of indices which characterize each spiral has to be zero. Moreover, we demonstrate that topological constraints imposed by the spherical geometry and inhomogeneity in excitability lead to the formation of pairs of spirals, with distinctive transient dynamics or as stable attractors, in which one spiral acts as a source and a second as a sink leading to a source-sink pair under a broad range of conditions. Our results explain the experimental observations of spirals on spherical beads Maselko and Showalter (1989, 1990). While we do not consider detailed models of cardiac wave propagation, our results may apply to some generic aspects of atrial arrhythmias because the thin walls of the atria can be described as two-dimensional (2d) inhomogeneous excitable media with specific geometrical features. num

## II An Index Theorem for Phase Singularities

The mathematical description of spiral waves is based on the notion of phase which in turn allows one to characterize spiral waves by an index. From this description, a number of topological results placing restrictions on spiral wave dynamics can be derived. Glass (1977); Winfree and Strogatz (1983); Sumners (1987); Winfree (2001); Cruz-White (2003)

With the exception of a finite number of singular points, with each point in an orientable and compact two-dimensional differentiable manifold $`M`$ we identify a unique phase lying on the unit circle, $`\mathrm{\Phi }S^1`$. The resulting phase map or phase field is assumed to be continuously differentiable, except at the singular points. The manifold can be triangulated Hocking (1988) (subdividing it into a set of polygons), where none of the edges or vertices of the polygons pass through a singularity. The index, $`I`$ (sometimes also called the topological charge or winding number) of a curve $`C`$ bounding a polygon is found by computing the line integral

$$2\pi I=_C\mathrm{\Phi }d𝐥,$$
(1)
where the polygon is always traversed in a clockwise orientation. By continuity of $`\mathrm{\Phi }`$, $`I`$ must be an integer. The index of a singular point is uniquely defined as the index of any curve $`C`$ provided that $`C`$ encircles the point but no other singular points. The index of a curve that does not enclose any singular points is obviously zero.

If the manifold $`M`$ has no boundaries, each edge of the triangulation is an edge of two polygons. Since the line integral adds up the change in phase along the various edges of the polygon, the sum of the indices of the singular points for a phase field in $`M`$ is

$$\mathrm{\Sigma }_MI=0,$$
(2)
where the sum is over all the singular points. This follows since the contribution of the change in phase of each edge to the total integral is counted twice, but since the edge is traversed in opposite directions each time, the net contribution of each edge is zero. This index theorem is applicable to tori, and other surfaces of genus different from 0. However, unlike the more familiar Poincaré index theorem (see Ref. Guillemin and Pollack (1974) p. 74 for vector fields) the sum of the indices of the singular points does not depend on the genus of the surface.

This index theorem for manifolds without boundaries can be extended to manifolds with boundaries. In the following, we will consider the case of structures that arise from puncturing orientable and compact two-dimensional differentiable manifolds. The index of a hole can be uniquely defined as the index of a curve $`C`$ provided that $`C`$ encircles the hole but no other holes or singular points and $`C`$ is positively oriented with respect to the domain which contains the hole and is bounded by $`C`$. Applying this definition and taking the summation in Eq. (2) over the singularities and the hole, or the holes if there is more than one hole, the index theorem can be proven by the same line of arguments as for the case of manifolds without boundaries.

This extension is important in the heart, for example, where the atrium is punctured by valves and veins. In such cases one is led to consider manifolds with holes; for example, a sphere with a hole. A sphere with a hole is topologically equivalent to a disk, and, indeed, the results for disks and for spheres with holes are consistent: For the disk $`D^2`$, bounded by a curve $`C`$,

$$\mathrm{\Sigma }_{D^2}I=_C\mathrm{\Phi }d𝐥,$$
(3)
so that the sum of the indices of the singular points in the disk is equal to the index of the curve $`C`$ bounding the disk. If there is a single singular point on the disk, with an index of +1, the index of the curve bounding the disk will also be +1. Imagine now the boundary of the disk to be brought together (like a draw-string bag) so that the boundary of the disk now defines a hole in the sphere. In this geometry the curve $`C`$ will be traversed in an opposite orientation (the hole is now inside $`C`$) from the direction it was traversed when it was the boundary of the disk. Now if there is a singular point with an index of +1 on the sphere, the index of the hole is -1, so that the sum of the indices is again zero.

Since it is necessary to conserve the sum of the indices, singularities of index $`\pm 1`$ usually arise and are destroyed in pairs of opposite sign Winfree and Strogatz (1983). An exception occurs when singularities are destroyed by collision into a boundary, so that the index of the singular point and the index of the boundary both change simultaneously. Another exception occurs if there are singularities with index different from one. In such cases interactions between different singularities can lead to destruction or creation of odd numbers of singularities Zaritski and Pertsov (2002).

## III The FitzHugh-Nagumo Equation

The FitzHugh-Nagumo (FHN) equation, num

$`{\displaystyle \frac{u}{t}}`$ $`=`$ $`ϵ^1\left({\displaystyle \frac{u^3}{3}}+uv\right)+D_u^2u,`$
$`{\displaystyle \frac{v}{t}}`$ $`=`$ $`ϵ\left(u\alpha v+\beta \right)+D_v^2v,`$ (4)
where $`u(𝐫,t)`$ and $`v(𝐫,t)`$ are two scalar fields, $`ϵ^2`$ is the ratio of the time scales associated with the two fields, and $`D_u`$ and $`D_v`$ are the constant diffusion coefficients, is a prototypical model for excitable media. We choose $`D_u=2`$ and $`D_v=0`$. The real parameters $`\alpha `$ and $`\beta `$ characterize the local dynamics and, hence, the local excitability. Whenever $`0<\alpha <1`$, $`\alpha ϵ^2<1`$ and $`|\beta |>\beta _H(1\alpha ϵ^2)^{1/2}(1/3(2\alpha +\alpha ^2ϵ^2)1)`$, the FHN system is excitable. At $`\beta _H`$ a Hopf bifurcation occurs such that for $`|\beta |<\beta _H`$ the system exhibits oscillations. In the following we take $`\alpha =0.2`$, $`ϵ=0.2`$ and $`\beta >\beta _H=0.863\mathrm{}`$.

We consider a spherical shell whose outer and inner radii are $`R_e`$ and $`R_i`$, respectively, and focus on thin spherical shells where $`R_i=40`$, $`R_e=42`$. The radii are large enough to avoid a curvature-dependent loss of excitability Davydov et al. (2002), and the shell is sufficiently thin that the dynamics is effectively 2d corresponding to the dynamics on a sphere. det  The initial condition is a domain of “excited” state, adjacent to a domain of the “refractory” state. Both domains have the forms of slices of the same size oriented from the north to the south pole det  and yield a pair of counter-rotating spirals.

In order to apply the topological results to the FHN equation, it is necessary to first define the phase. We define a phase, $`\mathrm{\Phi }(𝐫,t)`$, based on the equation $`\mathrm{tan}\mathrm{\Phi }(𝐫,t)=v(𝐫,t)/u(𝐫,t)`$ if $`v(𝐫,t)0`$ and $`u(𝐫,t)0`$. Thus, singular points at given $`t`$ are points $`𝐫`$ in the medium for which $`v(𝐫,t)=0`$ and $`u(𝐫,t)=0`$. For each $`t`$, we obtain a continuously differentiable phase map $`^t=\mathrm{\Phi }(:,t)|_{𝒟^t}`$ that associates to each point in a well-defined domain $`𝒟^t`$ a phase lying on the unit circle, $`\mathrm{\Phi }S^1`$. For our FHN medium, the domain is the surface of a (punctured) sphere reduced by a finite number of points where the phase is singular at fixed $`t`$.

Rotating spiral waves in the FHN equations are obviously associated with a singular point which is called the spiral core. In what follows, we assume that there are only single-armed spirals so that a clockwise rotating source has an index of +1, and a counterclockwise rotating source has index -1. A clockwise rotating sink has an index of -1, and a counterclockwise rotating sink has index +1. From Eq. (2), it is impossible to have a single rotating spiral wave on a sphere; in addition, there must be at least another singular point or a hole with nonzero index.

For excitable media, a non-zero index of a hole implies that wave fronts travel permanently around the hole such that the numbers of fronts travelling clockwise and counterclockwise are different. This includes the particularly simple case of a single wave front travelling around the hole which can be considered as a spiral wave associated with the hole.

## IV Spiral Wave Dynamics in Spherical Geometries

### IV.1 Dynamics in excitability gradients

First consider homogeneous FHN media with a constant $`\beta =\beta _{ex}0.9`$. The wave fronts emanating from the spiral cores with opposite index “annihilate” along the equator such that each spiral determines the dynamics on half of the sphere (see Fig. 1, left panel) — similar to what has been observed in Ref. Yagisita et al. (1998) for a different excitable system.

We have shown that this behavior is robust with respect to disorder in excitability with small amplitude and correlation length. If random, uncorrelated spatial variations in $`\beta `$ exist on length scales much smaller than the diameter of the spiral core meander ran , the dynamics is able to average over such small-scale inhomogeneities. The robustness explains why such states have been experimentally observed in some chemical reactions on spherical beads which are intrinsically inhomogeneous Maselko and Showalter (1990).

Applying a constant gradient in $`\beta `$ as sketched in the right panel of Fig. 1 leads to a different scenario. The time evolution of the spiral pair may be partitioned into four distinct regimes as shown in Figs. 2 and 3.

Because of the gradient, the frequencies of the two spirals differ since a higher value of $`\beta `$ corresponds to lower excitability, which generally implies a lower spiral frequency Winfree (1991). During a short transient T, the spiral with the higher frequency assumes control of the dynamics Krinsky and Agladze (1983) on the sphere. The location on the sphere, where wave fronts emanating from the two spiral cores annihilate, moves toward the core of the low-frequency spiral. Finally, the low-frequency spiral core is pushed farther from the high-frequency spiral core Krinsky and Agladze (1983); Ermakova et al. (1960)(see Fig. 3). After this short transient, the wave fronts travel from pole to pole leading to the creation of a source-sink pair. This (intermediate) state is shown in Fig. 4 and corresponds to regime I in Figs. 2 and 3. Viewed from the low excitability end of the sphere, the waves wind into a small region about the core, reminiscent of the structure of antispirals, i.e., inward moving spirals seen in oscillatory media Gong and Christini (2003). However, the origin of this inward spiral motion in oscillatory media differs and is distinct from that observed here. In oscillatory media, either spirals or antispirals are stable depending on system parameters and the wave length diverges on the border in the parameter space between these two regimes. Thus, antispirals exist independently of spirals. This is not the case here because the generation of an inward moving spiral relies on the presence of a spiral source and spherical geometry. For example, consider the FHN system with a disc geometry and a radial gradient in excitability such that the maximum value of $`\beta `$ is located in the center of the disc. In this case, a source-sink pair cannot occur because the high-frequency spiral, acting as a source, would push the other (low-frequency) spiral out of the system, excluding the presence of any strong random inhomogeneities in excitability which may pin the low-frequency spiral and prevent its motion (see, e.g., Biktashev et al. (2001)). The lower panel of Fig. 4 shows that remnants of the low frequency spiral persist in a small area around the core. It has been speculated that antispiral waves might occur in cardiac tissue Vanag and Epstein (2001). While excitable media cannot support a regime of exclusive antispirals due to their excitable character, our results show that source-sink pairs with similar characteristics could form in the heart where the underlying dynamics is excitable, the medium is inhomogeneous and the topology is similar to a (punctured) sphere.

Spiral dynamics of the type described above has been observed by Maselko and Showalter Maselko and Showalter (1989) in experiments on the excitable Belouzov-Zhabotinsky reaction on spherical beads. They attributed the generation of spiral source-sink pairs to inhomogeneities in the medium related to differing chemical environments. This is consistent with our findings for systems with gradients and is further confirmed by the work reported in Ref. Yagisita et al. (1998). While the generation of source-sink pairs due to a gradient in the FHN medium investigated here is only an intermediate state (but one that persists for approximately 240 spiral periods in our simulations), random spatial variations of the excitability with a correlation length comparable to the diameter of the spiral core meander or larger can lead to a final state consisting of such a source-sink pair ran . This is due to the fact that the source can be trapped in a region of depressed local excitability.

Not only does the gradient in the FHN medium change the local excitability but it also induces a drift of the spiral cores Biktashev and Holden (1995). For our model, the drift is rather slow compared to the transition to the source-sink pair which takes place during the transient regime T. This can be seen in Fig. 2. (The fluctuations in $`\beta (𝐫(t))`$ are due to spiral meandering.) In regime I, the dominating spiral drifts toward lower excitability and its wave fronts continuously push the other core in the opposite direction, thus, keeping the distance $`d`$ between the cores approximately constant (see Fig. 3). The fluctuations in $`d(t)`$ are again due to meandering of the spiral cores. Because of this drift, the local excitabilities at the two spiral cores approach each other until they become equal.

At this point, regime II in Figs. 2 and 3 is entered. The dynamics change drastically: the slaved spiral reverses its drift direction and regains control over its own dynamics. Both spirals drift toward lower excitability. Due to the geometric constraints imposed by the spherical geometry, the spirals approach each other until they form a bound pair (at $`t7500`$).

They finally reach a stable state (at $`t9500`$) corresponding to regime III in Figs. 2 and 3. Neither the *average* distance between the spirals nor the *average* local excitability changes further. Yet, on top of the persisting unsynchronized meandering of the two spiral cores, the bound pair slowly moves along a (closed) equi-$`\beta `$ curve on the surface of the sphere. The direction of the motion depends on the initial condition, i.e., whether the spiral with positive index was initially closer to the point of lowest excitability or to the point of highest excitability than the spiral with negative index. The wave dynamics generated by this bound pair is shown in the inset of Fig. 2.

While kinematic theory applies only to spirals with large cores, it is instructive to note that this theory predicts that the direction of the drift due to gradients depends on the model system and its parameters Mikhailov et al. (1994). Although our simulations have shown the same direction of the drift for a range of parameters in the meander region of the FHN phase diagram, it is conceivable that, under different circumstances favoring a drift toward higher excitability, one spiral could act as a permanent source and a source-sink pair could be the final attractor. Such a scenario would also be consistent with the experimental findings in Ref. Maselko and Showalter (1989).

### IV.2 Punctured spheres

Next, we consider a homogeneously excitable sphere with a single hole. Two scenarios can be observed depending on the location of the hole with respect to the spiral pair. If a spiral wave is not permanently attached to the hole, the dynamics is very similar to the case without any hole. If one of the spirals is permanently attached to the hole, the frequency of this spiral is lowered. The size of the hole determines the frequency of the spiral because the wave front has to travel around the hole. The transient dynamics is similar to that in regime T for the case with a gradient; however, no drift of the spiral cores is induced and the final state is a spiral source-sink pair as shown in Fig. 5. Not only is the net index conserved during the transition to a spiral source-sink pair but so is the index of the individual spirals: during the transition an outgoing counterclockwise (clockwise) oriented spiral is converted into an ingoing clockwise (counterclockwise) oriented spiral. Thus, the formation of spiral source-sink pairs conforms with the topological constraints.

If a gradient as well as a hole is present, the spiral drift discussed earlier also determines the final state, which depends on the hole’s location in the gradient field. For instance, if the point of lowest excitability in the medium is on the hole boundary, simulations show that the spiral which is not attached to the hole will stabilize close to this point. In this final state odd numbers of wavefronts are attached to the hole, again conserving the topological charge of the hole.

## V Concluding Remarks

Inhomogeneities due to spatially-varying excitability on (punctured) spherical shells lead to complex spiral wave dynamics and the formation of source-sink spiral pairs in excitable media. The results presented here are immediately applicable to excitable media in more complicated geometries such as tori or multi-holed tori and to situations in which multi-armed spirals are found. This includes mathematical modelling of cardiac tissue. The approach taken in this paper stresses constraints and aspects that apply to, and must be observed in, all realistic models of the heart satisfying certain criteria of continuity. There are also implications for the treatment of cardiac arrhythmias. In cardiology it is sometimes possible to develop maps showing the timing of the excitation over limited regions of heart Stein et al. (2002). In this case, a sink might be confused for a source (of the arrhythmia), and this might have implications for the diagnosis of the mechanism and the choice of therapy. The current work shows how partial knowledge about what is happening in some regions that could be observed, might be helpful in establishing properties of dynamics that could not be observed. While the types of sinks we have described here have only been observed in chemical media Maselko and Showalter (1989, 1990) so far, we certainly expect their existence in the cardiological domain.

We thank F. Chavez, M. Bär, S. G. Whittington and D. Sumners for helpful discussions and G. Rousseau for providing numerical tools. This work was supported in part by a grant from MITACS.


text:
# Critical Magnetic Prandtl Number for Small-Scale Dynamo

Present-time address: \]DAMTP/CMS, University of Cambridge, Wilberforce Rd., Cambridge CB3 0WA, UK.

## Abstract

We report a series of numerical simulations showing that the critical magnetic Reynolds number $`\text{Rm}_\text{c}`$ for the nonhelical small-scale dynamo depends on the Reynolds number Re. Namely, the dynamo is shut down if the magnetic Prandtl number $`\text{Pr}_\text{m}=\text{Rm}/\text{Re}`$ is less than some critical value $`\text{Pr}_{\text{m,c}}1`$ even for Rm for which dynamo exists at $`\text{Pr}_\text{m}1`$. We argue that, in the limit of $`\text{Re}\mathrm{}`$, a finite $`\text{Pr}_{\text{m,c}}`$ may exist. The second possibility is that $`\text{Pr}_{\text{m,c}}0`$ as $`\text{Re}\mathrm{}`$, while $`\text{Rm}_\text{c}`$ tends to a very large constant value inaccessible at current resolutions. If there is a finite $`\text{Pr}_{\text{m,c}}`$, the dynamo is sustainable only if magnetic fields can exist at scales smaller than the flow scale, i.e., it is always effectively a large-$`\text{Pr}_\text{m}`$ dynamo. If there is a finite $`\text{Rm}_\text{c}`$, our results provide a lower bound: $`\text{Rm}_\text{c}220`$ for $`\text{Pr}_\text{m}1/8`$. This is larger than Rm in many planets and in all liquid-metal experiments.

preprint: E-print astro-ph/0308336

The simplest description of a conducting fluid is in terms of equations of magnetohydrodynamics (MHD):

$`_t𝐮+𝐮𝐮`$ $`=`$ $`\nu \mathrm{\Delta }𝐮p+𝐁𝐁+𝐟,`$ (1)
$`_t𝐁+𝐮𝐁`$ $`=`$ $`𝐁𝐮+\eta \mathrm{\Delta }𝐁,`$ (2)
where $`𝐮`$ is velocity, $`𝐁`$ is magnetic field, $`𝐟`$ is the external-force density, $`\nu `$ is viscosity and $`\eta `$ is magnetic diffusivity. The pressure gradient $`p`$ is determined by the incompressibility condition $`𝐮=0`$. We have rescaled $`p`$ and $`𝐁`$ by $`\rho `$ and $`(4\pi \rho )^{1/2}`$, respectively ($`\rho `$ is density).

A fundamental property of Eqs. (12) is the ability of the velocity and magnetic fields to exchange energy. In three-dimensional turbulent flows (and in many other chaotic flows), this can take the form of net amplification of magnetic field with time, a process referred to as MHD dynamo. There are two kinds of dynamo. The first is the mean-field dynamo (growth of $`𝐁`$), which usually requires a net flow helicity Moffatt (1978). It is a large-scale effect that must be considered in conjunction with such system-specific properties as geometry, rotation, mean shear etc. The second kind is the small-scale dynamo: amplification of the magnetic energy $`B^2`$ due to random stretching of the field by the turbulent flow, requiring no net helicity Batchelor (1950); Zeldovich et al. (1984); Childress and Gilbert (1995). The stretching is opposed by the resistive diffusion, so the dynamo is only possible when the magnetic Reynolds number $`\text{Rm}=u^2^{1/2}\mathrm{}_0/\eta `$ exceeds a certain critical value ($`\mathrm{}_0`$ is the system scale).

In this Letter, we study the existence of small-scale dynamo in homogeneous incompressible turbulence with magnetic Prandtl number $`\text{Pr}_\text{m}=\nu /\eta <1`$ (i.e., $`\text{Re}=u^2^{1/2}\mathrm{}_0/\nu >\text{Rm}`$). This is an important issue because $`\text{Pr}_\text{m}`$ is small in stars ($`\text{Pr}_\text{m}10^2`$ at base of the Sun’s convection zone), planets ($`\text{Pr}_\text{m}10^5`$ Roberts and Glatzmaier (2000)), and in liquid-metal laboratory dynamos Gailitis et al. (2002); Forest et al. (2002); Bourgoin et al. (2002).

In three dimensions, most types of turbulence at scales much smaller than the system size are predominantly vortical and well described by Kolmogorov’s dimensional theory Frisch (1995). In this theory, the fastest field stretching is done by the small-scale velocities. The essential physics of small-scale dynamo should thus be contained within our homogeneous, isotropic, incompressible model.

The small-scale dynamo is most transparent in the limit of $`\text{Pr}_\text{m}1`$. Straightforward estimates show that, while velocity is dissipated at the viscous scale $`\mathrm{}_\nu \text{Re}^{3/4}\mathrm{}_0`$, magnetic field can occupy smaller scales down to the resistive $`\mathrm{}_\eta \text{Pr}_\text{m}^{1/2}\mathrm{}_\nu `$. The dynamo is driven by the fastest eddies: the viscous-scale ones, which are spatially smooth. The growing fields are organized in folds, with direction reversals at the resistive scale $`\mathrm{}_\eta `$ and field lines remaining approximately straight up to the flow scale $`\mathrm{}_\nu `$ Ott (1998); Chertkov et al. (1999); Schekochihin et al. (2002, 2003). Why this is a winning configuration is best seen on the example of a linear velocity field (locally uniform rate of strain) Zeldovich et al. (1984); Chertkov et al. (1999): the field aligns with the stretching direction of the flow but reverses along the “null” direction, so that compression cannot lead to resistive annihilation of antiparallel fields canceling the effect of stretching. For this mechanism to apply, it is essential that (i) the flow be spatially smooth, so fluid trajectories separate exponentially in time leading to exponential stretching, (ii) a scale separation between the field scale (reversals) and the flow be achievable. The large-$`\text{Pr}_\text{m}`$ turbulent dynamo satisfies both conditions, as do all deterministic chaotic dynamos Childress and Gilbert (1995). Thus, the small-scale dynamo, as it is usually understood, is essentially the large-$`\text{Pr}_\text{m}`$ dynamo. The often simulated case of $`\text{Pr}_\text{m}=1`$ belongs to the same class: the magnetic energy is amplified at scales somewhat smaller than the viscous scale and the field structure is similar to the case of $`\text{Pr}_\text{m}1`$ Schekochihin et al. (2003). When $`\text{Pr}_\text{m}<1`$, the field scale is resistively limited to be comparable to, or larger than, the viscous cutoff. The field interacts with inertial-range motions, which are spatially rough and cannot be thought of as having a locally uniform rate of strain. Is there still a small-scale dynamo?

In order to address this question, we have carried out a series of numerical simulations. Equations (12) were solved in a triply periodic box by the pseudospectral method. We used a random nonhelical forcing $`𝐟`$ applied at the box scale and white in time. The average injected power $`ϵ=𝐮𝐟`$ was kept fixed. The code units are based on box size 1 and $`ϵ=1`$. Defining $`\text{Rm}=u^2^{1/2}/\eta k_0`$ with $`k_0=2\pi `$ the box wave number, we have found that the dynamo existed for $`\text{Pr}_\text{m}1`$ provided $`\text{Rm}80`$ (cf. Fig. 1b). Once the existence of the dynamo at $`\text{Pr}_\text{m}=1`$ for a given $`\eta `$ was ascertained, $`\eta `$ was fixed. The viscosity $`\nu `$ was then decreased. The results are summarized in Fig. 1. In runs with $`\eta =2\times 10^3`$ (series A, $`\text{Rm}110`$), starting from weak field, the dynamo persisted at $`\text{Pr}_\text{m}=1/2`$, but was shut down at $`\text{Pr}_\text{m}=1/4`$. In runs with $`\eta =10^3`$ (series B, $`\text{Rm}=220\pm 10`$), there was dynamo at $`\text{Pr}_\text{m}=1/2`$ and $`1/4`$, but not at $`\text{Pr}_\text{m}=1/8`$. Figure 2 shows the time-averaged normalized magnetic-energy spectra for series A and B, as well as velocity spectra multiplied by $`k^2`$. The latter characterize the turbulent rate of strain and peak at the viscous scale. We see that as this scale drops below the resistive scale, the dynamo shuts down. Note that there is no initial-condition dependence: Runs A3 and B4, which decay starting from weak field, also decay if initialized in the saturated state of their $`\text{Pr}_\text{m}=1`$ counterparts (Fig. 1a).

Our main result is, thus, that $`\text{Pr}_{\text{m,c}}`$ exists even as Rm is kept approximately fixed at a value for which small-scale dynamo is possible at larger $`\text{Pr}_\text{m}`$ (Fig. 1b). In other words, the critical magnetic Reynolds number for growth $`\text{Rm}_\text{c}`$ increases with Re. Because of resolution constraints, we cannot afford a parameter scan to produce the dependence $`\text{Rm}_\text{c}(\text{Re})`$. What a priori statements about this dependence can be made on physical grounds?

Consider first the asymptotic case $`\text{Re}\text{Rm}1`$. The resistive scale then lies in the inertial range, $`\mathrm{}_0\mathrm{}_\eta \mathrm{}_\nu `$. As, in Kolmogorov turbulence, $`u_{\mathrm{}}/\mathrm{}\mathrm{}^{2/3}`$, most of the stretching is done by the eddies at the resistive scale $`\mathrm{}_\eta \text{Rm}^{3/4}\mathrm{}_0`$, where stretching is of the same order as diffusion Moffatt (1961); Vainshtein (1982). Since the inertial range is self-similar, the existence of the dynamo should not depend on the exact location of $`\mathrm{}_\eta `$, and it is the local (in $`k`$ space) properties of the turbulent velocity that determine its propensity to amplify magnetic energy. Therefore, if the dynamo fails, it does so at all $`\text{Pr}_\text{m}`$ below some critical value $`\text{Pr}_{\text{m,c}}`$ of order unity. The effective transition from the “large-$`\text{Pr}_\text{m}`$” to the “small-$`\text{Pr}_\text{m}`$” regime occurs at $`\text{Pr}_\text{m}=\text{Pr}_{\text{m,c}}`$. In this case, $`\text{Rm}_\text{c}/\text{Re}\text{Pr}_{\text{m,c}}=\text{const}<1`$ as $`\text{Re}\mathrm{}`$. Thus, if our results are asymptotic, then the turbulent small-scale dynamo is always, in essence, a large-$`\text{Pr}_\text{m}`$ one, and the folded direction-reversing fields are the only type of magnetic fluctuations that can be self-consistently generated and sustained by nonhelical turbulence. As the separation between parallel and transverse scales of the field diminishes at $`\text{Pr}_\text{m}<1`$ (Fig. 3), no steady fluctuation level can be maintained.

The second possibility is that $`\text{Rm}_\text{c}`$ asymptotes to some constant value for Re above those we are able to resolve: $`\text{Rm}_\text{c}\text{const}220`$ and $`\text{Pr}_{\text{m,c}}0`$ as $`\text{Re}\mathrm{}`$. Our results do not rule out this outcome, whereby asymptotically in Rm and Re, the dynamo persists at low $`\text{Pr}_\text{m}`$, but very large Rm is needed to achieve it in practice (numerically or experimentally). In stellar convective zones, Rm is, indeed, very large ($`10^6\mathrm{}10^9`$ for the Sun). On the other hand, in planets and in laboratory dynamos, $`\text{Rm}10^2`$, which is comparable to Rm in our simulations.

Note that the arguments above assume scale invariance of the inertial range, i.e., neglect the effects of intermittency. An intermittent velocity field will exhibit large coherent fluctuations of the rate of strain, which might be locally effective in stretching the magnetic field in a way similar to the large-$`\text{Pr}_\text{m}`$ dynamo 111This point was made to us by M. Vergassola.. Whether these fluctuations can provide enough stretching on the average to make a workable dynamo cannot be settled qualitatively. Note that an intermittent growth by rare strong bursts is evident in $`\text{Pr}_\text{m}<1`$ runs where the dynamo is suppressed but not shut down (most vividly in Run A2, see Fig. 1a).

No theory of the dynamo shutdown at low $`\text{Pr}_\text{m}`$ exists at present. Invoking turbulent diffusion (mixing) of magnetic fields by the subresistive-scale motions as the suppression mechanism makes heuristic sense, but does not provide an unambiguous verdict on the existence of the dynamo. Indeed, when $`\mathrm{}_0\mathrm{}_\eta \mathrm{}_\nu `$, the dominant contributions to both stretching and mixing are from the resistive scale $`\mathrm{}_\eta `$, and the outcome of their competition is impossible to predict on a heuristic basis. Many aspects of small-scale dynamo have received successful theoretical treatment in the framework of the Kazantsev model Kazantsev (1967), which assumes a Gaussian white-in-time velocity field, $`u^i(t,𝐱)u^j(t^{},𝐱^{})=\delta (tt^{})\kappa ^{ij}(𝐱𝐱^{})`$. The correlator can be expanded, $`\kappa ^{ij}(𝐲)\kappa ^{ij}(0)y^\alpha `$, when $`y\mathrm{}_B`$, the magnetic-field scale. When $`\text{Pr}_\text{m}1`$, $`\mathrm{}_B\mathrm{}_\nu `$, so $`\alpha =2`$, corresponding to the spatially smooth viscous-scale eddies. On the other hand, when $`\text{Pr}_\text{m}1`$, $`\mathrm{}_B\mathrm{}_\eta \mathrm{}_\nu `$ and magnetic field interacts with rough inertial-range velocities, which must be modeled by $`\alpha <2`$. The Kazantsev velocity is not a dynamo if it is too spatially rough, viz., when $`\alpha <1`$ Kazantsev (1967); Vergassola (1996); Rogachevskii and Kleeorin (1997); Vincenzi (2002). If we set aside fundamental objections to the $`\delta `$-correlated model and try to compare it to real turbulence, we still face the difficulty of interpreting the $`\delta `$ function. If we write equal-time velocity correlators by replacing the $`\delta `$ function by inverse correlation time $`1/\tau _\text{c}`$, then the relation between $`\alpha `$ and the spectral exponent of the turbulence depends on how we choose $`\tau _\text{c}`$. The usual choice for Kolmogorov turbulence is $`\tau _\text{c}y^{2/3}`$, the eddy-turnover time. Then $`\alpha =4/3>1`$ and there is dynamo Vainshtein (1982)222Note that the EDQNM closure, which effectively makes the same assumption about the correlation time, gives dynamo at low $`\text{Pr}_\text{m}`$, with $`\text{Rm}_\text{c}`$ independent of Re \[J. Léorat et al., J. Fluid Mech. 104, 419 (1981)\].. Note that $`\text{Rm}_\text{c}`$ in this case is typically much larger than for $`\alpha =2`$ Rogachevskii and Kleeorin (1997); Vincenzi (2002); Boldyrev and Cattaneo (2003). Although specific values of $`\text{Rm}_\text{c}`$ calculated from the Kazantsev model cannot be considered as quantitative predictions for real turbulence, they appear to point to the second possibility mentioned above (finite but unresolvably large $`\text{Rm}_\text{c}`$). We emphasize that all these results depend on the heuristic choice of $`\tau _\text{c}`$ (e.g., if $`\tau _\text{c}\text{const}`$, $`\alpha =2/3`$ and there is no dynamo) and on the universality of the physically nonobvious condition $`\alpha >1`$. It is fair to observe that our simulations are still too viscous to have a well-developed Kolmogorov scaling (Fig. 2). Thus, if the existence of the dynamo depends on the exact inertial-range scaling of the velocity field and/or only manifests itself at very large Rm, neither the Kazantsev theory nor simulations at current resolutions can lay claim to a definitive answer. Obviously, the Kazantsev theory also cannot capture any role the intermittency of the velocity field might play and, more generally speaking, it is doubtful that a $`\delta `$-correlated Gaussian flow is a suitable model of the inertial-range turbulence.

While, as far as we know, ours is the first systematic study of the small-scale dynamo suppression in homogeneous isotropic MHD turbulence with low $`\text{Pr}_\text{m}`$, indications of this effect have been reported in the literature in two previous instances. Dynamo suppression at low $`\text{Pr}_\text{m}`$ was seen by Christensen et al. Christensen et al. (1999) in their simulations of convection in a rotating spherical shell and by Cattaneo Cattaneo (2002) in simulations of Boussinesq convection. These cases of failed low-$`\text{Pr}_\text{m}`$ dynamos in inhomogeneous convection-driven turbulence are likely to be related to the same universal mechanism that made the dynamo inefficient in our simulations. Our key conclusion is that the dynamo suppression is a generic effect unrelated to the particular type of driving or other large-scale features of the system.

The mean-field dynano, which, in contrast, does depend on large-scale features such as helicity and rotation Moffatt (1978), may then be the only type of self-sustained field generation for low-$`\text{Pr}_\text{m}`$ systems. If a mean field is present, it gives rise to a source term $`𝐁𝐮`$ in the induction equation (2) and thus induces small-scale magnetic fluctuations. They have a $`k^{11/3}`$ spectrum at $`kk_\eta `$ Golitsyn (1960); Moffatt (1961); Vainshtein (1982), which has been seen in the laboratory Bourgoin et al. (2002) and in large-eddy simulations Ponty et al. (2003).

###### Acknowledgements.

We thank S. Boldyrev and F. Cattaneo for discussions of our results and of the Kazantsev model, and M. Vergassola for very useful comments. Our work was supported by PPARC (PPA/G/S/2002/00075), UKAEA (QS06992), NSF (AST 00-98670). Simulations were run at UKAFF (Leicester) and NCSA (Illinois).


text:
# Millimeter VLBI and Variability in AGN Jets

## 1. Present Status of 3mm-VLBI

The main advantage of VLBI observations at short millimeter wavelengths (mm-VLBI) is the very high angular resolution and the ability to study emission regions, which are self-absorbed at longer wavelengths. The small observing beam (up to $`50`$$`\mu `$as at 86 GHz) not only provides a spatial resolution of a few hundred to thousand Schwarzschild radii for a typical $`z>0.1`$ quasar, but also facilitates to locate and trace the motion of moving jet components with unprecedented accuracy and at early times after initial flux density outbursts.

To date, global VLBI observations at 86 GHz are performed regularly, however, with regard to the relatively fast motion seen in many quasars ($`0.11`$ mas/yr), not yet with dense enough time sampling (typically 2 observations per year). The worldwide combination of stations in the US – like the VLBA (presently 7 antennas equipped with receivers) and Haystack (37 m) – and sensitive antennas in Europe (Effelsberg 100 m, Pico Veleta 30 m Onsala 20 m and Metsahovi 14 m) forms the most powerful global mm-VLBI array. With a typical single baseline detection sensitivity ($`7\sigma `$) of typically $`0.2`$ Jy between Effelsberg and Pico Veleta and $`0.5`$ Jy on VLBA baselines, sources with $`\mathrm{S}_{86\mathrm{GHz}}12`$ Jy can be reliably imaged with a dynamic range of up to a few hundred. In Figure 1 we show examples of 3 mm-VLBI images for the quasar 3C 345 and for BL Lac.

## 2. VLBI at Higher Frequencies

Although sensitivity limitations and remaining calibration uncertainties still restrict the imaging capabilities of VLBI in the 3 mm band a bit, such VLBI observations nowadays can be regarded as more or less standard.

Efforts to push VLBI to even shorter wavelengths started in the late 1980’s and led to various successful detections of bright AGN at 1.3 mm, however, only on relatively short (continental) VLBI baselines (e.g. Greve et al. 1995, Krichbaum et al. 1997). Between 1995 – 2000 several attempts with various telescopes were made to achieve fringe detections also on the longer transatlantic baselines. These experiments were performed in the 2 mm and 1.3 mm bands, but failed due to technical difficulties. The recent promising detection of 3C 273 and 3C 279 at 147 GHz on the 3100 km (1.5G$`\lambda `$) baseline between Pico Veleta (Spain) and Metsähovi (Finland) in spring 2001 (Greve et al. 2002), and the availability of VLBI equipment and a new 2 mm receiver at the Heinrich Hertz telescope (HHT) on Mt. Graham (Arizona), stimulated a new transatlantic VLBI experiment at 147 GHz, which was performed in April 2002. We now are able to report fringe detections also on the very long inter-continental VLBI baselines (Krichbaum et al. 2002). The detection of the 3 sources NRAO 150, 3C 279, and 1633+382 on the VLBI baselines between Arizona and Spain at a fringe spacing of 4.2 G$`\lambda `$ (corresponding to $`49`$$`\mu `$as resolution) with signal-to-noise ratios of up to 75 (for 3C 279) demonstrates for the first time that global VLBI at 2 mm-wavelength is possible and that the previous technical and sensitivity limitations can be overcome.

The present 2 mm and 1 mm VLBI observations were performed in snapshot mode (only few VLBI scans available per source) and it is therefore not yet possible to make maps of the detected sources. However, some estimates for the source sizes can be made. In Figure 2 we show the normalized visibility amplitude (=correlated flux / total flux) of 3C 279 plotted versus the projected baseline length (data for April 2002: filled circles). The comparison of the flux density between long and short VLBI baselines (assuming a simple Gaussian point source structure) gives a size of the VLBI core of 3C 279 of $`35\pm 5`$$`\mu `$as at 147 GHz. For comparison, we also show visibility points from the previous 2 mm VLBI observation of 2001 (filled square) and the 1.3 mm VLBI observation of 1995 (open diamond). These data points indicate that the source structure is most likely variable in time and probably not pointlike, both of which is not unexpected knowing the source properties at longer wavelengths.

It is possible that the factor of two difference between the visibility amplitudes near and below 1 G$`\lambda `$ reflects a temporal variation of the core size between $`2001`$ and 2002. It is also possible that the lower amplitude in this uv-range, is due to structural beating introduced by the known sub-mas jet (1 G$`\lambda `$ corresponds to 0.2 mas). A second estimate of the source size therefore can be made, if we combine the visibility points from all 3 experiments (dotted line). This yields a more compact source size of $`25`$$`\mu `$as.

For an incoherently radiating compact synchrotron source with a brightness temperature close to the inverse Compton limit of $`10^{12}`$ K, the theoretically expected size is 10–20 $`\mu `$as, which is about a factor of 1.5–2 smaller than the two sizes derived above. It is therefore possible that one sees at 2 mm wavelength the VLBI core component of 3C 279 (the jet base) already spatially resolved. In this case, its intrinsic size would be of order of 0.2–0.3 pc. Of course, this estimate depends critically on the accuracy of the calibration of the VLBI antennas. Future experiments with an improved gain calibration will be necessary to check, if the rather low visibility amplitudes at G$`\lambda `$ uv-spacings (only $`1030`$ % of the total flux density is seen) are real.

## 3. Variability and new jet components in 3C 273

For many AGN, a correlation between flux density variability and ejection of VLBI components is suggested. From our multi-frequency VLBI monitoring of 3C 273 (1987 – 1997, Krichbaum et al. 2002), we identified 13 jet components (C6 – C18). These components move with apparent speeds of 4–8 c and seem to accelerate as they separate from the VLBI core. Our high frequency VLBI data (15-100 GHz) allow to trace the component motion back to the ejection from the VLBI core and by this to determine the times of zero separation from the core. The typical measurement uncertainty for the ejection times $`t_0`$ ranges between 0.2 – 0.5 yr. In Figure 3 (left), we plot $`t_0`$ and the millimeter-variability (22 – 230 GHz, data: H. Ungerechts, H. Teräsranta). We also add the Gamma-ray detections of 3C 273 from EGRET. In Figure 3 (right) we plot the onset times of the mm-flares derived from these light-curves (Türler et al. 1999) together with the VLBI ejection time ($`t_0`$) and the Gamma-ray fluxes. Except for component C15, which may be not related to a flux density outburst, a clear correlation between component ejection and onset times of mm-flares is seen.

For the Gamma-rays, the search for a correlation with mm-flares and VLBI component ejections is complicated by the coarse time sampling of the Gamma-ray light-curve, which does not exclude presence of more rapid variations. However, since at least 3 of the Gamma-ray detections shown in Figure 3 (left) are followed by lower Gamma-ray fluxes (with maxima near 1991.5, 1993.9, 1997.0), it is tempting to identify these with local maxima in a heavily undersampled light-curve. With this assumption, it is possible to assign the ejection time of component C12 to the Gamma-flare of 1991.5, C13 to a flare peaking somewhere between 1993.0 and 1993.9 and C18 to the flare peaking near 1997.0. For the Gamma-flare of 1997.0 it is very likely that the true maximum of the light-curve has been detected, since this measurement is surrounded by two measurements of lower intensity. For the components C14 and C15–C17 the situation is less clear. We note, however, that the components C15–C16–C17 were ejected nearly at the same time and are probably physically related to each other (shock-post-shock structure). Therefore, one of the 3 components is likely to be related with the Gamma-flare peaking around 1996.0. For C14, the ejection time correlates well with a mm-flare (peak: 1994.9) and a nearly contemporaneous Gamma-ray detection of moderate strength.

From a more detailed analysis (Krichbaum et al. 2002) we obtain for the time lag between component ejection and onset of a mm-flare: $`t_0t_0^{\mathrm{mm}}=0.1\pm 0.2`$ yr. If we assume that the observed peaks in the Gamma-ray light-curve are located near the times $`t_0^\gamma `$ of flux density maxima, we obtain $`t_0^\gamma t_0^{\mathrm{mm}}=0.3\pm 0.3`$ yr. Although the Gamma-ray variability may be faster, this result is fully consistent with the more general finding of enhanced Gamma-ray fluxes mainly during the rising phase of millimeter flares. Thus, we suggest the following sequence of events: $`t_0^{\mathrm{mm}}t_0t_0^\gamma `$ – the onset of a millimeter flare is followed by the ejection of a new VLBI component and, either simultaneously or slightly time-delayed, an increase of the Gamma-ray flux. If we focus only on those VLBI components, which were ejected close to the main maxima of the Gamma-ray light-curve in Figure 3, we obtain time lags of $`t_0^\gamma t_0`$ of $`0.5`$ yr for C12, $`0.9`$ yr for C13, $`0.2`$ yr for C16 and $`0.1`$ yr for C18. In all cases the Gamma-rays seem to peak a little later than the time of component ejection. With $`\beta _{app}4`$ near the core, the Gamma-rays would then escape at a radius $`r_\gamma 0.1`$ mas. This corresponds to $`r_\gamma 2000`$ Schwarzschild radii (for a $`10^9`$ M$`_{\mathrm{}}`$ black hole) or $`610^{17}`$ cm, consistent with theoretical expectations, in which Gamma-rays escape the horizon of photon-photon pair production at separations of a few hundred to a few thousand Schwarzschild radii.

## 4. Summary and Future Outlook

In summary, global mm-VLBI observations of Blazars (QSO’s, BL Lac’s) generally show a very compact core radiating near the inverse Compton limit, but also a considerable ($`3050`$ %) amount of flux and jet-like sub-structure on the sub-mas- to mas scale. VLBI observations at 2 & 1 mm indicate that the brightness temperatures of the VLBI cores of mm-bright blazars is not significantly different from those, observed at longer wavelengths. The VLBI-jets often are helically bent and in many cases the jet curvature increases towards the VLBI-core. There is evidence that the helicity of the jets and the component motion along spatially curved paths is related to internal rotation of the jet and precession at the jet base. In 3C 273, 3C 345, BL Lac, 0716+714 and many other sources, new superluminally moving jet components appear weeks to months after the $`\mathrm{𝑂𝑁𝑆𝐸𝑇}`$ of flux density outbursts, which propagate through the whole electromagnetic spectrum. In 3C 273 Gamma-ray flares and the ejection of relativistic jet components are closely related.

The relatively fast motion and the complexity of the source structures requires for the future more sensitive mm-VLBI antennas, which observe with denser (weekly to monthly) time sampling. Present day mm-VLBI suffers from the lack of short uv-spacings. As a consequence, it is presently not possible to reliably image all of the source flux visible at mm-wavelengths. The relatively high surface brightness of the sources on the short (100-1000 km long) baselines facilitates easy VLBI detection and therefore gives room also for less sensitive and smaller antennas to play a significant role in future high resolution mm-VLBI imaging.

The addition of more collecting area through sensitive antennas specially designed for the millimeter and sub-millimeter bands is most important for future mm-VLBI, particularly at the shorter wavelengths (2 & 1 mm). This includes existing antennas, which are not yet participating in mm-VLBI (eg. JCMT, SMA, NRO), but also new antennas like e.g. the 50 m LMT in Mexico and the ALMA prototype antenna APEX. Owing to their outstanding sensitivity, phased interferometers will play a particular important role. For the near future (2002/3), the participation of the IRAM interferometer at Plateau de Bure (France) is planned (first at 3 & 1 mm, later also at 2 mm). This will increase the present sensitivity by a factor of $`23`$ to the $`0.1`$ Jy level. In the more distant future CARMA - the merger of the BIMA and OVRO interferometers - and ALMA should lower the detection threshold to 1–10 mJy and by this largely enhance the observational possibilities. In parallel, higher data recording rates and observing bandwidths (GBit/s recording using the MKV system) and correction of atmospheric phase degradation (water vapor radiometry, phase referencing) will help to further improve the sensitivity.

All this will facilitate the imaging of Quasars, nearby Radio-Galaxies and of the Galactic Center source with the fascinating angular resolution of only $`1020`$$`\mu `$arcseconds !

## References

Greve, A., et al., 1995, A&A, 299, L33.

Greve, A., et al., 2002, A&A, 390, 19.

Krichbaum, T.P., et al. 1997, A&A, 323, L17.

Krichbaum, T.P., et al. 2002, in: Proc. of the 6th European VLBI Network Symposium, eds. E. Ros, R.W. Porcas, A.P. Lobanov, & J.A. Zensus (Bonn: MPIfR), 125.

Krichbaum, T.P., et al. 2002, A&A, submitted.

Türler, M., et al. 1999, A&A, 349, 45.


text:
# Biopolymer Elasticity

## Abstract

In recent years molecular elasticity has emerged as an active area of research: there are experiments that probe mechanical properties of single biomolecules such as DNA and Actin, with a view to understanding the role of elasticity of these polymers in biological processes such as transcription and protein-induced DNA bending. Single molecule elasticity has thus emerged as an area where there is a rich cross-fertilization of ideas between biologists, chemists and theoretical physicists. In this article we present a perspective on this field of research.

Polymers of different flexibilities are vital to life. Some polymers are as stiff as needles, others are as flexible as thread and still others are like twine, neither too stiff nor too floppy. The cytoskeletal structure of the cell, which gives rigidity to the cell, consists of semiflexible polymers - like Actin, Microtubules and intermediate filaments. The polymer which carries the genetic code, the DNA, is a semiflexible polymer (See Box $`1`$). A typical human DNA in a cell is about a meter long and it is rather amazing that such a long molecule is packed into a cell nucleus of the size of a few microns. This efficient packaging is a mystery and raises questions about the elasticity of DNA. DNA also assumes various different forms depending on the function it needs to perform \- for instance, it is unwound during replication so that the genetic code can be accessed and copied. At other times the DNA assumes a supercoiled configuration which protects and preserves the genetic code intact till it is needed. Such conformationally distinct configurations are also dictated by its elasticity.

Single Molecule Experiments:

Experimental studies of biopolymer molecules such as DNA have traditionally been limited to samples containing large numbers of molecules. This made it hard to probe the elastic properties of individual biopolymers which are of vital importance to biological processes such as protein-induced DNA bending. It is only quite recently, due to advances in technology that single molecule studies became feasible. It is now possible to design experiments in which single molecules are pulled, stretched and twisted to measure elastic properties. In a typical single molecule experiment, one considers a polymer molecule suspended between a fixed surface and a force sensor of some kind (See Fig. $`1`$). The force sensor is a bead in a laser trap or a flexible cantilever of an atomic force microscope and the molecule under consideration can range from a DNA in a simple random coil configuration to a globular protein in a unique three-dimensionally folded form. For instance, one can study the “equation of state” of a semiflexible polymer by measuring its extension as a function of applied force. One can also tag the ends with fluorescent dye and determine the distribution of end-to-end distances. We learn a lot from such experimental studies about biologically relevant mechanical properties of these polymers.

Such single molecule micromanipulation studies have opened up a vast area of cross-fertilization of ideas between biologists, chemists and theoretical physicists. Results of single molecule experiments have posed challenges to theoretical physicists who have combined ideas from statistical mechanics, quantum mechanics, differential geometry and topology to get a proper understanding of the field (See Box $`2`$).

Force-Extension Curves: Experiment And Theory Carlos Bustamante and his colleagues directly measured the elasticity of a single DNA molecule by attaching one end of the molecule to a glass slide and the other end to a micron size magnetic bead. A stretching force was applied to the DNA molecule by applying a magnetic field gradient and the corresponding extension was measured by observing the position of the bead under an optical microscope. The mechanical response of such a system was thus probed in the form of a force-extension curve.

The experimental curves were compared first against the simplest theory available for a polymer- the freely jointed chain (FJC) in which the polymer is viewed as a chain of rigid rods connected by revolving pivots. While the experimental curves agreed well with the FJC model at low forces, the comparison showed considerable discrepancies between theory and experimental curves in the intermediate and large force regime. This pointed towards a deficiency in the theory which was subsequently remedied by Marko and Siggia who used the Worm Like Chain (WLC) model. In the WLC model one takes into consideration the bending energy ($`=\frac{1}{2}A_0^L𝑑s\kappa ^2`$ with $`A=L_pk_BT`$ the bending modulus) proportional to the square of the curvature $`\kappa `$ of the space curve representing the polymer of contour length $`L`$ kept at a temperature $`T`$. The WLC model provided excellent quantitative agreement with the experimental force-extension curves (See Fig. $`2`$). This pointed to the fact that a DNA molecule needs to be viewed as a semiflexible (partly flexible and partly stiff) polymer rather than a completely flexible one.

Twisting Single Molecules:

In the experiments of Strick and collaborators the ends of a single molecule of double stranded DNA are attached to a glass plate and a magnetic bead. Magnetic fields are used to control the orientation of the bead and magnetic field gradients to apply forces on the bead. By such techniques the molecule is stretched and twisted and the extension of the molecule is monitored by the location of the bead. A typical experimental run measures the extension of the molecule as a result of the applied twist and force. These experiments are important to understanding the role of twist elasticity in biological processes like DNA replication.

In laboratory and biological conditions the thermal fluctuations of the DNA molecule are important and it is necessary to take these into account to understand the experiments of Strick et al. To determine the response of a DNA molecule to twist and stretch one has to calculate its partition function: sum over all possible configurations of the molecule with Boltzmann weight $`\mathrm{exp}[E/kT]`$. As Strick et al remark that is a tall order, not likely to be filled anytime soon. Computer simulations of the system are possible and have been done. These simulations, which incorporate realistic features like self-avoidance, agree well with laboratory experiments and can be viewed in two ways. From the point of view of a laboratory experimenter they are like theoretical models. From a theorist’s standpoint simulations can be viewed as controlled experiments. The main theoretical difficulty is the incorporation of self-avoidance, which is present in a real DNA molecule: the conformation of the molecule must never intersect itself. Such a constraint is virtually impossible to handle theoretically, and one is forced to resort to simpler models. Nevertheless, researchers have attempted to model such a system by incorporating a twist term in the energy functional and trying to explain the results of the experiments . Although the agreement with the experimental curves is good there are some issues related to self-avoidance of the polymer which are still poorly understood and one needs to introduce an adjustable cutoff parameter to fit the experimental curves. The significance of this cutoff parameter is not entirely clear and there are open issues yet to be understood.

Finite Size Effects in Single Molecule Experiments

In the past, experiments on biopolymers were confined to studying their bulk properties, which involved probing large numbers of molecules. The results of these experiments could be analyzed by using the traditional tools of thermodynamics. However, when one micromanipulates single polymer molecules the usual rules of thermodynamics applicable to bulk systems no longer hold. Single molecule experiments thus provide physicists with a concrete testing ground for understanding some of the fundamental ideas of statistical mechanics. In particular, fluctuations about the mean value of a variable play an important role due to the finite extent of the molecule. For instance, it turns out, that an experiment in which the distance between the ends of a polymer molecule is fixed (an isometric setup) and the tension fluctuates yields a different “equation of state” from one in which the tension between the ends is held fixed (an isotensional setup) and the end-to-end distance fluctuates. This asymmetry can be traced to large fluctuations about the mean value of the force or the extension, depending on the experimental setup. These fluctuations vanish only in the thermodynamic limit of very long polymers. A similar effect is seen in experiments involving twisting DNA molecules. If the DNA is short enough, experiments with a fixed torque and a fluctuating twist give a different Torque-Twist Relation from those with a fixed twist and a fluctuating torque.

Equilibrium and Nonequilibrium Statistical Mechanics: A Connection

We will now come to an interesting connection between equilibrium and nonequilibrium statistical mechanics which has emerged in the past few years and its relevance to single molecule experiments.

In recent years Jarzinsky has proved a remarkable equality relating work done in a series of finite time, nonequilibrium measurements and the equilibrium free energy difference between two given configurations of a system. What is most remarkable is, unlike linear response and the fluctuation-dissipation theorem (See Box $`3`$) which relate equilibrium statistical mechanics to nonequilibrium processes close to equilibrium, this relation is valid arbitrarily far from equilibrium.

Let us consider a finite classical system that depends on an external parameter $`\lambda `$. Let the system come to equilibrium with a reservoir at a temperature $`T`$. If we now switch the external parameter from an initial value $`\lambda =0`$ to a final value $`\lambda =1`$ infinitely slowly then the work $`W_{\mathrm{}}`$ done on the system is given by:

$$W_{\mathrm{}}=\mathrm{\Delta }F=F_1F_0$$

where $`F_\lambda `$ is the Free energy of the system at a temperature $`T`$ at a fixed value $`\lambda .`$ Let us now consider switching the parameter $`\lambda `$ at a finite rate. In such a situation the system is no longer in quasistatic equilibrium with the reservoir at every stage. The work done in such a case depends on the microscopic initial conditions of the system and the reservoir. One now needs to consider an ensemble of such switching measurements and the total work done is no longer equal to the free energy difference but always exceeds it.

It turns out, however, that a remarkable equality still holds connecting the work done in such a process and the difference between the initial and final values of the free energy. The equality, known as Jarzynski’s equality is given by: $`\overline{exp(\beta W)}=exp(\beta \mathrm{\Delta }F)`$. The overbar on the left pertains to an ensemble average over all possible finite time measurement realizations.

Recently Liphardt and collaborators have tested the validity of this remarkable equality via single molecule experimental measurements. They tested the equality by comparing work done in reversible and irreversible unfolding of a single RNA molecule.

Small beads were attached at the ends of a single RNA molecule. The bead at one end was held in an optical trap. By measuring the displacement $`\delta `$ of the bead in the trap and the stiffness constant $`k`$ of the trap, the force $`F=k\delta `$ on the molecule was measured. The bead at the other end was held by a micropipette and attached to a piezoelectric device which enabled the position of the bead to be controlled by varying the voltage. The position $`z`$ of the second bead relative to a given reference point was used to measure the extent of unfolding of the molecule. The work $`W`$ done was measured by integrating the product of the force and the displacement and plotted as a function of the extension $`z`$.

The experiment was repeated at varying rates of pull. There were essentially two classes of measurements: quasistatic, near equilibrium measurements where the work done was reversible and given by the free energy difference between the two end values of the control parameter. The other class consisted of measurements carried out at faster rates of pull (truly nonequilibrium situation) where the molecule does not have time to equilibriate during pulling. The experimenters took an average over all possible realizations of the nonequilibrium measurements and evaluated the quantity $`\overline{exp(\beta W)}`$ and compared it against its equilibrium counterpart $`exp(\beta \mathrm{\Delta }F)`$ and found good agreement within experimental error corroborating Jarzynski’s equality.

Concluding Remarks:

In this article we have tried to give a perspective on the growing and exciting area of biopolymer elasticity where there has been a merging of ideas between biologists, chemists and theoretical physicists. On one hand, elasticity experiments on biopolymers triggered by a need to understand its role in complex biological processes such as transcription and gene regulation have inspired physicists to come up with models that capture the essence of the problem and at the same time give quantitative explanations for force-extension and twist-extension measurements. On the other hand, predictions made by theoretical physicists have suggested new experiments. Computer simulations have played a two fold role: they have provided researchers with results of “controlled experiments” which are ideal for theoretical model building. At the same time, the results of the simulations when tested against real experiments enable experimenters to identify the key factors which control the experimental results. Apart from biological relevance, biopolymer elasticity experiments have provided a testing ground for theoretical physicists for deeper issues related to connection between equilibrium and nonequilibrium statistical mechanics.

BOX I: Mean Square End-To-End Separation Of A Randomly Lying Cord

Typical chromosomal DNA molecules are long (contour length $`L50mm`$) and when viewed under a microscope look like random coils and the statistics of their ensemble of equilibrium configurations is like that of a random walk. The mean square distance between the ends of a molecule grows in proportion to its contour length $`L`$:

$$<R^2>=2L_pL.$$

The constant $`L_p`$ is the persistence length which is the length over which the DNA molecule can be considered to be approximately straight.

In general, for a semiflexible polymer, the mean square separation between its ends predicted by the most popular theoretical model of biopolymer elasticity, the Worm Like Chain model is:

$$<R^2>=2L_p^2[L/L_p+e^{L/L_p}1]$$

which in the limit of $`L/L_p>>1`$ as for a DNA, goes as

$$<R^2>=2L_pL.$$

In the other limit, that is, for a stiff biopolymer like an Actin filament the mean square separation between ends grows ballistically with the contour length

$$<R^2>L^2$$

.

There is an experiment that one can do using a fairly low-tech setup to check the statistics of end-to-end separation of polymers. This ‘scaled up’ experiment tells us that the statistics of separation between the ends of a cord is the same as one gets for a semiflexible polymer. The original experiment due to Lemons and Lipscombe (See Suggested Reading List) consisted of taking some short ($`0.1m`$) and some long ($`0.8m`$) segments of a bungee cord of length ($`=0.8m`$) and persistence length ($`0.25m`$). To mimic the effect of randomization, the cord was dropped onto the floor from a height of $`2m`$. Positions were marked off on the cord corresponding to different segment lengths $`L`$. After dropping the cord, the distance $`R`$ of separation between the $`L=0`$ end and one of these marks was measured and for each such segment length the procedure was repeated about $`40`$ times. $`<R^2>`$ was then plotted against $`L`$ and the results agreed well with the theoretical expression

$$<R^2>=2L_p^2[L/L_p+e^{L/L_p}1].$$

The persistence length or the rigidity parameter $`L_p`$ was determined from the fit of the data against the theoretically predicted form.

BOX 2: Link=Twist+Writhe

DNA can form a closed loop. Since it is a double stranded helix, the linking number $`Lk`$ or the number of times one strand wraps around the other is fixed. In 1971 Fuller noticed an interesting connection: the linking number $`Lk`$ of a closed ribbon can be decomposed into the “Writhe” $`Wr`$ of its backbone plus a locally determined “Twist”$`Tw`$ :

$$Lk=Tw+Wr.$$

Link $`Lk`$ counts the number of times that one strand winds around the other. It is an integer and clearly a topological quantity and does not change unless you cut and paste the strands. Writhe $`Wr`$ is a nonlocal quantity which measures the wandering of the tangent vector to the polymer backbone and it depends only on the central curve or the backbone of the ribbon representing the polymer. Twist $`Tw`$ is a local quantity that measures the integrated angular velocity of a strand about the backbone of the polymer. Both $`Tw`$ and $`Wr`$ can change so as to keep $`Lk`$ unchanged. Fuller noticed this relation in the context of a question related to supercoiled double stranded DNA rings raised by J. Vinograd, a molecular biologist. This is an example where molecular biology has stimulated a mathematical line of enquiry.

BOX 3: Fluctuation-Dissipation Theorem

When you turn on a radio you can sometimes hear a noise which can be traced to the irregular motion of electrons. Nyquist was the first to recognize the connection between such thermal noise and the impedance of a resistor across which an irregular voltage difference is induced due to thermal motion of the electrons. One sees a similar connection in the context of Brownian motion ( See S. Ramaswamy’s article mentioned in the Suggested Reading List). The generalized theorem relating the fluctuations of a physical variable of a system in equilibrium to dissipative processes in the system subjected to an external force driving it slightly away from equilibrium, is called, the Fluctuation-Dissipation Theorem. The Fluctuation-Dissipation Theorem works in the linear response regime in which the response of a system to an external stimulus is directly proportional to the strength of the stimulus.

Suggested Reading:

1. F. B. Fuller, Biophysical Journal 68, 815 (1971).
2. M. Doi and S. F. Edwards, The Theory of Polymer Dynamics, (Clarendon, Oxford, 1992).
3. J.Marko and E.D.Siggia, Macromolecules 28, 8759(1995).
4. T. R. Strick et al, Science 271, 1835 (1996).
5. C. Bouchiat et al, Biophysical Journal 76, 409 (1999).
6. C. Bustamante et al, Current Opinion in Structural Biology 10, 279 (2000).
7. S. Ramaswamy, Resonance 5, 16 (2000).
8. D. S. Lemons and T. C. Lipscombe, Am. J. Phys 70, 570 (2002).
9. J. Liphardt et al, Science 296, 1832 (2002).
10. D. Keller, D. Swigon and C. Bustamante, Biophysical Journal 84, 733 (2003).