text:
# The Riemann tensor for nonholonomic manifolds

## Introduction

Let $`W_k^n`$ be the space of germs of irreducible $`k`$-dimensional distributions at the origin of $`^n`$. The orbits of $`\mathrm{Diff}^n`$-action on $`W_k^n`$ are usually of infinite codimension. As shown in \[VG1\], there are 3 types of exceptions: (1) $`W_k^n`$ is a point, (2) codimension 1 distributions determined as zero set of $`\alpha =0`$ for the forms $`\alpha `$ of maximal class (recall that the class of the form is the invariant equal to the least number of variables actually involved)

$$\alpha =\{\begin{array}{cc}dt+\underset{1ir}{}(p_idq_iq_idp_i)\hfill & \text{for }n=2r+1\hfill \\ dt+\underset{1ir1}{}(p_idq_iq_idp_i)\hfill & \text{for }n=2r\hfill \end{array}$$

called (usually, for $`n`$ odd) a contact structure (see \[Z\]) and (3) the Engel distribution locally given by the system of Pfaff equations $`\alpha _1=0`$, $`\alpha _2=0`$, where

$`(E)`$
$$\alpha _1=dx_4x_3dx_1,\alpha _2=dx_3x_2dx_1.$$

Any 4-dimensional manifold with the Engel distribution (E) will be called an Engel manifold.

Both contact and Engel manifolds are examples of nonholonomic manifolds (i.e., manifolds with nonintegrable distributions). Recall that the term nonholonomic was coined by H. Hertz \[H\] in his attempts to geometrically describe motions in such a way as to eliminate the notion of “force”. For a historical review see \[VG2\]; see also a very interesting paper by Vershik \[V\] with first rigorous mathematical formulations of nonholonomic geometry and indications to various possible applications in various, sometimes unexpected, areas (like optimal control or macro-economics). In \[V\] Vershik summarizes about 100 years of studies of nonholonomic geometry (Hertz, Carathéodory, Vrănceanu, Wagner, Schouten, Faddeev, Griffiths, Godbillon, to name a few; to his list I’d add that MathSciNet returns about 700 entries for nonholonomic, and for synonyms: 70 more for sub-Riemannian,1700 for Finsler, several dozens for “cat’s problem” and autoparallel; physics have their own simialr lists with thousands of entries because there seems to be more, actually, nonholonomic dynamical systems than holonomic ones; finally supergravity is also a nonholonomic structure, albeit on supermanifolds — ca 1800 more entries). In \[V\] Vershik also declares that “it is probably impossible to construct an analog of the Riemann tensor for the general nonholonomic manifold” though in some particular cases of small dimension they had been constructed (Martinet) (and, one can add now, in papers on supergravity).

Here I introduce the invariants that characterize “nonflatness” of nonholonomic manifolds, i.e., positively solve this old problem. This definition appeared around 1989 but publication of this paper was delayed and the definition first (as far as I know) appeared in contetext of supergravity in \[LP\]. For application to SUGRA (supergravity equations) and other instances where the N-curvature tensors (the left hand side of SUGRA) are calculated see \[GL1\], \[LPS\], \[LSV\].

The general definition is illustrated with contact and Engel manifolds. As expected, for the contact structure these invariants vanish, thus giving one more proof of the Darboux theorem on the canonical form of the contact form; one can show that they vanish for $`n`$ even as well. The Engel manifolds are more rigid: there could be nonflat ones among them.

Observe that our invariants, though natural analogs of the curvature and torsion tensors, do not coincide on nonholonomic manifolds with the classical ones and bearing the same name. Indeed, by definition on any nonholonomic manifold, there is a nonzero classical torsion while every contact manifold is flat in our sense. To avoid confusion, we should always add adjective “nonholonomic” for the invariants introduced below. Since this is too long, we will briefly say N-torsion, N-curvature, etc.

For manifolds with a $`G`$-structure, the prototypes of the invariants we consider are well-known: these are structure functions of the $`G`$-structure \[St\]; an important particular case is $`G=O(n)`$ when the structure function is called the Riemann tensor.

A theorem of Lie states that there are two types of differential equations, considered as submanifolds in the jet bundles: (1) those whose symmetries are induced by point transformations and (2) those whose symmetries are induced by contact transformations. For equations of the first type Goldschmidt gave criteria for formal integrability \[BCG\]. The notions introduced below — generalizations of Spencer cohomology for nonholonomic case — provide with a straightforward generalization of Goldschmidt’s criteria for formal integrability to the second half of differential equations.

Finally, observe that “non-flatness” of distributions with non-linear constraints (such as fields of cones \[A\], foreseen by Vershik in \[V\], or fields of spheres, like the ones constructed by any driver who switched the cruise control ON, can also be tested by invariants I introduce; in the case of the fields of solids, however, the Lie algebras involved are of infinite dimension (they correspond to “curved Grassmannians” of submanifolds, rather than linear subspaces) and feasibility to obtain the final result becomes very doubtful.

The waves in plasma is another example of a nonholonomic system with infinite dimensional algebras involved. This is a pity, since the curvature tensor is known to be responsible for stability of the dynamical system expressed by geodesics; so, in principle, we could just have computed the N-curvature in this case and settle an intriguing question: “is cold fusion possible or not?”.

## §1. Structure functions

Let $`M`$ be a manifold of dimension $`n`$ over a field $`𝕂`$. Let $`F(M)`$ be the frame bundle over $`M`$, i.e., the principal $`GL(n;𝕂)`$-bundle. Let $`GGL(n;𝕂)`$ be a Lie group. A $`G`$-structure on $`M`$ is a reduction of the principal $`GL(n;𝕂)`$-bundle to the principal $`G`$-bundle. Another formulation seems to be clearer: if one can select transition functions from one coordinate patch to another so that they belong to $`G`$ for every intersection pair of patches, then we have a $`G`$-structure on $`M`$.

The simplest $`G`$-structure is the flat $`G`$-structure defined as follows. Let $`V`$ be $`𝕂^n`$ with a fixed frame. Fix a frame in $`V`$. Having identified $`V`$ with $`T_vV`$ via parallel translation along $`v`$, we get a fixed frame in every $`T_vV`$. The flat $`G`$-structure is the bundle over $`V`$ whose fiber over $`vV`$ consists of all frames obtained from the fixed one under the $`G`$-action.

In textbooks on differential geometry (e.g., in \[St\]) the obstructions to identification of the $`k`$-th infinitesimal neighborhood of a point $`mM`$ on a manifold $`M`$ with $`G`$-structure with that of a point of the manifold $`V`$ with the above flat $`G`$-structure are called structure functions of order $`k`$. It is also explained there that such an identification is possible provided all structure functions of lesser orders vanish.

To formulate the precise statement, we need some definitions. Let $`S^i`$ denote the operator of the $`i`$-th symmetric power. Set

$$𝔤_1=T_mM,𝔤_0=𝔤=\mathrm{Lie}(G);$$

for any $`i>0`$ and $`k=1`$, …, $`i`$ set:

$$𝔤_i=\{X\mathrm{Hom}(𝔤_1,𝔤_{i1})X(v_1)(v_2,v_3,\mathrm{},v_i)=X(v_2)(v_1,v_3,\mathrm{},v_i)\text{ for any }v_k𝔤_1\}.$$

Now, set $`(𝔤_1,𝔤_0)_{}=\underset{i1}{}𝔤_i`$.

Suppose that the $`𝔤_0`$-module $`𝔤_1`$ is faithful. Then, clearly,

$$(𝔤_1,𝔤_0)_{}𝔳𝔢𝔠𝔱(n)=𝔡𝔢𝔯𝕂[[x_1,\mathrm{},x_n]],\text{ where}n=dim𝔤_1.$$

It is subject to an easy verification that the Lie algebra structure on $`𝔳𝔢𝔠𝔱(n)`$ induces same on $`(𝔤_1,𝔤_0)_{}`$. The Lie algebra $`(𝔤_1,𝔤_0)_{}`$ will be called the Cartan’s prolong (the result of Cartan’s prolongation) of the pair $`(𝔤_1,𝔤_0)`$.

Let $`\mathrm{\Lambda }^i`$ be the operator of the $`i`$-th exterior power; set $`C_{(𝔤_1,𝔤_0)_{}}^{k,s}=𝔤_{ks}\mathrm{\Lambda }^s(𝔤_1^{})`$; we usually drop the subscript of $`C_{(𝔤_1,𝔤_0)_{}}^{k,s}`$ or indicate only $`𝔤_0`$. Define the differential $`_s:C^{k,s}C^{k,s+1}`$ setting for any $`v_1,\mathrm{},v_{s+1}V`$ (as usual, the slot with the hatted variable is ignored):

$$(\delta _sf)(v_1,\mathrm{},v_{s+1})=\underset{i}{}(1)^i[f(v_1,\mathrm{},\widehat{v}_{s+1i},\mathrm{},v_{s+1}),v_{s+1i}].$$

As expected, $`\delta _s\delta _{s+1}=0`$. The homology of this bicomplex is called Spencer cohomology of the pair $`(𝔤_1,𝔤_0)`$ and denoted by $`H_{(𝔤_1,𝔤_0)_{}}^{k,s}`$.

###### Proposition .

(\[St\]) Structure functions constitute the space of the $`(k,2)`$-th Spencer cohomology of the Cartan prolongation of the pair $`(T_mM,𝔤)`$ where $`𝔤=\mathrm{Lie}(G)`$.

The Riemann tensor on an $`n`$-dimensional Riemannian manifold is the structure function from $`H_{(𝔤_1,𝔬(n))_{}}^{k,2}`$.

###### Remark .

(cf. \[G\].)

$$\underset{k}{}H_{(𝔤_1,𝔤_0)_{}}^{k,s}=H^s(𝔤_1;(𝔤_1,𝔤_0)_{}).$$

This remark considerably simplifies calculations, in particular, if the Lie algebra $`(𝔤_1,𝔤_0)_{}`$ is simple and finite dimensional: we can apply the Borel-Weil-Bott- …(BWB) theorem. In the nonholonomic case considered in what follows the remark provides with a compact definition111Cf. with the problems encountered in the pioneer paper \[T\], where the case $`d=2`$ for Lie algebras was independently introduced. of structure functions. We can recover the bigrading (explicitly preserved in \[T\]) at any moment but to work with just one grading is much simpler, as with generating functions of series instead of individual terms. Moreover, there are a few general theorems that help to calculate Lie algebra cohomology and no theorems on Spencer cohomology.

## §2. Structure functions of nonholonomic manifolds

To embrace contact-like structures we have to slightly generalize the notion of Cartan prolongation: the tangent space to a point of a nonholonomic manifold is naturally rigged with a nilpotent Lie algebra structure (e.g., for any odd dimensional contact manifolds this is the Heisenberg algebra), cf. \[VG\].

### 2.1. Nonholonomic manifolds

(\[VG1\], \[VG2\]) Let $`M^n`$ be an $`n`$-dimensional manifold with a nonintegrable distribution $`D`$. Let

$$D=D_1D_2D_3\mathrm{}D_d$$

be the sequence of strict inclusions, where the fiber of $`D_i`$ at a point $`xM`$ is

$$D_{i1}|_x[D_1,D_{i1}]|_x$$

(here $`[D_1,D_{i1}]=\{[X,Y]XD_1,YD_{i1}\}`$) and $`d`$ is the least number such that

$$D_d|_x[D_1,D_d]|_x=D_d|_x.$$

In case $`D_d=TM`$ the distribution is called completely nonholonomic. The number $`d=d(M)`$ is called the nonholonomicity degree. A manifold $`M`$ with a distribution $`D`$ on it will be referred to as nonholonomic one if $`d(M)1`$.

Here we will consider nonholonomic manifolds only when they are completely nonholonomic. This does not diminish generality, actually, since if the distribution in question is not completely nonholonomic, we can take the integral manifold corresponding to it; the restriction of the same distribution onto this integral manifold is completely nonholonomic, all right. If, nevertheless, we have to, for some reason, consider the general case it is easy mutatis mutandis; it corresponds to nontransitive Shchepochkina prolongs, cf. below, and the general machinery is easily adjusted to them).

Let $`n_i(x)=dimD_i|_x`$. The distribution $`D`$ is called a regular one if all the dimensions $`n_i`$ are constants on $`M`$. In what follows we will only consider regular distributions. This is also not a restriction if everything is smooth.

The tangent bundle over a nonholonomic manifold $`(M,D)`$ can be naturally endowed with a bundle of nilpotent Lie algebras as follows. Let $`L`$ be the Lie algebra of vector fields preserving the distribution in a neighborhood of some point $`xM`$. Naturally, $`L`$ is filtered and let $`𝔤`$ be the associated graded Lie algebra. Set

$$𝔤_{}=\underset{di1}{}𝔤_i.$$

Clearly, $`d`$, the depth of $`𝔤_{}`$, coinsides with the nonholonomicity degree of $`D`$, see \[VG2\].

An analog of $`G`$-structure for nonholonomic manifolds is defined to be the pair $`(𝔤_{},𝔤_0)`$, where $`𝔤_0`$ is a subalgebra of the Lie subalgebra which preserves the $``$-grading of $`𝔤_{}`$, i.e., $`𝔤_0(𝔡𝔢𝔯𝔤_{})_0`$. If $`𝔤_0`$ is not indicated explicitely we assume that $`𝔤_0=(𝔡𝔢𝔯𝔤_{})_0`$, i.e., is the largest one.

Thus, we see that with every nonholonomic manifold $`(M,D)`$ a natural $`G`$-structure is associated: $`\mathrm{Lie}(G)=(𝔡𝔢𝔯𝔤_{})_0`$. But the structure functions of this $`G`$-structure do not reflect the nonholonomic nature of $`M`$. Indeed, when we compute the structure functions (be they “traditional” $`H_{(𝔤_1,𝔤_0)_{}}^{k,2}`$, or introduced above $`H^s(𝔤_1;(𝔤_1,𝔤_0)_{})`$) we use the fact that $`𝔤_1`$ is a commutative Lie algebra (partial derivatives commute).

To take the nonholonomic nature of $`M`$ into account we need an analog of the above Proposition for the case when the natural basis of the tangent space consists not of partial derivatives but rather of covariant derivatives corresponding to the connection determined by the same Pfaff forms whose vanishing determines the distribution. Therefore, instead of $`𝔤_1`$ we have a nilpotent algebra $`𝔤_{}`$.

So we have to define

(1) the simplest nonholonomic structure — the “flat” one, and

(2) the analog of Cartan prolong when $`𝔤_1`$ is replaced with $`𝔤_{}`$,

(3) (we have already understood) what is the analog of $`𝔤_0`$ if it is not mentioned as is usually the case (only the distribution is given), but it remains to figure out

(4) what is the analog of $`(𝔤_1,𝔤_0)_{})`$, and, finally,

(5) what is the analog of $`H_{(𝔤_1,𝔤_0)_{}}^{k,2}`$.

Let us answer these remaining questions.

Fix a frame $`F`$ in $`𝕂^n`$. Having identified $`T_vV`$ with $`V`$ by means of the transformation that preserves $`D`$, we fix a frame — the image of $`F`$ — in each $`T_vV`$. The nonholonomic structure of the nonholonomic manifold $`(𝕂^n,D)`$ is the pair of bundles: the frame bundle and the distribution $`D`$. The nonholonomic structure of $`(𝕂^n,D)`$ will be considered to be flat if the frames over $`v`$ are obtained from the fixed one by means of the $`G`$-action, where $`G`$ is a group whose Lie algebra is $`(𝔡𝔢𝔯𝔤_{})_0`$.

### 2.2. Structure functions of a nonholonomic manifold

Given a pair $`(𝔤_{},𝔤_0)`$ as above, define its $`k`$-th Shchepochkina prolong (first published by I. Shchepochkina \[Sh\] for Lie superalgebras; see also more accessible \[Sh5\]) for $`k>0`$ to be:

$$𝔤_k=(i(S^{}(𝔤_{})^{}𝔤_0)j(S^{}(𝔤_{})^{}𝔤_{}))_k,$$

where the subscript singles out the component of degree $`k`$ and where $`i:S^{k+1}(𝔤_1)^{}𝔤_1S^k(𝔤_1)^{}𝔤_1^{}𝔤_1`$ is a natural embedding and $`j:S^k(𝔤_1)^{}𝔤_0S^k(𝔤_1)^{}𝔤_1^{}𝔤_1`$ a natural map.

Similarly to the case when $`𝔤_{}`$ is commutative, define $`(𝔤_{},𝔤_0)_{}`$ to be $`\underset{kd}{}𝔤_k`$; then, as is easy to verify, $`(𝔤_{},𝔤_0)_{}`$ is a Lie algebra. By the same reasons as for the manifolds with $`G`$-structure (\[St\]), the space $`H^2(𝔤_{};(𝔤_{},𝔤_0)_{})`$ is the space of structure functions — obstructions to the identification of the infinitesimal neighborhood of a point on the manifold with a nonholonomic structure (given by $`𝔤_{}`$ and $`𝔤_0`$) with same on a flat nonholonomic manifold with the same $`𝔤_{}`$ and $`𝔤_0`$.

The space $`H^2(𝔤_{};(𝔤_{},𝔤_0)_{})`$ naturally splits into homogeneous components whose degree will be called the order of the structure function; the minimal order of structure function is equal to $`2d`$. As in the case of a commutative $`𝔤_{}=𝔤_1`$, the structure functions of order $`k`$ can be interpreted as obstructions to flatness of the nonholonomic manifold with the $`(𝔤_{},𝔤_0)`$-structure provided the obstructions of lesser orders vanish. Observe that for a nonholonomic manifold the order of structure functions is no more in direct relation with the numbers of the infinitesimal neighborhoods of the points we wish to identify: distinct partial derivatives bear different “weights”.

When I told this to A. Beilinson in 1989 he winced and, as I got this far, he said that though he conceeded the idea is of interest, he doubted it can work in the general case: all my constructions are graded, whereas with every distribution only filtered algebras $`𝔤_{}`$ are naturally associated. This incredulity is justified but fortunately, different filtrations on $`𝔤_{}`$ with the same graded $`𝔤_{}`$ are governed precisely by the structure functions of lesser orders and all of them must vanish for the N-curvature to be well-defined. So, given a filtered $`𝔤_{}`$, just take the associated graded algebra and proceed as indicated above. This indicates a way to study the stable normal forms of Pfaff equations listed in \[Z\].

## §3. The Lie algebra $`𝔢`$ of Engel vector fields

Let us describe the Lie algebra $`𝔢`$ of vector fields that preserve the Engel distribution. Since each of the forms $`\alpha _1`$ and $`\alpha _2`$ is a familiar contact form, we know enough to demonstrate (we leave it as an easy exercise for the reader) the following statement.

###### Lemma .

A weight basis of the $``$-graded Lie algebra $`𝔢`$ is as follows (indicated are the degrees):

| $`3`$ | $`2`$ | $`1`$ | $`0`$ | $`n1`$ |
| --- | --- | --- | --- | --- |
| $`_4`$ | $`_3+x_1_4`$ | $`_1`$ | $`x_1_1+x_2_2+2x_3_3+3x_4_4`$ |  |
|  |  | $`_2+x_1_3+\frac{1}{2}x_1^2_4`$ | $`x_1_1x_2_2+x_4_4`$ | $`x_1^{n+1}_2+\frac{x_1^{n+2}}{n+2}_3+\frac{x_1^{n+3}}{n+3}_4`$ |
|  |  |  | $`x_1_2+\frac{x_1^2}{2}_3+\frac{x_1^3}{3}_4`$ |  |

The structure of the Engel algebra $`𝔢`$ seem to have never been described; so let us do it. First, introduce a shorthand for some elements of $`𝔢`$:

$$\begin{array}{cc}\hfill D_3& =_3+x_1_4;\hfill \\ \hfill E& =x_1_1+x_2_2+2x_3_3+3x_4_4;\hfill \\ \hfill H& =x_1_1x_2_2+x_4_4;\hfill \\ \hfill X_n& =x_1^{n+1}_2+\frac{x_1^{n+2}}{n+2}_3+\frac{x_1^{n+3}}{n+3}_4\text{ for }n1.\hfill \end{array}$$

Now, it is clear that $`𝔢`$ is solvable with

$$𝔢/[𝔢,𝔢]=Span(E,H);𝔢^{(1)}/[𝔢^{(1)},𝔢^{(1)}]=Span(_1),\text{ where }𝔢^{(1)}=[𝔢,𝔢].$$

The negative part, $`𝔢_{}`$ plays for $`𝔢`$ the same role the Heisenberg algebra plays for the contact algebra. This negative part was, of course calculated earlier, see, e.g., \[G1\] where it also bears the name of Engel algebra; perhaps one should say the “small” and the “total” ones. It seems very strange that the total algebra had never beed desribed; I am sure such a description is buried in works of classics.

In review to \[G1\] \[MR 96k:58010\] it is written “There is a Darboux theorem for Engel structures: any two are locally diffeomorphic”. Vershik also quoted this statement as a well-known. But I can not place it. Indeed, by Darboux theorem for the contact structure I understand the statement that says that the contact form can be locally reduced to a canonical form , cf. \[Z\]. (In contradistinction, a related statement that the nondegenerate differential 2-form $`\omega `$ can only be reduced to a canonical form if and only if $`d\omega =0`$, cf. \[P\], so the statement on reducibility of the nondegenerate closed 2-form, sometimes encountered in textbooks on Differential Geometry, is a tautology.)

Concerning Engel structures, there is, for example, paper \[G2\] with examples of several non-equivalent Engel structures on $`^4`$. Doesn’t it contradict the statement on an analog of Darboux theorem? So let us calculate the N-curvature of the contact and Engel structures: to find out if these invariants vanish, hence, there are analogs of Darboux theorem for them.

## §4. Structure functions of the contact and Engel structures

Éli Cartan worked before “(co)homology” was defined, still, he introduced a notion that helps to calculate some of these cohomologies. This notion is involutive Lie algebra of vector fields. Serre gave a nice reformulation of this notion (\[St\]); thanks to Serre’s statement for the case when $`(𝔤_{},𝔤_0)_{}`$ is a simple vectorial Lie algebra instead of computing $`H^2(𝔤_1;(𝔤_1,𝔤_0)_{})=\underset{k}{}H_{(𝔤_1,𝔤_0)_{}}^{k,2}`$, it only suffices to compute the summand with the least $`k=1`$: only this term might be nonzero because $`(𝔤_{},𝔤_0)_{}`$ is involutive. This really helps, see examples in \[LPS\].

Regrettably, Serre’s theorem is only formulated so far for algebras of depth 1. So I calculated $`H^2(𝔤_{};(𝔤_{},𝔤_0)_{})`$ by hands; later the result was confirmed by independent computer calculations up to convincingly high degree:

###### Theorem .

The space of nonholonomic structure functions on any Engel manifold is of dimension $`2`$. All structure functions are of order $`2`$, their representatives are:

$$_1^{}X_1^{}H+_1^{}X_1^{}E\text{ and}_4^{}D_3^{}_4.$$

The first of these cycles is “pure” in the sence that it does not represent its cohomology class modulo any coboundary, the coboundaries vanish. The other one is mixed: one can add any coboundary to it.

The following is another form of the classical theorem (see \[LP\]):

###### Theorem .

(Darboux theorem) The space of nonholonomic structure functions on any contact manifold vanishes.


format: text

text:
# Effects of Dynamical Quarks in UKQCD Simulations

## 1 Introduction

The UKQCD Collaboration has embarked upon a study of unquenching effects in QCD lattice simulations particularly in the static quark potential, light hadron spectrum and topological sector. The philosophy that UKQCD uses in this study is to simulate at a variety of sea quark masses, $`m^{sea}`$, but with a fixed the lattice spacing, $`a`$. This is in order to maintain, as far as possible, constant lattice systematic effects due to the finite-ness of $`a`$ and the volume.

It is well known that the lattice cut-off, $`a`$, is a function of both the gauge coupling, $`\beta `$, and the dynamical quark mass, $`m^{sea}`$ (see e.g. ). For this reason, the philosophy we have chosen is to simulate at points in the $`(\beta ,m^{sea})`$ plane which have constant $`a`$. We term this the “matched” trajectory. The hope is that this procedure will disentangle lattice systematic effects (due to the finite-ness of $`a`$ and volume) with unquenching effects. Any variation of a physical quantity along this matched trajectory can be attributed to unquenching effects rather than, e.g. $`𝒪(a^2)`$ effects.

This paper presents an overview of this work. See  for a full description.

## 2 Simulation Details

The standard Wilson gauge action was used together with the Clover $`𝒪(a)`$-improved fermion action. The coefficient, $`c_{SW}`$, used in the Clover term was non-perturbatively determined by the Alpha Collaboration .

A range of $`(\beta ,\kappa ^{\mathrm{sea}})`$ values were chosen in order to maintain a constant value of the lattice spacing. This used the technology of . The parameter values for the simulations are displayed in Table 1. The lattice volume used was $`16^3\times 32`$ throughout.

Table 1 lists also the lattice spacing obtained from the Sommer scale, $`r_0`$ . The first error listed is statistical, and the second (shown as $`{}_{y}{}^{}{}_{}{}^{+x}`$) is the systematic error from variations in the fit used for the static quark potential. The central values quoted were obtained using all potential data satisfying $`\sqrt{2}𝒓8`$. It can be seen that the last four simulations listed (i.e. those at $`\kappa ^{\mathrm{sea}}=0.1350,0.1345,0.1340`$ and $`0`$) are matched to within errors, with the $`\kappa ^{\mathrm{sea}}=0`$ simulation corresponding to a quenched run.

The simulations at $`\kappa ^{\mathrm{sea}}=0.1355`$ and $`0.13565`$ were performed in order to study lighter sea quark masses and do not lie on the above matched trajectory.

Full details of autocorrelation times and other algorithmic issues can be found in .

## 3 The static potential

The static quark potential was determined for all our ensembles using the method originally proposed in  and using a variational basis of “fuzzed” links .

Figure 1 plots the static potential in units of $`r_0`$. The zero of the potential has been set at $`r=r_0`$. The data are well described by the universal bosonic string model potential  which predicts

$`\left[V(r)V(r_0)\right]r_0`$ $`=`$ $`(1.65e)\left({\displaystyle \frac{r}{r_0}}1\right)`$ (1)
$``$ $`e\left({\displaystyle \frac{r_0}{r}}1\right).`$
Of course, the fact that the scaled potential measurements all have the same value and slope at $`r=r_0`$ simply reflects the definition of $`r_0`$. Figure 2 shows the deviations from this model potential for the matched ensembles. Here $`e=\pi /12`$ .

Overall there is no discernable difference between the ensembles at distance $`rr_0`$. Furthermore, the data follow the string model very well. However, at shorter distances, $`r<0.5r_0`$, there is a deviation from the string model, and a variation amongst the ensembles. This seems to be systematic in the sea quark mass – the deviation from the string model increases as $`m^{sea}`$ decreases. In fact, careful correlated fits of the potential show that the parameter $`e`$ in Eq.(1) increases by 18$`{}_{10}{}^{}{}_{}{}^{+13}`$% in going from the quenched to $`\kappa ^{\mathrm{sea}}=.13500`$ data. We note that the ensembles in Figure 2 are matched, and therefore we can exclude lattice systematics from the effects that we have observed. (Similar findings in the case of two flavours of Wilson fermions have been reported by the SESAM-T$`\chi `$L collaboration  where an increase of $`1633\%`$ was found.) Finally we note that there is no evidence of string breaking in our static quark potential. However the distance scales covered $`r<\mathrm{\hspace{0.33em}1.3}fm`$ is not large.

## 4 Hadronic Spectrum

In this section one of our main aims will be to uncover unquenching effects in the light hadron spectrum. Because we have a matched data set, any variation amongst our ensembles can be attributed to unquenching effects. However, the task of identifying variations is likely to be hard for those quantities which are primarily sensitive to physics at the same scale as that used to define the matching trajectory in the $`(\beta ,\kappa ^{\mathrm{sea}})`$ parameter space ($`r_0`$ in this case). This is expected to be the case for the hadron spectrum considered here where the quark masses are still relatively heavy.

Two-point hadronic correlation functions were produced for each of the ensembles using interpolating operators for pseudoscalar, vector, nucleon and delta channels as described in . Mesonic correlators were constructed using both degenerate and non-degenerate valence quarks, whereas only degenerate valence quarks were used for the baryonic correlators.

In the following, we review the main fitting procedures which were used to obtain the light hadron spectrum results. Further details of the fitting procedure can be found in .

We used the fuzzing procedures of  to generate correlators of the type LL, FL and FF where F denotes fuzzed, and L local operators.

Correlated fits were used throughout the fitting analysis of the correlation functions with the eigenvalue smoothing technique of of  employed. A factorising fit was performed which combined the three fuzzed correlators LL, FL and FF together to extract a better estimate of the ground state parameters.

Full details of the fitting procedure can be found in .

### 4.1 The J parameter

In Figures 3 and 4 the vector meson masses and hyperfine splittings are plotted against the corresponding pseudoscalar masses for all the datasets. It is difficult to identify an unquenching signal from these plots - the data seem to overlay each other. Note that in , it was reported that there was a tendency for the vector mass to increase as the sea quark mass decreases (for fixed pseudoscalar mass). The observations for the present matched dataset imply that this may have been due to either an $`𝒪(a)`$ effect (since the dataset in  was not fully improved at this level) or a finite volume effect. The conclusion therefore is that it is important to run at a fixed $`a`$ in order to disentangle unquenching effects from lattice artifacts or finite volume effects.

A possible explanation for the lack of convincing unquenching effects in our meson spectrum is the following. Our matched ensembles are defined to have a common $`r_0`$ value, so any physical quantity that is sensitive to this distance scale (and the static quark potential itself) will also, by definition, be matched. Our mesons, because they are composed of relatively heavy quarks, are examples of such quantities, and this is a possible reason why there is no significant evidence of unquenching effects in the meson spectrum.

A further point regarding hyperfine splitting in Figure 4 is that the lattice data for the matched ensembles tends to flatten as the sea quark mass decreases. (The quenched data has a distinctly negative slope, whereas the $`\kappa ^{\mathrm{sea}}=0.1350`$ data is flat.) Thus the lattice data is tending towards the same behaviour as the experimental data which lies on a line with positive slope (independent of the value used for $`r_0`$). This behaviour is apparently spoiled by the unmatched run with $`\kappa ^{\mathrm{sea}}=0.1355`$ (see Figure 4) which has a clear negative slope. However, the $`\kappa ^{\mathrm{sea}}=0.1355`$ data does not satisfy the finite volume bound of , and therefore we can attribute the trend in this data to a lattice artefact.

We now study the $`J`$-parameter defined as

$$J=M_V\frac{dM_V}{dM_{PS}^2}|_{K,K^{}}.$$
(2)
In the context of dynamical fermion simulations, this parameter can be calculated in two ways. The first is to define a partially quenched $`J`$ for each value of the sea quark mass. In this case, the derivative in (2) is with respect to variations in the valence quark mass (with the sea quark mass fixed). The second approach is to define $`J`$ along what we will term the ‘unitary’ trajectory, i.e. along $`\kappa ^{\mathrm{sea}}=\kappa ^{\mathrm{val}}`$. In Table 2, the results from both methods are given. These values of $`J`$ are around 25% lower than the experimental value $`J_{expt}=0.48(2)`$.

Finally we note that the physical value of $`J`$ (i.e. that which most closely follows the procedure used to determine the experimental value of $`J_{expt}=0.48(2)`$), should be obtained from extrapolating the results from the first approach to the physical sea quark masses. We call this the third approach. In order to perform this extrapolation, we extrapolate the three matched dynamical $`J`$ values obtained from the first approach linearly in $`(M_{PS}^{\mathrm{unitary}})^2`$ to $`(M_{PS}^{\mathrm{unitary}})^2=0`$. $`M_{PS}^{\mathrm{unitary}}`$ is the pseudoscalar meson mass at the unitary point. The value for $`J`$ from the third approach is presented in Table 2 and we note that it is approaching the experimental value for $`J`$.

The results from all three approaches are plotted in Figure 5, together with the experimental result. There is some promising evidence that the lattice estimate of $`J`$ increases towards the experimental point as the sea quark mass decreases (see the $`J`$ value from approaches 1 and 3).

Recently there has been a proposed ansatz for the functional form of $`M_V`$ as a function of $`M_{PS}^2`$ . Although this ansatz would be interesting to pursue, all our data have $`M_{PS}/M_V>\mathrm{\hspace{0.33em}0.6}`$, and for this region, the ansatz of  is linear to high precision. Therefore we have chosen to interpolate our data with a simple linear function and await more chiral data before using the ansatz of .

### 4.2 Lattice spacing

This subsection presents a determination of $`a`$ from the meson spectrum which complements that from $`r_0`$.

A common method of determining $`a`$ from the meson spectrum uses the $`\rho `$ mass. However, this requires the chiral extrapolation of the vector meson mass down to (almost) the chiral limit which, as was discussed in the previous subsection, may be problematic. An alternative method of extracting the lattice spacing using the vector meson mass at the simulated data points (i.e. without any chiral extrapolation) was given in . Using this method, we obtain the lattice spacing values as shown in Table 1 labelled $`a_J`$. Note that these are in general 10-15% larger than the values from Section 3 where the lattice spacing was determined from $`r_0`$. A possible explanation for this discrepancy is that the potential and mesonic spectrum are contaminated with different $`𝒪(a^2)`$ errors (or that the value $`r_0=0.49`$ fm is 10-15% too small!).

In order to investigate unquenching effects in the meson spectrum, we define the quantity

$$\delta _{i,j}(\beta ,m^{sea})=1\frac{a_i(\beta ,m^{sea})}{a_j(\beta ,m^{sea})},$$
(3)
where $`a_i`$ is the lattice spacing determined from the physical quantity $`i=\{M_\rho ,M_K,f_\pi \mathrm{}\}`$. Note that when $`\delta _{i,j}=0`$, the lattice prediction of $`M_i`$ with scale taken from $`M_j`$ agrees with experiment. Thus $`\delta `$ is a good parameter to study unquenching effects. We expect that $`\delta _{i,j}(\beta ,m^{sea})=𝒪(a^2)`$ since we are working with a non-perturbatively improved clover action.

In Figure 6, $`\delta _{i,j}`$ is plotted against $`(aM_{PS}^{\mathrm{unitary}})^2`$ for the matched datasets. In this plot we have fixed $`j=r_0`$ and the various physical quantities $`i`$ are $`\sqrt{\sigma }`$ (the string tension) and the hadronic mass pairs $`(M_K,M_K)`$ & $`(M_\rho ,M_\pi )`$. The method used to determine the scale $`a_i`$ from these mass pairs is that of . It is worth noting that the experimental point on this same plot would occur at an $`x`$co-ordinate (depressingly) of $`(aM_{PS}^{\mathrm{unitary}})^2=(aM_\pi )^2200`$.

Figure 6 does not show signs of unquenching for quantities involving the hadronic spectrum i.e. the mass pairs $`(M_K,M_K)`$ & $`(M_\rho ,M_\pi )`$. However, there is evidence of unquenching effects when comparing the scale from $`r_0`$ with that from $`\sqrt{\sigma }`$ since the quenched value of $`\delta _\sqrt{\sigma }`$ is distinct from the dynamical values.

One may wonder if the the $`\delta `$ values may have been distorted by not choosing the simulation parameters $`(\beta ,m^{sea})`$ exactly on the matched trajectory. In order to obtain a rough estimate of the effect of mis-matched value of $`\beta `$, we use the renormalisation group inspired ansatz for $`a_i`$ :

$$a_i(g_0^2)=\mathrm{\Lambda }^1f_{PT}(g_0^2)\times \left[1+X_if_{PT}(g_0^2)^{n_i}\right],$$
(4)
where $`f_{PT}(g^2)`$ is the usual asymptotic scaling function obtained from integrating the $`\beta `$-function of QCD and $`X_i`$ is the coefficient of the $`𝒪(a^n)`$ lattice systematic. The functional form for $`a(g_0^2)`$ was originally applied for the quenched theory, but let us assume that it can also be applied in the unquenched case. Using Eq.(4), we see that a mismatch in $`\beta `$ of $`\mathrm{\Delta }\beta `$ would lead to a relative error in $`\delta `$ of

$$\frac{\delta (\beta +\mathrm{\Delta }\beta )\delta (\beta )}{\delta (\beta )}3\mathrm{\Delta }\beta .$$
(5)
This shows that even an error in $`\beta `$ of as much as $`\mathrm{\Delta }\beta 0.01`$ introduces a relative error in $`\delta (\beta )`$ of only 3%, ruling out the possibility that a significant distortion in $`\delta `$ could have occurred due to a mismatching in $`\beta `$.

### 4.3 Edinburgh Plot

In Table 1 the ratios $`M_{PS}/M_V`$ are displayed for the unitary case $`\kappa ^{\mathrm{sea}}`$= $`\kappa ^{\mathrm{val}}`$. As can be seen the simulated data is a long way from the experimental value $`M_\pi /M_\rho =0.18`$. Figure 7 shows the ‘Edinburgh plot’ ($`M_N/M_V`$ v.s. $`M_{PS}/M_V`$) for all the data sets. There is no significant variation within the dynamical data as the sea quark mass is changed, but the dynamical data does tend to lie above the (matched) quenched data. This latter feature may be indicative of finite volume effects since these are expected to be larger in full QCD compared to the quenched case .

## 5 Chiral Extrapolations

In  we used 3 approaches to perform the chiral extrapolations of hadron masses: “Pseudo-Quenched”; “Unitary Trajectory”; and a “Combined Chiral Fit”. In this publication we will refer only to the last approach. We take

$$\widehat{M}(\kappa ^{\mathrm{sea}};\kappa ^{\mathrm{val}})$$

$$=A(\kappa ^{\mathrm{sea}})+B(\kappa ^{\mathrm{sea}})\widehat{M}_{PS}(\kappa ^{\mathrm{sea}};\kappa ^{\mathrm{val}})^2$$

$$=A_0+A_1\widehat{M}_{PS}(\kappa ^{\mathrm{sea}};\kappa ^{\mathrm{sea}})^2$$

$$+\left[B_0+B_1\widehat{M}_{PS}(\kappa ^{\mathrm{sea}};\kappa ^{\mathrm{sea}})^2\right]\widehat{M}_{PS}(\kappa ^{\mathrm{sea}};\kappa ^{\mathrm{val}})^2,$$

using the nomenclature $`\widehat{M}aM`$. The first argument of $`M(\kappa ^{\mathrm{sea}};\kappa ^{\mathrm{val}})`$ refers to the sea quark and the second to the valence quark.

The results of these extrapolations are shown in Table 3. We stress that this functional form for the extrapolation is not motivated by theory, but is used as a numerical analysis technique in order to test for evidence of unquenching effects. As can be seen from Table 3, the parameters $`A_1`$ and $`B_1`$ are compatible with zero (to $`2\sigma `$) and therefore we conclude that there is no evidence of unquenching effects (i.e. there is no statistical variation in the fit parameters $`A`$ and $`B`$ with sea quark mass $`\kappa ^{\mathrm{sea}}`$).

## 6 Conclusions

This paper attempts to uncover unquenching effects in the dynamical lattice QCD simulations at a fixed (matched) lattice spacing (and volume) and various dynamical quark masses. This approach allows a more controlled study of unquenching effects without the possible entanglement of lattice and unquenching systematics.

We see some sign of unquenching effects in the static quark potential at short distance, but no significant sign of unquenching effects in the meson spectrum. This is presumably since our dynamical quarks are relatively massive, and so the meson spectrum is dominated by the static quark potential. This potential is, by definition, matched amongst our ensembles at the hadronic length scale $`r_0`$, and so any variation of the meson spectrum within our matched ensemble must surely be a “higher” order unquenching effect which is beyond our present statistics.

It is interesting to note that other work (using the same ensembles) has shown interesting unquenching effects in the glueball and topological sector .

## 7 Acknowledgements

The author wish to thank all of his collaborators in UKQCD. The support of the Particle Physics and Astronomy Research Council is gratefully acknowledged.


text:
# References

IFT preprint UFIFT-MATH-02-1

On a low energy bound in a class of chiral field theories with solitons

Sergei V. Shabanov 111on leave from Laboratory of Theoretical Physics, JINR, Dubna, Russia

Institute for Fundamental Theory, Departments of Physics and Mathematics,
University of Florida, Gainesville, FL- 32611, USA

## Abstract

A low energy bound in a class of chiral solitonic field theories related the infrared physics of the SU(N) Yang-Mills theory is established.

1. The model. Consider $`N1`$ smooth fields $`n_a=n_a(x)`$ in spacetime taking their values in the Lie algebra of $`SU(N)`$. The fields are chosen to be commutative $`[n_a,n_b]=0`$ and orthonormal $`(n_a,n_b)=\delta _{ab}`$ with respect the Cartan-Killing form in the Lie algebra. For any two Lie algebra elements $`\xi `$ and $`\eta `$, the Cartan-Killing form is defined as $`(\xi ,\eta )=\mathrm{tr}(\widehat{\xi }\widehat{\eta })`$ where the operator $`\widehat{\xi }`$ acts on the Lie algebra as a Lie derivative $`\widehat{\xi }\eta =[\xi ,\eta ]`$. There can only be $`N1`$ mutually commutative and linearly independent elements in the Lie algebra of $`SU(N)`$ because the rank of $`SU(N)`$ is $`r=N1`$ (the dimension of the Cartan subalgebra). If $`h_a`$ form an orthonormal basis in the Cartan subalgebra in a matrix representation of SU(N), then

$$n_a(x)=U^{}(x)h_aU(x),$$
(1)
where $`U(x)SU(N)`$. In Eq. (1) $`U(x)`$ is defined modulo the left multiplication by elements from the Cartan subgroup generated by $`h_a`$ (the maximal Abelian subgroup $`T=U(1)^{N1}`$). So, in fact, $`U(x)SU(N)/T`$ since any group element can be represented as a product of an element of $`T`$ and an element of the quotient $`SU(N)/T`$. Under the condition that $`n_a`$ approach fixed constant values at the spatial infinity, $`n_a(x)h_a`$, i.e., $`U(x)`$ approaches the group unity, the fields $`n_a`$ define a map of a spatial three-sphere $`𝐒^3`$ into the manifold $`SU(N)/T`$ for every moment of time. The third homotopy group of this map is nontrivial $`\pi _3(G/T)Z`$, $`G=SU(N)`$. When $`N=2`$, the only field $`n_1`$ can be regarded as a unit 3-vector. It is a Hopf map: $`𝐒^3𝐒^2SU(2)/U(1)`$. The corresponding topological number is the Hopf invariant which can also be interpreted as a linking number of two curves in $`𝐒^3`$ being preimages of two distinct points of $`𝐒^2`$. The two-forms $`F^a=F_{jk}^adx^jdx^k`$, $`j,k=1,2,3`$, where

$$F_{jk}^a=iN_b(n_a,[_jn_b,_kn_b]),$$
(2)
are closed, that is, $`F_{jk}^a=_jC_k^a_kC_j^a`$. This is proved at the end of next section. The forms $`F^a`$ may not be exact. This follows from the fact that the cohomology ring $`H^{}(G/T)`$ is rationally generated by $`H^2(G/T)`$ . The topological number of the map $`𝐒^3G/T`$ should be constructed out of the 2-forms $`F^a=iN(n_a,_b[dn_b,dn_b])`$ on $`G/T`$. Introducing the field $`B_i^a=\frac{1}{2}ϵ_{ijk}F_{jk}^a`$ with $`ϵ_{ijk}`$ being the Levi-Chevita tensor, the topological number of the above map can be written as

$$Q=(16\pi ^2N)^1𝑑x_{i,a}C_i^aB_i^a.$$
(3)
For SU(2), $`Q`$ is a Hopf invariant. Since $`SU(2)/U(1)G/T`$, the normalization coefficient in (3) can be chosen so that $`Q`$ is an integer when $`n_a`$ realize a Hopf map.

An explicit realization of the Hopf map by the fields $`n_a`$ is as follows. Consider the Cartan-Weyl basis in the Lie algebra. Let $`\alpha `$ be a positive root. For every positive root $`\alpha `$, there are two basis elements $`e_\alpha `$ and $`e_\alpha =\overline{e}_\alpha `$ such that for any element $`h`$ from the Cartan subalgebra

$`[h,e_\alpha ]`$ $`=`$ $`(h,\alpha )e_\alpha ,`$ (4)
$`[e_\alpha ,e_\alpha ]`$ $`=`$ $`\alpha ,[e_\alpha ,e_\beta ]=N_{\alpha ,\beta }e_{\alpha +\beta },`$ (5)
where $`N_{\alpha ,\beta }0`$ if $`\alpha +\beta `$ is a root. Note that the elements $`\alpha ,e_\alpha `$ and $`e_\alpha `$ form a basis of an SU(2) subalgebra (associated with the root $`\alpha `$). Let $`U_\omega (x)SU(2)/U(1)SU(N)/T`$ where the subgroup SU(2) is associated with a simple root $`\omega `$. One can always choose $`h_1=N^{1/2}\omega `$. The norm of any root of SU(N) is $`1/N`$ with respect to the Cartan-Killing form (see next section). Then $`n_1(x)=N^{1/2}U_\omega ^{}\omega U_\omega `$ is a Hopf map. The other fields realize a trivial map, $`n_a=h_a,a>1`$. Indeed, $`U_\omega (x)=\mathrm{exp}[iu_\omega (x)]`$ where $`u_\omega (x)=\phi _\omega (x)e_\omega +\overline{\phi }_\omega (x)e_\omega `$. For $`a>1`$, it follows from (4) that $`n_a=U_\omega ^{}h_aU_\omega =h_a`$ because $`(\omega ,h_a)(h_1,h_a)=0`$. Now, if we introduce an orthonormal basis in the SU(2) subgroup, $`\tau _1=i(e_\omega e_\omega )/\sqrt{2}`$, $`\tau _2=(e_\omega +e_\omega )/\sqrt{2}`$ and $`\tau _3=\sqrt{N}\omega `$, then $`[\tau _j,\tau _k]=iN^{1/2}ϵ_{jkn}\tau _n`$. Let $`𝒏`$ be an isotopic unit three-vector whose components are $`(\tau _j,n_1)`$. It defines the Hopf map by construction. From (2) we infer that $`F_{jk}^a=\delta ^{a1}\sqrt{N}𝒏(_j𝒏\times _k𝒏)`$. Hence our $`B_j^a`$ and $`C_j^a`$ contain an extra factor $`\sqrt{N}`$ when the fields $`n_a`$ realize a Hopf map associated with an SU(2) subgroup of SU(N). This explains the normalization factor $`N^1`$ in (3). Since all the root have the same norm in SU(N), the normalization coefficient in (3) for any SU(2) subgroup has to be the same. The root system is invariant under the Weyl symmetry, and so should be $`Q`$. The sum over $`a`$ in (3) provides this invariance.

The dynamics of the fields $`n_a`$ is determined by the Lagrangian density

$$=m^2\underset{\mu ,a}{}(_\mu n_a,_\mu n_a)\frac{g}{4}\underset{\mu ,\nu ,a}{}F_{\mu \nu }^aF_{\mu \nu }^a,$$
(6)
$`\mu ,\nu =0,1,2,3`$; and $`_0`$ stands for the time derivative. In the case of SU(2), this Lagrangian density describes the Faddeev model  for knot solitons. The knot solitons have been extensively studied numerically . The model (6) has been introduced in  and may also have solitonic solutions. The Lagrangian (6) is believed to describe (in a certain approximation) the infrared physics of the SU(N) Yang-Mills theory . Recent analytical  and lattice  studies of this correspondence in the SU(2) case look promising.

Due to the Lorentz symmetry of the Langrangian density, a Lorentz transformation of a static solution is a time dependent solution of the Euler-Lagrange equations for (6). Solutions that describe interacting solitons are not static (even modulo Lorentz transformations). In this paper a low energy bound for static solitons with a topological number $`Q`$ is established:

$`E`$ $``$ $`c_N|Q|^{3/4},`$ (7)
$`c_N`$ $`=`$ $`8\pi ^23^{3/8}\left({\displaystyle \frac{2N^3}{N^2N1}}\right)^{1/4}\sqrt{m^2g},`$ (8)
where $`E`$ is the energy functional

$`E`$ $`=`$ $`m^2{\displaystyle 𝑑x\underset{j,a}{}(_jn_a,_jn_a)}+{\displaystyle \frac{g}{2}}{\displaystyle 𝑑x\underset{j,a}{}B_j^aB_j^a}`$ (9)
$``$ $`{\displaystyle 𝑑x\left(_2(x)+_4(x)\right)}E_2+E_4.`$ (10)
For the Faddeev-Hopf knot solitons the low energy bound was found in  and improved in  (meaning a larger constant $`c_2`$). Beyond conventional perturbation theory, the Yang-Mills quantum dynamics can be studied by the large $`N`$ expansion method with the purpose to establish a relation (duality) to a string theory on some manifold. Therefore it is of interest to investigate the $`N`$ dependence of the low energy bound for solitons in the model (6).

2. Notations and necessary facts. We would need the following algebraic inequalities. Let $`a_i,b_i0`$, $`a=_ia_i`$, $`b=_ib_i`$ and $`\gamma 1`$. Then

$`a_1^\gamma +a_2^\gamma +\mathrm{}+a_r^\gamma `$ $``$ $`a^\gamma ;`$ (11)
$`\sqrt{a}\sqrt{a_1}+\sqrt{a_2}+\mathrm{}+\sqrt{a_r}`$ $``$ $`\sqrt{r}\sqrt{a};`$ (12)
$`_ia_i^pb_i^q`$ $``$ $`a^pb^q,p+q=1.`$ (13)
Define $`p_i=a_i/a1`$. Then $`_ip_i=1`$. The inequality (11) follows from an obvious inequality $`p_i^\gamma p_i`$ if one takes the sum over $`i`$. The second inequality is proved by squaring it and applying the basic algebraic inequality $`\sqrt{a_i}\sqrt{a_j}\frac{1}{2}(a_i+a_j)`$. The third inequality is an algebraic Hölder inequality (see, e.g., ).

An arrow is used to denote vectors in space, e.g., $`\stackrel{}{}\varphi =(_1\varphi ,_2\varphi ,_3\varphi )`$ for the gradient. The scalar product for two vector fields is

$$\stackrel{}{u},\stackrel{}{v}=𝑑x\stackrel{}{u}\stackrel{}{v}.$$

The $`L_p`$ norm of a vector field reads

$$\stackrel{}{u}_p=\left[𝑑x\left(\stackrel{}{u}\stackrel{}{u}\right)^{p/2}\right]^{1/p}.$$
(14)
The following functional inequalities are used in the sequel

$`|\stackrel{}{u},\stackrel{}{v}|`$ $``$ $`\stackrel{}{u}_p\stackrel{}{v}_q,p^1+q^1=1,`$ (15)
$`\stackrel{}{u}_{6/5}`$ $``$ $`\stackrel{}{u}_1^{2/3}\stackrel{}{u}_2^{1/3},`$ (16)
$`\stackrel{}{u}_6`$ $``$ $`\lambda _1\mathrm{𝖼𝗎𝗋𝗅}\stackrel{}{u}_2,\lambda _1=(48)^{1/6}(3\pi )^{2/3}.`$ (17)
The first two inequalities are Hölder type inequalities . The third one follows from Rosen’s result for scalar functions  (cf. also )

$$\varphi _6\lambda _1\stackrel{}{}\varphi _2,$$
(18)
where the $`L_p`$ norm for scalar functions is defined by (14) for one-dimensional vectors. Let $`\varphi =(\stackrel{}{u}\stackrel{}{u})^{1/2}`$. We have

$$\stackrel{}{}\varphi \stackrel{}{}\varphi =\varphi ^2_j(_j\stackrel{}{u}\stackrel{}{u})^2_j_j\stackrel{}{u}_j\stackrel{}{u}.$$

Making use of this inequality, we infer (17) from (18):

$$\stackrel{}{u}_6=\varphi _6\lambda _1\stackrel{}{}\varphi _2\lambda _1\left[𝑑x_{i,j}(_ju_i)^2\right]^{1/2}=\lambda _1\mathrm{𝖼𝗎𝗋𝗅}\stackrel{}{u}_2.$$
(19)
The last equality in (19) is true if $`\mathrm{𝖽𝗂𝗏}\stackrel{}{u}=0`$ and $`\stackrel{}{u}`$ decreases sufficiently fast at spatial infinity, which we require in (17). The coefficient $`\lambda _1`$ is the least possible coefficient in inequality (17) .

Let $`\omega _a`$, $`a=1,2,\mathrm{},r`$, be simple roots of SU(N). They have the same norm $`(\omega _a,\omega _a)\gamma ^2`$. The angle between $`\omega _a`$ and $`\omega _{a\pm 1}`$ is $`2\pi /3`$, and otherwise the roots are perpendicular. Any positive root can be written in the form $`\alpha =\omega _a+\omega _{a+1}+\mathrm{}+\omega _{a+q}`$ for $`a+qr`$. From this it is easy to deduce that all roots have the same norm with respect to the Kartan-Killing form, $`(\alpha ,\alpha )=\gamma ^2`$. To find the actual norm $`\gamma `$, one should compute, say, the matrix $`\widehat{\omega }_1`$ in the Cartan-Weyl basis and take the trace of its square. From (4) it follows that $`\widehat{\omega }_1`$ is block diagonal. The block associated with the Cartan subalgebra is zero because $`\omega _a`$ commute amongst each other. The nontrivial blocks come from the subspaces spanned by $`e_\alpha `$ and $`e_\alpha `$ where the positive root $`\alpha `$ is either equal to $`\omega _1`$ or contains $`\omega _2`$ or $`\omega _1+\omega _2`$. There are $`r1`$ roots containing $`\omega _2`$ and $`r1`$ roots containing $`\omega _1+\omega _2`$. Then $`\gamma ^2=\mathrm{tr}(\widehat{\omega }_1^2)=\gamma ^4N`$ as is deduced from (4). Hence

$$(\omega _a,\omega _a)=N^1.$$
(20)
As a consequence of (20), the following identity holds for any Lie algebra element $`v`$

$$v=\underset{a}{}n_a(n_a,v)+N\underset{a}{}[n_a,[n_a,v]]v_{}+v_{}.$$
(21)
The proof is based on the following observation. Relation (21) is covariant under the adjoint action of SU(N). So, according to (1), $`n_a`$ can be replaced by $`h_a`$ after a corresponding adjoint rotation of $`v`$. Decomposing $`v`$ in the Cartan-Weyl basis, one can see that the first term in (21) is the Cartan component of $`v`$. The double commutator in the second term can be computed by means of (4) and gives rise to the factor $`_a(\alpha ,h_a)^2=(\alpha ,\alpha )=N^1`$ for every basis element $`e_\alpha `$. Thus the second term in (21) is nothing but a projector onto the subspace orthogonal to the Cartan subalgebra spanned by $`n_a`$.

By differentiating (1) one finds

$`_\mu n_a+i[A_\mu ,n_a]=0,`$ (22)
$`i_\mu U^{}UA_\mu {\displaystyle \underset{a}{}}n_aC_\mu ^a,(n_a,A_\mu )=0.`$ (23)
Equation (22) can be interpreted as: The fields $`n_a`$ are transported parallel with respect to the connection $`A_\mu `$. Taking a commutator of (22) with $`n_a`$, summing over $`a`$ and making use of the identity (21), the connection $`A_\mu `$ can be explicitly written via $`n_a`$,

$$A_\mu =iN\underset{a}{}[_\mu n_a,n_a].$$
(24)
The connection (24) has been introduced by Cho to study monopoles in the Yang-Mills theory for SU(2) and SU(3) . By multiplying (22) by $`n_b`$ using the Cartan-Killing form, one deduces that the derivatives of $`n_a`$ are orthogonal to the fields themselves

$$(_\mu n_a,n_b)=0.$$
(25)
Now we show that the tensor (2) is an Abelian gauge field tensor (cf. ), that is, the two-forms $`F^a`$ are closed, $`dF^a=0`$. Consider the following algebraic transformations

$`_\mu A_\nu _\nu A_\mu `$ $`=`$ $`2iN_a[_\nu n_a,_\mu n_a]=2iN_a[[A_\mu ,n_a],[A_\nu ,n_a]]`$ (26)
$`=`$ $`2iN_a\left\{[A_\nu ,[n_a,[n_a,A_\mu ]]]+[n_a,[[n_a,A_\mu ],A_\nu ]]\right\}`$ (27)
$`=`$ $`2i[A_\nu ,A_\mu ]+iN_a[n_a,[n_a,[A_\mu ,A_\nu ]]]`$ (28)
$`=`$ $`i[A_\mu ,A_\nu ]i_an_a(n_a,[A_\mu ,A_\nu ]).`$ (29)
In (26) we have used (24); next, $`_\mu n_a`$ has been transformed via (22); Eq. (27) follows from the Jacobi identity; to derive (28), the first term in (27) has been transformed by means of the algebraic identity (21), while the second one via the Jacobi identity; finally, by applying the algebraic identity (21) to (28), Eq. (29) has been deduced. Introducing the Yang-Mills field strength tensor

$$F_{\mu \nu }=_\mu A_\nu _\nu A_\mu +i[A_\mu ,A_\nu ],$$
(30)
it follows from (29) that

$$F_{\mu \nu }=i_an_a(n_a,[A_\mu ,A_\nu ])=_an_aF_{\mu \nu }^a.$$
(31)
The last equality in (31) is deduced by multiplying (29) and the middle of (26) by $`n_a`$ using the Cartan-Killing form. Now observe that the field strength (curvature) of the pure gauge connection (23) is zero. Making use of the decomposition (23) of a pure gauge connection we obtain

$$0=F_{\mu \nu }_an_a\left(_\mu C_\nu ^a_\nu C_\mu ^a\right),$$

where the identity (22) has been used again for algebraic transformations. Thus, $`F_{\mu \nu }^a=_\mu C_\nu ^a_\nu C_\mu ^a`$. Note that (23) allows one to determine $`C_\mu ^a`$ via the group element $`U`$ explicitly. In (3) the vector potential $`C_i^a`$ can always be chosen to satisfy the Coulomb gauge, $`_iC_i^a=0`$, thanks to the gauge freedom $`C_\mu ^aC_\mu ^a+_\mu \xi ^a`$.

3. A key algebraic inequality. In this section the following inequality is proved

$$_4\kappa _N_2^2,\kappa _N=\frac{gN}{4m^4}\left(1\frac{1}{N^2N}\right).$$
(32)
It is used in the next section to establish the low energy bound. Consider the $`(N^21)\times (N^21)`$ matrix

$$G=\underset{i,b}{}_in_b_in_b.$$

It can be regarded as a linear operator on the Lie algebra, i.e., $`G\eta =_{i,b}_in_b(_in_b,\eta )`$ for any Lie algebra element $`\eta `$. It has $`N1`$ zero eigenvalues because $`Gn_a=0`$. The matrix $`G`$ satisfies

$$\mathrm{tr}G^2\frac{1}{N^2N}\left(\mathrm{tr}G\right)^2.$$
(33)
The proof is simple. Let $`g_k`$, $`k=1,2,\mathrm{},n=N^2N`$, be nonzero eigenvalues of $`G`$. They are real since $`G=G^T`$ with respect to the Cartan-Killing form. Consider a function of one real variable $`\xi `$, $`f(\xi )=_k(g_k\xi )^2`$. Computing the sum explicitly, one finds that $`f(\xi )=\mathrm{tr}G^22\xi \mathrm{tr}G+\xi ^2n`$. The function attains its absolute minimum at $`\xi =\xi _0=\mathrm{tr}G/n`$. Since $`f(\xi )0`$ for all $`\xi `$’s, the inequality (33) follows from $`f(\xi _0)0`$.

Consider a local Cartan-Weyl basis which is obtained by an adjoint transformation of the basis (4), (5) with the group element $`U(x)`$ defined in (1). Denoting $`n_\alpha =U^{}e_\alpha U`$ and $`n_\alpha =U^{}e_\alpha U=\overline{n}_\alpha `$ we have

$$[n_a,n_\beta ]=(h_a,\beta )n_\beta ,[n_\alpha ,n_\alpha ]=U^{}\alpha U\alpha _U,$$
(34)
where $`(\alpha _U,n_a)=(\alpha ,h_a)`$ and $`(n_a,n_\beta )=0`$, $`(n_\alpha ,n_\alpha )=1`$, $`(n_\alpha ,n_\alpha )=(n_\alpha ,n_\alpha )=0`$.

To establish a relation between $`_4`$, $`_2`$ and $`G`$, we decompose the connection (24) in the local Cartan-Weyl basis

$$A_i=\underset{\alpha >0}{}(A_i^\alpha n_\alpha +c.c.)H_i+\overline{H}_i.$$
(35)
Then we obtain

$$_2=m^2\mathrm{tr}G=m^2\underset{i,b}{}([A_i,n_b],[A_i,n_b])=\frac{m^2}{N}\underset{i}{}(A_i,A_i)=\frac{2m^2}{N}\underset{i}{}(\overline{H}_i,H_i).$$
(36)
The second equality follows from (22); the third one is a consequence of (21) and (23). Making use of (31) we also get

$`_4`$ $`=`$ $`{\displaystyle \frac{g}{4}}{\displaystyle \underset{a,i,j}{}}(n_a,[A_i,A_j])^2`$
$`=`$ $`{\displaystyle \frac{g}{4}}{\displaystyle \underset{\alpha ,\beta >0}{}}{\displaystyle \underset{i,j}{}}(A_i^\alpha \overline{A}_j^\alpha c.c.)(\alpha ,\beta )(A_i^\beta \overline{A}_j^\beta c.c.)`$
$``$ $`{\displaystyle \frac{g}{4N}}{\displaystyle \underset{i,j}{}}[(H_i,\overline{H}_j)c.c.]^2.`$ (38)
Note that the local Cartan component of $`[A_i,A_j]`$ can only come from the second commutation relation in (34). Hence the sum over $`a`$ in (S0.Ex5) yields the factor $`_a(\alpha ,h_a)(h_a,\beta )=(\alpha ,\beta )=N^1\mathrm{cos}\theta _{\alpha \beta }N^1`$ for any two positive roots $`\alpha `$ and $`\beta `$ with the angle $`\theta _{\alpha \beta }`$ between them. In a similar fashion we derive

$`\mathrm{tr}G^2`$ $`=`$ $`{\displaystyle \underset{a,b}{}}{\displaystyle \underset{i,j}{}}([A_i,n_a],[A_j,n_b])^2={\displaystyle \underset{a,b}{}}{\displaystyle \underset{i,j}{}}(n_a,[[n_b,A_i],A_j])^2`$
$`=`$ $`{\displaystyle \underset{\alpha ,\beta >0}{}}{\displaystyle \underset{i,j}{}}(A_i^\alpha \overline{A}_j^\alpha +c.c.)(\alpha ,\beta )^2(A_i^\beta \overline{A}_j^\beta +c.c.)`$
$``$ $`{\displaystyle \frac{1}{N^2}}{\displaystyle \underset{i,j}{}}[(H_i,\overline{H}_j)+c.c.]^2.`$ (40)
Combining (40) and (38), we infer

$$4Ng^1_4+N^2\mathrm{tr}G^24\underset{i,j}{}(H_i,\overline{H}_j)(H_j,\overline{H}_i)4\left(\underset{i}{}(H_i,\overline{H}_i)\right)^2=4\left(\mathrm{tr}G\right)^2,$$
(41)
where the Schwartz inequality

$$\left|(H_i,\overline{H}_j)\right|^2(H_i,\overline{H}_i)(H_j,\overline{H}_j)$$

has been used. The inequality (32) immediately follows from (41), (36) and (33).

4. The low energy bound. Let $`\lambda _0=(16\pi ^2N)^1`$. Then Eq. (3) can be written as $`Q=\lambda _0_a\stackrel{}{C}^a,\stackrel{}{B}^a`$. Making use of (15) one gets (cf. the case $`N=2`$ in )

$`|Q|`$ $``$ $`\lambda _0{\displaystyle \underset{a}{}}\stackrel{}{C}^a_p\stackrel{}{B}^a_p^{}`$ (42)
$``$ $`\lambda _0\lambda _1{\displaystyle \underset{a}{}}\mathrm{𝖼𝗎𝗋𝗅}\stackrel{}{C}^a_2\stackrel{}{B}^a_{6/5}`$
$`=`$ $`\lambda _0\lambda _1{\displaystyle \underset{a}{}}\stackrel{}{B}^a_2\stackrel{}{B}^a_{6/5}`$
$``$ $`\lambda _0\lambda _1{\displaystyle \underset{a}{}}\stackrel{}{B}^a_2\stackrel{}{B}^a_2^{1/3}\stackrel{}{B}^a_1^{2/3}.`$ (44)
To get (S0.Ex8), Eq. (17) has been used, which dictated the choice $`p=6`$ in (42), and also $`_iC_i^a=0`$; then the substitution $`\stackrel{}{B}^a=\mathrm{𝖼𝗎𝗋𝗅}\stackrel{}{C}^a`$ has been made; the last inequality (44) is a consequence of (16). The energy can be written as

$$E=m^2\underset{a}{}\stackrel{}{}n_a_2^2+\frac{g}{2}\underset{a}{}\stackrel{}{B}^a_2^2=\underset{a}{}(E_{2a}+E_{4a})=E_2+E_4.$$
(45)
Hence, continuing (44) we get

$`|Q|`$ $``$ $`\lambda _0\lambda _1\left(2/g\right)^{2/3}{\displaystyle \underset{a}{}}\left(E_{4a}^{4/3}\right)^{1/2}\left(\stackrel{}{B}^a_1^{4/3}\right)^{1/2}`$ (46)
$``$ $`\lambda _0\lambda _1\left(2/g\right)^{2/3}\left[\left({\displaystyle \underset{a}{}}E_{4a}\right)^{4/3}\right]^{1/2}\left[\left({\displaystyle \underset{b}{}}\stackrel{}{B}^b_1\right)^{4/3}\right]^{1/2}`$
$`=`$ $`\lambda _0\lambda _1\left(2/g\right)^{2/3}E_4^{2/3}\left[{\displaystyle \underset{b}{}}{\displaystyle 𝑑x\sqrt{\stackrel{}{B}_b\stackrel{}{B}_b}}\right]^{2/3}`$
$``$ $`\lambda _0\lambda _1g^1\left(2E_4\right)^{2/3}(N1)^{1/3}\left[{\displaystyle 𝑑x\sqrt{2_4}}\right]^{2/3}`$ (48)
$``$ $`\lambda _0\lambda _1g^1\left(2E_4E_2\right)^{2/3}\left[2(N1)\kappa _N\right]^{1/3}`$ (49)
$``$ $`\lambda _0\lambda _1g^12^{2/3}\left[2(N1)\kappa _N\right]^{1/3}E^{4/3}`$
$`=`$ $`c_N^{4/3}E^{4/3},`$
where the constant $`c_N`$ is given in (8). To get (S0.Ex9), the Hölder inequality (13) for $`p=q=1/2`$ and then (11) for $`\gamma =4/3`$ have been applied; (48) follows from (12); the algebraic inequality (32) has been used to deduce (49); the final result (S0.Ex10), which is equivalent to (7), comes from the basic inequality $`E_4E_2E^2/4`$.

If the Lagrangian (6) defines an effective theory of the SU(N) gauge fields in some approximation, then the coefficients $`m`$ and $`g`$ should depend on $`N`$, the Yang-Mills coupling constant and a mass scale that determines the energy range in which the approximation is valid.

The result (7) is trivially generalized to the case when the mass scales and coupling constants are different for each mode $`n_a`$, that is, $`m^2`$ and $`g`$ are replaced by $`m_a^2`$ and $`g_a`$, respectively, and inserted into the corresponding sums over $`a`$ in (6) (cf. ). In this case, $`m^2=\mathrm{max}_a\{m_a^2\}`$ and $`g=\mathrm{max}_a\{g_a\}`$ in (8). Indeed, all the inequalities in sections 3 and 4 still hold as a consequence of $`m_a^2m^2`$ and $`g_ag`$. Equality (36) becomes an inequality $`m^2\mathrm{tr}G_2`$, which, however, does not affect our derivation of (32) from (41).

Acknowledgments. I wish to thank John Klauder and David Metzler for useful discussions, David Groisser for the reference , and Richard Ward for helpful comments and the reference .


text:
# Force and kinetic barriers to unzipping of the DNA double helix

## Abstract

A theory of the unzipping of double-stranded (ds) DNA is presented, and is compared to recent micromanipulation experiments. It is shown that the interactions which stabilize the double helix and the elastic rigidity of single strands (ss) simply determine the sequence dependent $`12`$ pN force threshold for DNA strand separation. Using a semi-microscopic model of the binding between nucleotide strands, we show that the greater rigidity of the strands when formed into dsDNA, relative to that of isolated strands, gives rise to a potential barrier to unzipping. The effects of this barrier are derived analytically. The force to keep the extremities of the molecule at a fixed distance, the kinetic rates for strand unpairing at fixed applied force, and the rupture force as a function of loading rate are calculated. The dependence of the kinetics and of the rupture force on molecule length is also analyzed.

Introduction: In cells, proteins apply forces to unzip and stretch DNA. These forces can be studied in single-molecule experiments (Fig. 1), and are of biophysical as well as biological interest. Our focus here is primarily on unzipping experiments where forces are applied across the double helix to adjacent 5’ and 3’ strands (Fig. 1A). In experiments, the control parameters may be the force $`f`$ itself, the distance between the last base pairs $`2r`$, or the rate of force increase or ‘loading rate’ (Fig. 1B). We discuss the results expected in all these situations.

We first use a thermodynamical equilibrium approach to show that the sequence-dependent force associated with unzipping of large DNAs, $`f_u12`$ piconewton (pN), can be simply deduced from the known free energy of DNA denaturation and the elasticity of single-stranded DNA. The unzipping experiments of Essevaz-Roulet et al (Fig. 1A)  and Rief et al (Fig. 1D)  are accurately described at this macroscopic level.

Other experimentally observable aspects of unzipping can only be investigated using a more detailed description of base-pairing interactions. We therefore present a semi-microscopic model which accounts for hydrogen bonds and stacking interactions . We show that a free-energy potential barrier originates from the greater range of conformational fluctuation of DNA strands when isolated, relative to when they are bound together to form dsDNA (Fig. 2). Our model can be investigated in detail and allows precise calculation of the effects of this barrier for the initiation of unzipping and the kinetics of strand dissociation.

We compute the force necessary to keep apart the two extremities of the DNA molecules at some distance $`2r`$, as well as the shape of the opening fork (Fig. 3). Due to the potential barrier, this force is much larger at small $`r`$ (and can reach some hundreds of pN) than the asymptotic value $`f_u`$ at large $`r`$. Analysis of unzipping in thermal equilibrium at the high level of precision possible in AFM experiments would allow unambiguous verification of this predicted force barrier.

The barrier makes strand dissociation an activated process with dynamics that can be analyzed using nucleation theory. Unzipping starts with a transition ‘bubble’ a few ($`4`$) bases long (Fig. 4). We calculate the free energy of this bubble, and determine how the dissociation rate depends on applied force and molecule length (Fig. 5). Results are compared to the experiments of Bonnet et al.  and of Pörschke .

Extending Evans’ theory for the breaking of single bonds  to the case of a one-dimensional polymer , we then calculate the most probable rupture force when the DNA molecule is subjected to a force which increases at a constant ‘loading rate’ (Fig. 6). The dependence of the rupture force upon loading rate and molecular length could be quantitatively tested by AFM unzipping experiments; these results also shed light on AFM DNA-stretching experiments of Struntz et al. and of Rief et al..

Thermodynamic Description of Unzipping: An unstressed double helix is stabilized against spontaneously dissociating into its two strands by the interaction free energy per base pair, which from a thermodynamic perspective we may take to be some average amount $`g_0`$. Although dependent on sequence, we may consider $`g_0=1.4k_BT`$, the value determined from single molecule experiments on an AT rich sequence in $`\lambda `$ phage , as a reference for the free energy difference between dsDNA and separated ssDNAs. Our emphasis is on an understanding of the free-energy balance in unzipping rather than to study inhomogenous sequence effects.

In the presence of applied torque $`\mathrm{\Gamma }`$ and unzipping force $`f`$ (Fig. 1B) the free energy difference per base pair between unzipped and base-paired DNA strands is

$$\mathrm{\Delta }g(\mathrm{\Gamma },f)=g_{\mathrm{ssDNA}}g_{\mathrm{dsDNA}}=g_0+\theta _0\mathrm{\Gamma }+2g_\mathrm{s}(f)$$
(1)
When $`\mathrm{\Delta }g<0`$, opening is thermodynamically favorable. The last two terms represent mechanical work done per base pair unzipping the double helix. In order, they are the work done by the torque ($`\theta _0=2\pi /10.5`$ is simply the change in strand winding angle during conversion of dsDNA to separated strands), and the stretching free energy of the unzipped single strands.

The function $`g_\mathrm{s}(f)`$ in (1) is the ssDNA stretching free energy per base at fixed force. The leading factor of 2 simply reflects the fact that two bases of ssDNA are created for each base-pair of dsDNA which is unzipped. The ssDNA elastic behavior is complicated by nucleotide-interaction effects , but experimental force-extension curves for $`\lambda `$ phage ssDNA in 150 mM Na+ are well described by a freely-jointed chain-like (FJCL) elastic response for forces $`>1`$ pN, with Kuhn length $`d=15`$Å . The corresponding free energy for forces up to $`20`$ pN is

$$g_\mathrm{s}^{FJCL}(f)=k_BT\frac{l_{ss}}{d}\mathrm{log}\left[\frac{\mathrm{sinh}(df/[k_BT])}{df/[k_BT]}\right]$$
(2)
where the contour length per base pair is $`l_{ss}=27`$ $`\mu `$m$`/48.5`$ kb $`5.6`$ Å.

At zero applied force, $`g_\mathrm{s}(0)=0`$, thus ssDNA is stable when $`\mathrm{\Gamma }<\mathrm{\Gamma }_u=g_0/\theta _0=2.4k_BT`$, in good agreement with an experimental estimate of the twisting torque needed to denature an AT rich sequence in a $`\lambda `$ phage DNA  (the sign indicates a left-handed dsDNA-unwinding torque). In this case, the work done by the torque during opening is simply $`g_0=1.4k_BT`$ per base pair.

In the opposite case where torque is zero (i.e. for dsDNA with no constraint of its strand linking number), the critical unzipping force (at which $`\mathrm{\Delta }g(f)=0`$) is $`f_u^{FJCL}=11`$ pN. These results are in good agreement with the mid-range of unzipping forces encountered with experiments on inhomogeneous-sequence DNAs by Essevaz-Roulet et al. (Fig. 1A: 12 pN threshold to start, then 10 to 15 pN during unzipping $`\lambda `$-DNA). Data of Rief et al. (Fig. 1D; $`20\pm 3`$ pN for poly(dG-dC), $`9\pm 3`$ pN for poly(dA-dT)) gives via (1) denaturation free energies of 1.1 and 3.5 $`k_BT`$ per AT and GC base pair respectively, in good agreement with thermodynamical data. Finally, the projected length of one ssDNA nucleotide along the force direction at the unzipping transition is given as $`d_u^{FJCL}4`$Å. Bockelman et al. used a similar theory to analyze unzipping force dependence on sequence.

The inset of Fig. 3 shows the curve in the torque-force plane on which $`\mathrm{\Delta }g=0`$, which is the ‘phase boundary’ separating dsDNA and unzipped ssDNAs. This boundary is predicted to have the shape $`f_u(\mathrm{\Gamma }\mathrm{\Gamma }_u)^{1/2}`$ for small $`f_u`$.

For forces up to 15 pN, (2) is approximated to $`0.15k_BT`$ accuracy by the simple quadratic form

$$g_\mathrm{s}^G(f)=\frac{f^2}{C}.$$
(3)
where the ssDNA effective elastic constant is $`C=0.12k_BT/`$Å2. Using this form allows analytical solution for the unzipping force, $`f_u^G=(C|g_0|/2)^{1/2}=`$ 12 pN ; at this force the projection of ssDNA monomer length along the force direction is $`d_u^G=(2|g_0|/C)^{1/2}5`$Å. This quadratic approximation is quantitatively nearly equivalent to the nonlinear model (2) at forces up to $`15`$ pN (e.g. note the accord between the torque-force ‘phase boundaries’ in Fig. 3, inset); this will be key to the continuum theory below.

Semi-Microscopic Model of Strand Binding: Features at the nucleotide scale relevant to the onset of unzipping are ignored in models like (1). We therefore move to a model which uses the distances $`2r(n)`$ between corresponding $`n`$th base pairs of the two strands, as degrees of freedom. The energy of the DNA strands is:

$$H=_0^N𝑑n\left\{\frac{1}{2}m\left(r(n)\right)\left(\frac{dr}{dn}\right)^2+U\left(r(n)\right)\right\}$$
(4)
This model is similar to models previously used to describe thermal denaturation. The precise form of (4) follows from our previously developed model for denaturation by torque  by integration over angular degrees of freedom, followed by continuum limit for the base index $`n`$.

The first term in (4) describes interactions between neighboring bases along each strand, and must depend strongly on the inter-base half-distance $`r(n)`$, since conformational fluctuations of the strands are highly quenched inside the double helix, relative to those along ssDNA. The strand rigidity is $`m(r)=Ee^{b(rR_0)}+C`$, where $`R_0=10`$ Å is the double helix radius, and where $`1/b=0.6`$ Å is the separation at which the strand rigidity changes from its double-helix value $`E+C`$, to the much smaller ssDNA value $`C`$. We use the value $`C=0.12`$ $`k_BT/`$Å2 from the previous section. The rigidity of the strands inside the double helix has been determined from Raman measurements of internal vibrations of dsDNA, to be $`E=58k_BT/`$Å2.

The second term in (4) is a potential, acting between the two strands, made up of the hydrogen-bonding energy between corresponding bases, plus a torque energy: $`U(r)=U_H(r)\mathrm{\Gamma }R_1/r`$ where $`R_1=6`$ Å. We use the Morse potential form for the hydrogen bonding interaction: $`U_H(r)=D[(e^{a(rR)}1)^21]`$, with $`D=5.84k_BT`$ and $`a=6.3`$ Å-1 (Fig. 2 inset).

In thermal equilibrium, strand unpairing is described by the partition function $`Z=𝒟r(n)e^{H/k_BT}`$. $`Z`$ can be computed by use of a continuum transfer matrix technique along the $`n`$ coordinate, leading to a Schrödinger-like equation:

$$\left[\frac{(k_BT)^2}{2m(r)}\frac{^2}{r^2}+V(r)\right]\psi (r)=g\psi (r).$$
(5)
The free-energy potential $`V(r)=U(r)+(k_BT/2)\mathrm{ln}[m(r)/m(\mathrm{})]`$ includes an entropic contribution due to the decrease in the rigidity $`m(r)`$ with strand unbinding. This entropic potential (Fig. 2) arises when going from the path integral to the Schrödinger equation  with non constant mass, and here it generates a large force barrier to the initiation of unzipping, which strongly affects the kinetics of strand separation.

The lowest eigenvalue $`g_0`$ is the equilibrium free energy, and the corresponding $`\psi _0(r)`$ describes thermal fluctuations of distance between the two strands. For $`g0`$ the eigenvalue spectrum becomes continuous, corresponding to the appearance of completely separated ssDNAs at $`g=0`$. When DNA is the thermodynamical favorable state $`g_0<0`$ represents the free energy per base pair of a long dsDNA, relative to separated ssDNAs. The fluctuations of $`r`$ are confined to the Morse well, so $`g_0`$ is well approximated if we take $`V(r)=V(R_0)+U_H(r)U_H(R_0)`$ and $`m(r)=m(R_0)`$ for which (5) is exactly soluble. Inside the well, $`\psi _0(r)`$ is the Morse ground state ; outside the well it can be computed using the WKB approximation. Using the parameters listed above, and at zero torque, the lowest eigenvalue is $`g_0=1.4k_BT`$, in accord with the corresponding number assumed in the previous section. Application of an unwinding torque gradually increases $`g_0`$ until at $`\mathrm{\Gamma }_u=2.4k_BT`$ it becomes zero and the DNA unwinds, exactly as occurs in the thermodynamic model of the previous section (Fig. 3, inset).

Force required to hold ssDNA ends at a given distance: This situation can be analyzed simply in terms of $`\psi _0(r)`$, without further computation. This is because, for a semi-infinite dsDNA, the function $`\psi _0(r)`$ corresponding to $`g_0`$ is the probability distribution for the two ssDNA ends to fluctuate a distance $`2r`$ apart. Therefore the free energy associated with a fluctuation which separates the two ssDNA ends by a distance $`2r`$ (or equivalently the total work done separating the two ends to a distance $`2r`$) is

$$W(r)=k_BT\mathrm{log}\psi _0(r)$$
(6)
up to an additive constant which is unimportant for our analysis. Note that, neglecting surface interactions, this is half the free energy associated to separate a distance $`2r`$ in the middle of a long dsDNA, since an interior ‘bubble’ is made of two ‘forks’.

The unzipping force that must be supplied to hold the two ssDNA ends a distance $`2r`$ (Fig. 1B) apart is thus just the derivative of (6), $`f(r)=dW(r)/d(2r)`$. This is in Fig. 3, which displays a large force barrier of $`270`$ pN as the strands are forced apart. The barrier peak occurs for a half-separation $`rR_0=0.5`$ Å, and then decays to the long-molecule unzipping force $`f_u12`$ pN by $`rR_04`$ Å. At large distances, the work done per base pair by unzipping is $`2f_ud_u3k_BT`$, twice as much as the denaturation free energy $`g_0`$; this is because the force must unzip the DNA, and extend the highly flexible ssDNAs. The peak force is large compared to the fluctuations in force associated with sequence. The force barrier will not be observable in large-scale unzipping experiments, but should be observable in AFM studies. A stiff cantilever with roughly 0.1 Å thermal noise should be used to measure the force barrier as a function of essentially fixed opening distance.

In this fixed-distance experiment, one might also measure the shape of the opening ‘fork’, by determining the relation between opening distance and base position, $`n(r)`$ (Fig. 1B). The most probable configuration $`n(r)`$ satisfies the equation of motion associated with (5), which expresses force balance along the chain,

$$m(r)\frac{d^2r}{dn^2}+\frac{1}{2}m^{}(r)\left(\frac{dr}{dn}\right)^2=V^{}(r)$$
(7)
We integrate (7) to obtain $`r(n)`$, the shape of the opening ‘fork’. Starting from the opening point where $`r=R_0`$ and $`dr/dn=0`$, $`n(r)=2[E(rR_0)/V^{}(R_0)]^{1/2}`$. Far from the opening point $`n(r)=r/d_u^Gr/(5`$Å).

Unzipping kinetics at fixed force: Many experiments on short (10-100 bp) dsDNAs ( and caption of Fig. 1) probe the kinetics of strand separation. The equilibrium results discussed above are a starting point for a kinetic theory of unzipping based on nucleation theory. The Schrödinger equation (5) with the fluctuation-corrected potential $`V(r)`$ describes dsDNA and ssDNA as locally stable molecular states. The general problem faced in unzipping kinetics is the transition from an initially metastable state (dsDNA or ssDNA, depending on the force, see below) to a final, stable (lower-free-energy) state. Strand dissociation requires the whole polymer chain to cross the free energy potential barrier of $`V(r)`$ (Fig. 2), which makes the transition rate strongly length and force dependent. Using the effective potential $`V(r)`$ (including the entropic barrier) corresponds to averaging over microscopic fluctuations of individual bases, restricting us to consider unzipping rates slow compared to those of these microscopic fluctuations (the experimentally relevant regime).

The transition rate (equivalently the inverse lifetime) has the form familiar from transition-state theory :

$$\nu =\nu _0e^{G^{}/k_BT}$$
(8)
It requires the activation free energy $`G^{}`$ of a transition state, relative to the initial metastable state. The transition state is the saddle-point configuration of the free energy, with one unstable direction leading monotonically down to the initial and final states, and is the dominant transition pathway.

Our transition states are just the partially unzipped configurations $`r^{}(n)`$ determined from (7), for force-dependent boundary conditions consistent with the initial ($`r_i`$) and final ($`r_f`$) states. The activation free energy is

$`G^{}`$ $`=`$ $`H[r^{}(n)]+(k_BT/2){\displaystyle 𝑑n\mathrm{ln}[m[r^{}(n)]/m(\mathrm{})]}2fr^{}(0)G_m`$ (9)
$`=`$ $`{\displaystyle _{r_i}^{r_f}}𝑑r\sqrt{2m(r)(V(r)g_m)}2f(r_fr_i),`$ (10)
i.e. the free energy associated with (7) minus the free energy of the metastable state $`(G_m=Ng_m)`$ from which the transitions occur. $`G^{}`$ is in practice the free energy of the few-base-pair ‘bubble’ portion of the transition state which separates the unzipped and double-stranded regions.

Finally, the rate prefactor $`\nu _0`$ is the linear growth rate of unstable perturbations around the saddle-point configuration. We assume viscous dynamics with a friction coefficient per base $`\zeta =6\pi \eta R_0`$ (water viscosity $`\eta =1\times 10^3\mathrm{kg}/(\mathrm{m}\mathrm{sec})`$). From a detailed calculation we find $`\nu _0Da^2/(4\zeta )=1\times 10^{12}`$ sec-1, essentially the ratio of the negative curvature of the Morse potential near the top of the well, to the friction coefficient. We now describe how to compute the dissociation rate and the nucleation bubble shape (i.e. boundary condition $`r_i,r_f`$) depending on the unzipping force.

Kinetics of unzipping where ssDNA is stable ($`f>f_u`$): If a steady force is applied which is slightly bigger than the equilibrium unzipping threshold $`f_u`$, then the initial dsDNA is metastable relative to separated strands: $`\mathrm{\Delta }g(f)=g_02f^2/C<0`$ (1), and formula (8) is directly applicable (Fig. 4A) to calculation of the dissociation rate $`\nu _{}`$. The free energy of the initial metastable state is $`g_m=g_0`$. The transition state in this case is a short ssDNA ‘nucleus’ of $`n^{}`$ bases at the open end ($`n=0`$) of the dsDNA, with boundary conditions at the ends of fixed force $`f`$ ($`dr/dn=2f/m[r(n)]`$), and free energy equal to the metastable dsDNA value ($`V(r(n))2f^2/m[r(n)]=g_0`$). This determines $`r_i=r(n^{})`$ and $`r_f=r(0)`$.

The size of the nucleation ‘bubble’ depends weakly on the force and is $`n^{}(f)4`$ bases, independent of the overall DNA molecule length. Therefore, the dissociation time $`t_{}=1/\nu _{}`$ (Fig. 5) is length-independent for forces above $`f_u=12`$ pN. The decrease in dissociation time with forces $`>f_u`$ is due to reduction of $`G^{}`$ by the applied force. Beyond $`f_b230`$ pN, the barrier in $`V(r)`$ is completely overcome, and unzipping is immediate ($`t_{}=1/\nu _0`$). Our computation addresses only the initial unzipping barrier-crossing event and does not include the time necessary to push the fork down the dsDNA, which would introduce a weak molecular-length dependence for forces $`>f_u`$ in Fig. 5. Fork motion near $`f_u`$ may also be retarded by additional barriers introduced by inhomogeneous sequence. Kinetics of unzipping where ssDNA is metastable ($`f<f_u`$): Below the equilibrium unzipping force threshold, an infinitely long molecule is stable as a dsDNA. However, relatively short ($`<20`$ bp) oligomers have a finite strand-dissociation time which can be accelerated by applied force. Here the final ssDNA state is metastable (although with a possibly long lifetime) relative to the initial dsDNA state. Therefore our calculation scheme applies not to unzipping (off-rate $`\nu _{}`$), but instead to annealing of ssDNAs (on-rate $`\nu _+`$, Fig. 4B). We compute the off-rate from the on-rate, using the equilibrium condition $`\nu _{}/\nu _+=e^{N\mathrm{\Delta }g(f)/k_BT}`$, where $`N`$ is the number of base-pairs in the molecule, and where $`\mathrm{\Delta }g(f)`$ is the force-dependent free energy difference per base pair between the paired and unpaired states .

The calculation of the on-rate requires computation of a saddle point configuration with a nucleation bubble of dsDNA at the end of ssDNAs to which the force is applied (Fig. 4B). The relevant boundary conditions are analogous to those above, apart from the requirement that the ssDNAs be in close proximity but not base-paired (i.e. slightly outside the potential barrier of Fig. 2) and the free energy of the metastable initial state equals $`g_m=2g_s(f)=2f^2/C`$. The dsDNA nucleation ‘bubble’ remains near 4 bp as the force $`f`$ decreases from $`f_u`$ to zero, leading to an essentially $`N`$-independent on-rate $`\nu _+`$ (again the time necessary for the actual ‘zipping up’ is not included). Then, the $`N`$-dependence from the energy difference of paired and unpaired DNA results in a strong molecular-length dependence of $`\nu _{}`$ as shown in Fig. 5 for forces less than $`f_u`$.

At zero force, this becomes a calculation of dissociation time for free dsDNAs in solution with result $`t_{}=10^{0.6N6.3}`$ sec. The inset of Fig. 5 shows this as a function of dsDNA length; a 10 bp dsDNA has a lifetime of roughly 1 sec; the exponential length-dependence results in a 30 bp DNA being stable for $`10^{12}`$ sec$`30000`$ yr. Our estimate of $`t_{}`$ for $`N=5`$, $`t_{}0.3`$ msec, is in agreement with the results of Bonnet et al.  (Fig. 1E). The prediction for the nucleation bubble size, $`n^{}=4`$, is close to the value $`n^{}=3`$ measured by Pörschke for a poly(A)poly(U) acid.

Rupture of dsDNA during gradual loading: The previous two sections discuss fixed-force experiments; an alternate experiment is to steadily increase force ($`f=\lambda t`$ where $`\lambda `$ is the ‘loading rate’ in pN/sec) and then to measure the force at which rupture occurs. Generally, rapid loading rates result in rupture at large forces. We are able to predict the most probable rupture force versus loading rate and molecular length, and our results display a rich range of possibilities (Fig. 6).

Using the calculations of off-rates presented above, the probability distribution for rupture with force is

$$P(f)=\nu _{}(N,f)\mathrm{exp}\left(\frac{1}{\lambda }_0^f𝑑f^{}\nu _{}(N,f^{})\right)$$
(11)
For a number of molecular lengths, Fig. 6 shows the location of the peak of this distribution, the most probable rupture force that would be measured experimentally. For sufficiently slow loading rate, rupture occurs at zero force simply by thermal dissociation; for molecules $`>20`$ bp thermal dissociation is practically unobservable.

At some loading rate $`\lambda _1`$, the peak in $`P(f)`$ jumps to finite force and then increases, with each length of molecule following a different curve. In this regime, the forces at rupture are typically below the equilibrium unzipping threshold $`f_u`$, since there is time for many thermal attempts at barrier-crossing to the metastable ssDNA state during loading. The length-dependence follows from the calculation of formation of metastable ssDNA discussed above, and has been qualitatively observed for 5’-5’ pulling experiments (Fig. 1C).

At a higher loading rate $`\lambda _2`$ the peak of the rupture force distribution hits $`f_u`$ which remains the most likely rupture force up to a loading rate $`\lambda _3`$. For $`\lambda >\lambda _2`$, dissociation is occuring at forces large enough that the ssDNA final state is stable, resulting in no $`N`$-dependence. Finally, beyond the very large loading rate $`\lambda _310^{5.5}`$ pN/sec, the rupture force gradually increases simply because the molecule is unable to respond to the force before it becomes very large. $`\lambda _3`$ separates equilibrium and nonequilibrium time scales for very long sequences ($`N\mathrm{}`$): the rupture force is independent on the loading rate and equal to $`f_u`$ when $`\lambda <\lambda _3`$, and increases above.

Our framework treats nonequilibrium rupture of a one dimensional object, a development of previous theory necessary for interpretation of unzipping experiments. In light of our detailed predictions, 3’-5’ AFM unzipping experiments (Fig. 1A) should yield interesting results, and would be of help for an accurate determination of the free-energy potential $`V(r)`$.

We thank C. Bouchiat, G. Bonnet, V. Croquette, F. Pincet for useful discussions. Work at UIC was supported by NSF Grant DMR-9734178, by the Petroleum Research Foundation of the American Chemical Society and by the Research Corporation. S. Cocco is partly funded by A. della Riccia grant. R. Monasson is supported in part by the MRSEC Program of the NSF under Award number DMR-9808595.

$`^{}`$