text:
# Phase diagram of force-induced DNA unzipping in exactly solvable models

## I Introduction

In recent years, micro and nanomanipulation of single biological macromolecules has become feasible due to a dramatic improvement of experimental techniques. By using devices such as optical tweezers , soft microneedles  and atomic force microscopes , it is now possible to access physical and mechanical properties of fundamental biological objects, namely proteins, nucleic acids, and molecular motors, on the scale of individual molecules. Special effort has been devoted to the measurement of force-elongation characteristics of double stranded DNA molecules (dsDNA), determining its response to external forces and torques in the absence of enzymes. The mechanical unzipping of dsDNA structure by a force pulling the end of one of the two strands, the end of the other strand being anchored to some physical support, has been studied by Bockelmann et al., who have measured the average force along the opening of the two strands. Mechanical forces are in fact exerted on the DNA molecule by different enzymes during the process of DNA replication or transcription . On the other hand, the double helical structure of dsDNA may be disrupted ’in vitro’ by changing pH, solvent conditions and/or temperature . This transition is known as melting denaturation, and it has been long studied by theoretical physicists . Instead, only very recently DNA mechanical denaturation has been the subject of theoretical studies. Most of these studies have considered so far a simple extension of the Poland-Sheraga model , in which the two DNA strands are homogeneous ideal polymer chains interacting with each other, introducing a constant force pulling apart the two strands by acting on one extremity of both strands. By using a mapping into a quantum mechanical problem, it has been shown  that the opening of the two strands occurs only if the pulling force is increased to a critical value. The unzipping transition turns out to be a first order phase transition, whereas the melting transition is second order, in the ideal case, in $`d=3`$ (recently both simulations and exact results have shown that the melting transition becomes first order, when considering mutually self-avoiding walks (SAWs) ). The heterogeneous case has also been studied with similar techniques . Very recently, Monte Carlo simulations of interacting pairs of self avoiding walks on a cubic lattice and of bead-and-string chains in the continuum space have determined the whole phase diagram in the force-temperature plane, revealing the existence of a re-entrant zipping-unzipping transition by decreasing the temperature for a finite range of forces .

In this paper, we obtain exact analytical results for a class of simple lattice models of interacting pairs of homogeneous directed self-avoiding walks. We extend the model introduced in ref.  in $`D=1+1`$ to the generic dimension case $`D=d+1`$, and analyze both the scaling behavior of the first order unzipping transition and the multicritical scaling laws at the melting transition in the force-temperature plane. We also consider different versions of the model depending on whether the crossing of the two walks is allowed or not, or simply penalized by an entropy cost. In all cases, we find that the critical pulling force increases with temperature at low $`T`$, implying thus the existence of a cold unzipping transition . This seemingly paradoxical property is due to the competition between the energy gained by increasing the open portion of DNA and the entropy lost with the full stretching of the separated strands. Cold unzipping will be exactly proved for both ideal and self-avoiding chains in the lattice, and for a discrete chain with constant distance between consecutive beads in the continuum space. We also discuss the role of denatured bubbles forming in dsDNA opening, as opposed to the end-opening of the strands induced by the pulling force. By comparing the phase diagram of the different models considered, we point out that the approximation of neglecting bubbles and considering only Y-shaped configurations shares many of the features of a mean-field type approximation, with an upper critical dimension $`d_c=4`$.

Our paper is organized as follows. In Sec.II we introduce the models of directed walks on the lattice in any dimension $`D`$ and compute their thermodynamical properties, deriving the phase diagram in the force-temperature plane. In Sec.III we analyze the scaling laws at both the thermal melting and mechanical unzipping transition, and we highlight the features that are physically more interesting in our exact results (namely the existence of a re-entrant cold unzipping transition, our main result, and the role of forbidding the mutual crossing of the two strands). In Sec.IVa we discuss the physical meaning of re-entrance in our (lattice) models, generalizing the proof of its occurrence to the more realistic case of SAWs. In Sec.IVb we prove that also a discrete chain in the continuum space undergoes cold unzipping, and then argue why such an effect had not been observed in previous calculations , performed in the limit of a continuum chain. Eventually, we discuss some related models (with generic SAWs) and give some additional technical details in the Appendices.

## II Exact solution of the lattice models

### A Introduction of the models and their behaviour for thermal melting

We consider a simple class of models in a $`D=d+1`$ dimensional hypercubic lattice $`𝐙^𝐃`$. The two strands of a homogeneous DNA molecule with $`N`$ base pairs are mimicked by two SAWs, directed along the $`(1,\mathrm{},1)`$ direction. The two chains have one end in common, while at the other end a force $`\stackrel{}{g}`$ is pulling in the $`(1,1,0,\mathrm{},0)`$ direction. The walks gain a binding energy $`ϵ`$ ($`ϵ>0`$) every time the bases with the same monomer index (or same projection upon $`(1,\mathrm{},1)`$) interact, i.e. wrong base pairing is forbidden.

The canonical partition function for two $`N`$-step directed self avoiding walks (DSAWs) is:

$$Z_N(\beta ,\stackrel{}{g})=\underset{\stackrel{}{x}𝐙^𝐃}{}p_N(\stackrel{}{x},\beta ϵ)\mathrm{exp}(\beta \stackrel{}{g}\stackrel{}{x}),$$
(1)
where $`p_N(\stackrel{}{x},\beta ϵ)`$ is the canonical partition function of two directed interacting strands whose last base pairs are at relative distance $`\stackrel{}{x}`$. The following recursion relation holds for $`p_N(\stackrel{}{x},\beta ϵ)`$:

$$p_{N+1}(\stackrel{}{x},\beta ϵ)=\underset{i,j=1}{\overset{D}{}}p_N(\stackrel{}{x}\stackrel{}{e}_i+\stackrel{}{e}_j,\beta ϵ)(1+(\mathrm{exp}(\beta ϵ)1)\delta _{\stackrel{}{x},\stackrel{}{0}}),$$
(2)
where $`\stackrel{}{e}_i`$, $`i=1,\mathrm{},D`$, are the canonical euclidean versors of the $`D`$-dimensional space and $`\delta `$ is the Kronecker delta.

These and similar equations have been intensively studied, at $`\stackrel{}{g}=\stackrel{}{0}`$, within simple models of DNA thermal melting , and, in $`D=2`$, within the context of random walk adsorption  and wetting problems . Note that the choice of the $`(1,\mathrm{},1)`$ direction along which the walks are directed, is crucial in allowing us to write local recursion relations. It is thus no surprise that the model belongs to the same universality class of random walks in $`d=D1`$ dimensions .

It is well known (see  for a recent review), that, in the absence of a pulling force, within the model defined by equation (2), dsDNA undergoes a phase transition between a low temperature double stranded phase and a high temperature denaturated state only for $`D>3`$ ($`d>2`$): in $`D=2,3`$ dsDNA remains double stranded at all temperatures.

However, we note that in (2) the two strands are allowed to cross each other, without any restriction. This seems rather unphysical, since every real chain has a finite ”hard core” distance. Thus, one should consider also the case in which crossing between different strands is forbidden (recent studies of both homogeneous  and heterogeneous  DNA melting transition have indeed considered the $`d=1`$ case with forbidden crossing). Let us focus for the moment on the $`D=2`$ case, with no pulling force, where the effect of forbidding crossing is most dramatic. Here, the direction $`(1/\sqrt{2},1/\sqrt{2})`$ can be identified as ”time”, and its normal $`(1/\sqrt{2},1/\sqrt{2})`$ as ”space” ($`x`$). The model with crossing (w.c.) implies no restriction on the relative distance $`x`$, whereas the one without crossing (w.o.c.) implies that $`x`$ cannot change sign, e.g. $`x0`$. The model w.o.c. is equivalent to surface adsorption models previously considered . While the thermal melting of the two strands w.c. takes place at $`T_c=\mathrm{}`$, in the model w.o.c. $`T_c=ϵ/\mathrm{log}\left(4/3\right)`$ as it can be deduced from ref.. To elucidate this point further, we have tackled an intermediate model, where we do not forbid crossing, but we make it disadvantageous, by assigning a cost $`V>0`$ each time the strands pass through one another. With calculations similar to those reported below, we find that melting takes place at the critical temperature (as a function of $`V`$):

$$T_c(V)=ϵ\left[\mathrm{log}\left(\frac{4e^V}{3e^V+1}\right)\right]^1.$$
(3)
As expected, as $`V\mathrm{}`$, equation (3) yields $`T_cϵ/\mathrm{log}\left(4/3\right)`$, whereas $`T_cϵ/V`$ as $`V0`$.

To implement the non-crossing constraint in the model in generic dimension, we require that, if the DSAWs join, one coming from direction $`\stackrel{}{e}_i`$ and the other from direction $`\stackrel{}{e}_j`$ (with $`ij`$), when they divide again, they cannot both proceed along their previous direction; i.e we forbid that the first walk goes along direction $`\stackrel{}{e}_i`$ and the second one along $`\stackrel{}{e}_j`$. This excludes one configuration out of $`D(D1)`$ at the splitting point (the remaining $`D`$ possibilities would lead to the zipping of the two strands). As regards the $`D=3`$ thermal melting, again the model w.c. has $`T_c=\mathrm{}`$, while, if crossing is forbidden as described above, the critical temperature is $`T_c=ϵ/\mathrm{log}\left(9/8\right)`$. This result was also found in a similar calculation by Rubin . For $`D>3`$ the models with and without crossing both undergo a denaturation transition at a finite, though different, critical temperature, so that the effect of forbidding crossing here is less important (the critical temperatures are the same at the leading order $`1/D`$ in the $`D\mathrm{}`$ expansion, as expected).

### B Behaviour of the models at non zero force

Let us now turn to the models with the pulling force $`\stackrel{}{g}`$ in generic dimension $`D`$, with $`\stackrel{}{g}=(g,g,0,\mathrm{},0)`$. We can find the asymptotic value of the canonical partition function (1) by locating  the singularity closest to the origin of its generating function:

$$Z(z,\beta ϵ,\beta g)=\underset{N=0}{\overset{\mathrm{}}{}}z^N\underset{\stackrel{}{x}}{}p_N(\stackrel{}{x},\beta ϵ)\mathrm{exp}(\beta \stackrel{}{g}\stackrel{}{x}),$$
(4)
which we will also refer to as the grand partition function. We stress that $`p_N(\stackrel{}{x},\beta ϵ)`$ is a generic partition function: its detailed form will depend on whether or not we allow crossing.

To proceed, it is useful to partition the DNA molecule in ds helices (with the strands attached to each other) and bubbles, sequences in which the two chains share just the first and the last base pairs; we also single out the contribution of the unzipped end of the two DSAWs, i.e. the part from the last contact to the end. In this way, the grand partition sum (4) can be expressed as:

$$Z(z,\beta ϵ,\beta g)=\frac{1}{1Dz\mathrm{exp}(\beta ϵ)}\frac{\mathrm{exp}(\beta ϵ)}{1\frac{\mathrm{exp}(\beta ϵ)}{1Dz\mathrm{exp}(\beta ϵ)}B(z)}S(z,\beta g),$$
(5)
where we have defined:

$`B(z)={\displaystyle \underset{N=2}{\overset{\mathrm{}}{}}}z^Nb_N,S(z,\beta g)=1+{\displaystyle \underset{N=1}{\overset{\mathrm{}}{}}}z^N{\displaystyle \underset{\stackrel{}{x}}{}}\left[c_{N1}(\stackrel{}{x})c_{N1}(\stackrel{}{0})\delta _{\stackrel{}{x},\stackrel{}{0}}\right]\mathrm{exp}(\beta \stackrel{}{g}\stackrel{}{x}).`$ (6)
In (6), $`b_N`$ is the number of $`2N`$-step bubbles, and $`c_N(\stackrel{}{x})`$ the partition function of two $`N`$-step DSAWs never touching each other and having their last sites at a mutual distance $`\stackrel{}{x}`$, and their initial sites at a relative distance $`\stackrel{}{e}_{i_0}\stackrel{}{e}_{j_0}`$ for some $`i_0j_0`$. By summing over all possibilities for initial conditions (see equation (8) below), $`b_N=_{ij;i,j=1}^Dc_{N2}(\stackrel{}{e}_i\stackrel{}{e}_j)`$, so that we need to find an explicit expression for the $`c_N`$’s. The equations they obey are:

$$c_{N+1}(\stackrel{}{x})=\underset{i,j=1}{\overset{D}{}}c_N(\stackrel{}{x}\stackrel{}{e}_i+\stackrel{}{e}_j)c_N(\stackrel{}{0})\underset{ij;i,j=1}{\overset{D}{}}\delta _{\stackrel{}{x},\stackrel{}{e}_i\stackrel{}{e}_j}.$$
(7)
Notice that in eq.(7) $`\stackrel{}{0}`$ acts as a sink or absorbing state, i.e. once the two walks join, they never leave. In this way, $`c_N(\stackrel{}{x})`$ with $`\stackrel{}{x}\stackrel{}{0}`$ counts the number of pairs of walks that never touch each other (and with the last monomers at a relative distance $`\stackrel{}{x}`$), while $`c_N(\stackrel{}{0})`$ counts the number of remaining pairs of walks that at some point come into contact and then remain stuck together. This last quantity plays the role of an arbitrary constant and has to be subtracted away from the final counting as in (6) (see  for another example in which this prescription was used in a $`D=3`$ example).

In order to simplify the presentation from now on we restrict our calculations on models w.o.c. (detailed calulations for this case and the outline of calculations for the simpler models w.c. are deferred to Appendix A). The initial conditions for such models are:

$`c_0(\stackrel{}{x})`$ $`=`$ $`\left({\displaystyle \underset{ij;i,j=1}{\overset{D}{}}}\delta _{\stackrel{}{x},\stackrel{}{e}_i\stackrel{}{e}_j}\right)\delta _{\stackrel{}{x},\stackrel{}{e}_{i_0}\stackrel{}{e}_{j_0}},`$ (8)
where $`\stackrel{}{e}_{i_0}`$ and $`\stackrel{}{e}_{j_0}`$ are the directions forbidden as described in the previous subsection.

By performing a Fourier and a discrete Laplace transform on eq.(7), and by using the relation between $`b_N`$ and $`c_N`$ given above, we can derive an explicit expression for the bubble generating function $`B(z)`$ (see eq. (A6)). The singularities of the second denominator of the grand partition function (5) are necessary to find the melting temperature. We find that the equations locating these two singularities are :

$`z_1`$ $`=`$ $`{\displaystyle \frac{1}{D^2}},`$ (9)
$`\mathrm{exp}(\beta ϵ)Dz_2`$ $`=`$ $`B\left(z_2\right).`$ (10)
The first singularity $`z_1`$ is leading to the usual random walk behavior in the absence of any interaction. The second singularity $`z_2\left(\beta ϵ\right)`$ is a function of the strength of the attractive interaction between the two strands, determining thus the behaviour in the native zipped phase. At $`g=0`$, the critical temperature for thermal melting is obtained when the two singularities above coincide ($`z_1=z_2`$). Unlike the models w.c., which have $`T_c=\mathrm{}`$ in $`D3`$, those w.o.c. have a finite critical temperature in any $`D2`$ (see Appendix A for further details).

As regards the force-dependent third factor in (5), it reads:

$$S(z,\beta g)=zS_{sing}(z,\beta g)+1zc(z,\stackrel{}{0}).$$
(11)
A third new singularity, depending on the external force, and fundamental in our calculations arises when computing:

$$S_{sing}(z,\beta g)_{[\pi ,\pi ]^D}\frac{d^D\stackrel{}{q}}{(2\pi )^D}\stackrel{}{c}(z,\stackrel{}{q})\underset{\stackrel{}{x}}{}\mathrm{exp}[(\beta \stackrel{}{g}i\stackrel{}{q})\stackrel{}{x}],$$
(12)
where $`\stackrel{}{c}(z,\stackrel{}{q})`$ and $`c(z,\stackrel{}{0})`$ are given in eqs. (A4,A5). Notice that in eq.(12) one can immediately compute the integrals on $`dq_3\mathrm{}dq_D`$, by using the well known identity $`2\pi \delta \left(q\right)=_x\mathrm{exp}\left(iqx\right)`$, so that one is left with the double integral:

$$_\pi ^\pi \frac{dq_2}{(2\pi )}_\pi ^\pi \frac{dq_1}{(2\pi )}\stackrel{}{c}(z,q_1,q_2,0,\mathrm{},0)\underset{x_1,x_2}{}\mathrm{exp}[(\beta giq_1)x_1+(\beta giq_2)x_2].$$
(13)
We now give details for the evaluation of the inner integral only, the outer one being equivalent. Defining $`h(q_1,q_2)\left(\stackrel{}{c}(z,q_1,q_2,0,\mathrm{},0)\right)^1`$, we aim at proving that the complex translation $`q_1q_1i\beta g`$ can be performed in the integral $`_\pi ^\pi \frac{dq_1}{2\pi }\frac{\mathrm{exp}\left(iq_1x_1\right)}{h(q_1,q_2)}`$. Notice that $`h`$ is a periodic function of $`q_1`$. We extend $`q_1`$ to complex numbers and consider the contour $`\gamma `$ as shown in Fig. 1. By the residue theorem, one can write:

$`{\displaystyle _\pi ^\pi }{\displaystyle \frac{dq_1}{2\pi }}{\displaystyle \frac{\mathrm{exp}\left(iq_1x_1\right)}{h(q_1,q_2)}}`$ $`+`$ $`{\displaystyle _0^{\beta g}}{\displaystyle \frac{idy}{2\pi }}{\displaystyle \frac{\mathrm{exp}\left(i\pi x_1+x_1y\right)}{h(\pi +iy,q_2)}}`$ (14)
$`{\displaystyle _\pi ^\pi }{\displaystyle \frac{dq_1}{2\pi }}{\displaystyle \frac{\mathrm{exp}\left(iq_1x_1\beta gx_1\right)}{h(q_1i\beta g,q_2)}}`$ $`+`$ $`{\displaystyle _{\beta g}^0}{\displaystyle \frac{idy}{2\pi }}{\displaystyle \frac{\mathrm{exp}\left(i\pi x_1+x_1y\right)}{h(\pi +iy,q_2)}}={\displaystyle \underset{q_0}{}}\mathrm{Res}(,q_0),`$ (15)
where the sum runs over all poles $`q_0`$ of $`\frac{\mathrm{exp}\left(iq_1x_1\right)}{h(q_1,q_2)}`$ (at fixed $`q_2`$) inside the contour $`\gamma `$ and $`\mathrm{Res}(,q_0)`$ is the residue of $``$ at $`q_0`$. Note that due to the periodicity of $`h`$ the second and fourth terms in the previous equation cancel and thus the complex translation can be performed provided that no pole of $``$ can be found inside the contour $`\gamma `$. An equivalent conclusion can be reached for the outer integral on $`dq_2`$. The condition of having no poles inside the contours of integration is satisfied as long as $`z<z_3(\beta g)`$, where:

$$z_3(\beta g)=\frac{1}{\left(D2\right)^2+2+2\mathrm{cosh}\left(2\beta g\right)+4\left(D2\right)\mathrm{cosh}\left(\beta g\right)}$$
(16)
is found by solving $`h(i\beta g,i\beta g)(1zf(i\beta \stackrel{}{g}))=0`$, with $`f(\stackrel{}{q})`$ defined in eq. (A2). In this case, we find consistently that $`z_3(\beta g)`$ is just the singularity in the resulting expression $`S(z,\beta g)=1+z\left[\stackrel{}{c}(z,i\beta \stackrel{}{g})c(z,\stackrel{}{0})\right]`$.

As $`z_3`$ is always smaller than $`z_1`$ (for $`g0`$), the singularity closest to the origin has to be determined between $`z_2`$, controlled by the attractive energy $`ϵ`$, and $`z_3`$, controlled by the pulling force $`g`$. If $`z_2<z_3`$, the DNA molecule is zipped, otherwise it is unzipped. The free energy per monomer $`f`$ and the average distance between the two ends (projected onto $`\widehat{g}(1,1,0,\mathrm{},0)`$), $`<x_g>`$, are defined as:

$`f\underset{N\mathrm{}}{lim}T{\displaystyle \frac{\mathrm{log}\left[_\stackrel{}{x}p_N(\stackrel{}{x},\beta ϵ)\mathrm{exp}\left(\beta \stackrel{}{g}\stackrel{}{x}\right)\right]}{N}}`$ (17)
$`<x_g>\underset{N\mathrm{}}{lim}<\widehat{g}\stackrel{}{x}>={\displaystyle \frac{f}{g}}N`$ (18)
These quantities read in the thermodynamic limit:

$`f`$ $`=T\mathrm{log}z_2\left(\beta ϵ\right);{\displaystyle \frac{<x_g>}{N}}=0g<g_c(T,ϵ)`$ (19)
$`f`$ $`=T\mathrm{log}z_3\left(\beta g\right);{\displaystyle \frac{<x_g>}{N}}=z_3\left(\beta g\right)\left[4\mathrm{sinh}\left(2\beta g\right)+4\left(D2\right)\mathrm{sinh}\left(\beta g\right)\right]g>g_c(T,ϵ)`$ (20)
where the critical force $`g_c(T,ϵ)`$ is found by imposing $`z_2\left(\beta ϵ\right)=z_3\left(\beta g\right)`$, as given in eqs. (10,16). The above equations show the existence of a first order transition at $`g=g_c(T,ϵ)`$, if $`g_c(T,ϵ)>0`$. The phase diagram for the models w.c. can be found exactly in the same way (see Appendix A for some details on this). It is interesting to notice that the singularity $`z_3(\beta g)`$ does not depend on whether or not we allow crossing, whereas $`z_2`$ in the w.c. case is different from the w.o.c. case. This is enough of course to make the whole phase diagram different for the two kinds of models, as will be discussed in section III.

We stress here that the main interest in using directed walks is that they are a sub-class of SAWs in the same dimension. However, qualitatively similar results can be obtained also with simple random walks (RWs), but in that case it would be not physically meaningful to forbid crossing as done above.

Experimentally, however, it should not be hard to set up an experiment closely related to the calculation we have performed, by stretching the ds molecule in one given direction before applying the external force.

## III Scaling laws for thermal melting and unzipping

Let us focus for concreteness on the result for the model in $`D=2`$ w.o.c., in order to analyze the scaling laws of the system. The phase diagram, as found also in , explicitly reads (see Fig. 2):

$$g_c(T,ϵ)=\frac{T}{2}\mathrm{cosh}^1\left[\frac{1}{2}\frac{1}{\sqrt{1\mathrm{exp}(\beta ϵ)}(1\mathrm{exp}(\beta ϵ))}1\right].$$
(21)
The model exhibits a first-order “unzipping” transition if we move at a fixed value of $`T<ϵ/\mathrm{log}\left(4/3\right)`$, as shown in (19). As is shown in Fig. 2, $`lim_{T0}g_c(T,ϵ)=\frac{ϵ}{2}`$ and $`g_c(T,ϵ)`$ attains its maximum at $`T=T_M0.9ϵ`$ where $`g_c(T_M)0.68ϵ>\frac{ϵ}{2}`$. The transition is second order at $`g=0`$; as in  we find that close to the melting point

$$<n>\frac{}{(\beta ϵ)}\mathrm{log}Z_N\tau ^1f\left(\tau N^{\varphi _t}\right),$$
(22)
where $`<n>`$ is the mean number of native contacts, with $`\tau =\left(TT_c\right)/T_c`$, and $`\varphi _t=1/2`$ for the crossover exponent at thermal melting in the $`D=2`$-case. The scaling function $`f(x)`$ behaves as usual

$$\{\begin{array}{cc}f(x)0\hfill & \mathrm{as}x+\mathrm{},\hfill \\ f(x)x\hfill & \mathrm{for}|x|1,\hfill \\ f(x)\left|x\right|^{1/\varphi _t}\hfill & \mathrm{for}x\mathrm{},\hfill \end{array}$$
(23)
such that $`<n>N^{\varphi _t}`$ at the transition.

The exact expressions for $`<n>`$ and for $`<x_g>`$ can be found by inverse Laplace transforms of the quantities $`Z(z,\beta g,\beta ϵ)`$, $`\frac{}{(\beta ϵ)}Z(z,\beta g,\beta ϵ)`$,$`\frac{}{(\beta g)}Z(z,\beta g,\beta ϵ)`$. However, if we only need scaling relations in the thermodynamic limit, we can use the discrete Tauberian theorem (see  and Appendix B), which relates the critical behaviour of a series to the asymptotic behaviour of its coefficients. By using this method, we find, in the vicinity of the unzipping mechanical transition:

$`<n>{\displaystyle \frac{h_1\left[N^{\varphi _1}\left(A_1\gamma +A_2\tau \right)\right]}{A_1\gamma +A_2\tau }},`$ (24)
$`<x_g>{\displaystyle \frac{h_2\left[N^{\varphi _2}\left(A_1\gamma +A_2\tau \right)\right]}{A_1\gamma +A_2\tau }},`$ (25)
where $`A_1(T_c,g_c),A_2(T_c,g_c)`$ are determined in such a way that $`\frac{\mathrm{d}g_c(T,ϵ)}{\mathrm{d}T}=\frac{A_2}{A_1}\frac{g_c}{T_c}`$, $`\gamma =\left(gg_c\right)/g_c`$, $`\tau =\left(TT_c\right)/T_c`$, and $`\varphi _1=\varphi _2=1`$ for the two crossover exponents, consistently with the fact that the unzipping transition is first-order (the explicit form of $`A_1(T_c,g_c),A_2(T_c,g_c)`$ is worked out in Appendix B). Note that $`A_1>0`$, whereas $`A_2`$ is negative in the re-entrant part of the transition curve. The scaling functions $`h_{1,2}(x)`$ behave in a similar way to $`f(x)`$, and thus we get that $`<n>N^{\varphi _1}`$ and $`<x_g>N^{\varphi _2}`$ are both extensive at the transition point. The physical interpretation is that a macroscopic portion of the chain is still in the double stranded state, but the rest of the chain is unzipped; these are just the two phases coexisting at the first order transition. This result will be used later to justify the use of Y-shaped configurations.

It is instructive to compare this phase diagram with that of the same $`D=2`$-case when we allow crossing. The transition line (plotted in Fig. 2), obtained from $`z_2=z_3`$ (see Appendix A) is:

$$g_c(T,ϵ)=T\mathrm{tanh}^1\left[1\mathrm{exp}(\beta ϵ)\right]=\frac{ϵ}{2}+\frac{1}{2\beta }\mathrm{log}\left[2\mathrm{exp}(\beta ϵ)\right],$$
(26)
Both the models behave similarly near $`T=0`$, namely they yield a transition to the cold denaturated state described in . When we move at a constant force $`g`$ such that $`g_c(0)<g<Max_T\left[g_c(T,ϵ)\right]`$, at low enough temperatures the molecule is unzipped, and it zips by increasing $`T`$. This feature of the phase diagrams seems at first sight paradoxical, since one would expect the critical force to decrease monotonically as the temperature is increased. The physical explanation of this result will be given below. In the model w.c., moreover, $`g_c`$ always increases, approaching $`ϵ`$ as $`T\mathrm{}`$, and the two strands remain zipped for every temperature when $`g<ϵ`$. In the model w.o.c., instead, the two strands unzip again by further increasing the temperature. As regards scaling, on the other hand, the two models are identical so that we can say that the effect of forbidding crossing is irrelevant in the renormalization group sense, but has a dramatic effect on the form of the critical line.

Let us now discuss the results that we have obtained in higher dimensions (the critical lines for three and four dimensional models are shown in Fig. 3). Common to all models is the cold unzipping transition found in the $`D=2`$ case. Moreover, the scaling laws at the unzipping transition do not depend on dimensionality in the crossover exponents. On the other hand, the thermal melting has dimensionality-dependent critical exponents (see ). Consequently, the behaviour of the critical line $`g_c(T,ϵ)`$ near the end point $`T=T_c`$ also depends on $`D`$. As $`TT_c^{}`$ the result that we find in generic dimension is:

$`{\displaystyle \frac{g_c(T,ϵ)}{T}}\{\begin{array}{cc}ϵ/Tϵ/T_c(D)\hfill & D=2,4\hfill \\ \mathrm{exp}\left[\frac{a}{(ϵ/Tϵ/T_c(D))}\right]\hfill & D=3\hfill \\ \left[ϵ/Tϵ/T_c(D)\right]^{1/2}\left[\mathrm{log}(ϵ/Tϵ/T_c(D))\right]^{1/2}\hfill & D=5\hfill \\ \left[ϵ/Tϵ/T_c(D)\right]^{1/2}\hfill & D>5\hfill \end{array},`$ (31)
where $`a`$ is a constant, and $`T_c=\mathrm{}`$ in $`D=2,3`$ for models w.c.

## IV Physical explanation of re-entrance

### A Re-entrance in lattice models

To obtain some physical insight into our analytical results, we compute the free energy of a Y-shaped molecule of DNA with a force pulling at the extremities (see Fig. 4). In this way, we neglect all the configurations with bubbles. Such an approximation is valid at low T and yields the exact result in the limit $`T0`$. Since the configurations of the unzipped part are weighted by $`\mathrm{exp}\left(\beta \stackrel{}{g}\stackrel{}{x}\right)`$, in this limit only the completely stretched configuration will contribute to the free energy for the unzipped part of the Y. In the limit of low temperature, the free energy is then:

$$F(m,N)(Nm)(ϵ+T\mathrm{log}\mu )2gm,$$
(32)
where $`m`$ is the number of monomers in the unzipped part, $`N`$ is the total number of bases and $`\mu `$ is the connective constant of a single walk. Note that eq. (32) (with $`\stackrel{}{g}=g(1,0,\mathrm{},0)`$) is valid also in the case of two self-avoiding walks (SAWs), since the connective constant of SAWs constrained to avoid the fully stretched unzipped part, is the same as for ordinary SAWs  (for a more detailed discussion on SAWs we refer to Appendix C). By minimizing with respect to $`m`$, we find the critical force $`g_c(T,ϵ)`$ such that (32) is minimum for $`m=0`$ when $`g<g_c(T,ϵ)`$, and for $`m=N`$ when $`g>g_c(T,ϵ)`$:

$$g_c(T,ϵ)\stackrel{T0}{}\frac{1}{2}(ϵ+T\mathrm{log}\mu ).$$
(33)
This is indeed what we have found in our calculations, and confirms that re-entrance is a robust feature of lattice models, not depending on dimensionality or self-avoidance. In other words, (33) means that, at low $`T`$, it is more difficult to open a dsDNA helix as $`T`$ is increased, because the energy gain obtained through the unzipping is more than compensated by the entropy loss, since there is only one possible completely stretched configuration versus $`\mu ^{Nm}`$ possibilities for the double stranded portion of the chain (we will see in the following section that the entropy loss in the continuum space exhibits a power law correction). For high enough $`T`$, on the other hand, also other configurations will contribute to the open portion, increasing its entropy, and the energy gain will eventually favour the unzipped state, except that in $`D=2`$ w.c., where the presence of bubble enhances the entropy of the native portion at any $`T`$, as explained below.

We have also calculated the phase diagram obtained by considering only the Y configuration: this amounts to putting the generating function for bubbles $`B(z)=0`$ in (5). As expected, the Y-approximation gives the exact behaviour (33) in the limit $`T0`$, whereas it gets wronger and wronger as $`T`$ is increased (see Fig. 2 where we plot the two-dimensional case explicitly given by $`g_c(T,ϵ)=\frac{T}{2}\mathrm{cosh}^1\left[\mathrm{exp}(\beta ϵ)1\right]`$ ). For the critical line near the melting point, this approximation yields $`g_c(T,ϵ)(T_cT)^{1/2}`$ for all $`D`$, which from (31) is correct only for $`D>5`$ ($`d>4`$). This is because such an approach neglects bubbles, which are expected to be more and more relevant as $`d`$ decreases (for $`d>d_c=4`$ once the two walks have splitted away, they basically never meet again) and $`T`$ increases. In this view, it is clear that the puzzling behaviour found in the $`D=2`$ model with crossing, where $`g_c`$ increases for all $`T`$, is due to the growing presence of bubbles: before DNA can be unzipped, we have to disentangle all the bubbles that are forming, which is harder and harder as $`T\mathrm{}`$.

### B Re-entrance in the continuum

We now prove that also discrete chains in the continuum space undergo a cold denaturation for $`T0`$. Let us consider two $`N`$-monomer chains in $`𝐑^𝐝`$ with no constraint on the directedness, with a force $`\stackrel{}{g}`$ pulling at the extremities and with constant unitary distance between consecutive monomers. As $`\beta \mathrm{}`$, bubbles can be neglected and only Y-shaped configurations contribute to the partition sum. The partition function for a Y-configuration then reads:

$$Z_N(\beta ϵ,\beta g)=\underset{m=0}{\overset{N1}{}}Z_{Nm}^z(\beta ϵ)Z_{2m}^u(\beta g),$$
(34)
where $`Z_{Nm}^z(\beta ϵ)=_{𝐑^𝐝}d^d\stackrel{}{x}_1\mathrm{}d^d\stackrel{}{x}_{Nm1}\delta \left(|\stackrel{}{x}_1|1\right)\mathrm{}\delta \left(|\stackrel{}{x}_{Nm1}|1\right)\mathrm{exp}((Nm)\beta ϵ)`$ is the partition function of the zipped portion of the strands, and $`Z_{2m}^u(\beta g)=_{𝐑^𝐝}d^d\stackrel{}{y}_1\mathrm{}d^d\stackrel{}{y}_{2m}\delta \left(|\stackrel{}{y}_1|1\right)\mathrm{}\delta \left(|\stackrel{}{y}_{2m}|1\right)\mathrm{exp}\left(\beta \stackrel{}{g}_{i=1}^{2m}\stackrel{}{y}_i\right)`$ is that of the unzipped end. Chain discreteness is crucial in ensuring the validity of eq.(34) for $`\beta \mathrm{}`$, since when one bubble has formed, the length of the stretched portion of the chain can be at most $`_{i=1}^{2m1}\stackrel{}{y}_i=2(m1)\frac{\stackrel{}{g}}{|\stackrel{}{g}|}`$. The integrals in (34) can be performed in any dimensions, and the evaluation of $`Z_{2m}^u(\beta g)`$ involves the modified Bessel function of the first kind of order $`\frac{d2}{2}`$. In particular, $`Z_{2m}^u(\beta g)\stackrel{\beta \mathrm{}}{}\left(\frac{2\pi }{\beta g}\right)^{m(d1)}\mathrm{exp}(2m\beta g)`$, and there is a power-law correction with respect to the lattice result. We find that $`Z_N(\beta ϵ,\beta g)\stackrel{\beta \mathrm{}}{}_{m=0}^{N1}\mathrm{\Omega }_d^{Nm1}\mathrm{exp}((Nm)ϵ)Z_{2m}^u(\beta g)`$, where $`\mathrm{\Omega }_d`$ is the surface area of the unit sphere in $`d`$ dimensions. From this, the critical force is easily found:

$$g_c(T,ϵ)\stackrel{T0}{}\frac{ϵ}{2}\frac{d1}{2}T\mathrm{log}\frac{T}{ϵ}\frac{T}{2}\mathrm{log}\left[2^{2d3}\pi ^{\frac{d}{2}1}\mathrm{\Gamma }\left(\frac{d}{2}\right)\right],$$
(35)
where $`\mathrm{\Gamma }`$ is the Euler Gamma function. The critical line $`g_c(T,ϵ)`$ increases at low $`T`$ and re-entrance is present also in the continuum, even enhanced with respect to the lattice case. The leading term $`T\mathrm{log}T`$ in equation (35) (due to the power-law correction in $`Z_{2m}^u(\beta g))`$ is indeed not present in the $`d=1`$ case, when one correctly recovers equation (33) for a lattice random walk with $`\mu =2`$.

We finally wish to discuss the relation of the present treatment with previous work () done on the unzipping of homo-DNA in the limit of a continuum chain. The ds molecule with the pulling force had been described by means of an effective hamiltonian, which, apart from an irrelevant center of mass term, reads:

$$H=_0^N𝑑n\left(\frac{Td}{b^2}\left(\frac{d\stackrel{}{r}}{dn}\right)^2+V\left(\stackrel{}{r}(n)\right)\stackrel{}{g}\frac{d\stackrel{}{r}}{dn}\right),$$
(36)
where $`\stackrel{}{r}`$ is the relative separation between the strands, $`b`$ is the effective Kuhn length of single-stranded DNA and $`V`$ is a realistic short range attractive potential, namely a delta function or a potential well (the two cases should be equivalent according to the standard quantum theory as long as $`𝑑\stackrel{}{r}V(\stackrel{}{r})`$ is the same). The system described by (36) is equivalently represented by a quantum system with a non-hermitian hamiltonian. One finds  that there is a first order unzipping transition when the force reaches the critical value:

$$g_c(T,ϵ)=\sqrt{\frac{4ϵ_0(T)Td}{b^2}},$$
(37)
where $`ϵ_0(T)`$ is the ground state energy of the quantum hamiltonian obtained from (36) when $`\stackrel{}{g}=\stackrel{}{0}`$. In the quantum mechanical system, for $`g<g_c(T,ϵ)`$ the ground state is a bound state, while for $`g>g_c(T,ϵ)`$ the spectrum of (36) is continuous.

Let us focus on the $`d=1`$ case, corresponding to our two-dimensional models without crossing constraints. It is well known that at $`g=0`$ both a symmetric square well and a delta function always have at least one bound state, meaning that DNA remains ds at all $`T`$. In Fig. 5 we show the critical force for the two potentials (obtained simply from (37)). It can be seen from the figure that at high $`T`$ both the potentials saturate towards the same limit, as in our calculation for the model w.c. This is expected since the delta potential case can be recovered from our equation (2) in the limit of both continuum space and time (see Fig. 4). At low $`T`$, however, $`g_c(T,ϵ)\left(\frac{T}{\mathrm{\Delta }}\right)^{1/2}`$ for a square well of width $`\mathrm{\Delta }`$, and is constant for the delta potential, so that re-entrance is present only in the former case, but with a behaviour different from that found in the lattice. This already signals that the low temperature behaviour of the solution obtained through the quantum mapping is rather unphysical. This is due to the fact that in the limit of continuum chain the chain constraint is modelled in a ‘soft’ way, by using a harmonic potential between consecutive beads along the chain. This interaction is usually assumed to be entropic (as in eq. (36)) and thus effectively vanishes at $`T=0`$. In that limit eq. (36) is describing a set of $`N`$ ‘free’ particles moving in an effective potential determined by $`V`$ and $`g`$. Consequently, the quantum mapping is expected to describe the system well except in the low temperature limit. Indeed, the re-entrance found for the square well potential is due to the unboundedness of the potential felt by the “free particles” as soon as $`g0`$. As $`T0`$, fluctuations are not important, and the particles stay close to the minimum of such potential, which is $`\mathrm{}`$ for a force whatever small but finite. This is enough to unzip the strands .

As Zhou has observed in , one can improve this model by placing a hard core either inside the potential well or at a finite distance from the delta (this is analogous in spirit to introducing the crossing constraint in our lattice models). Once again, one finds that for T near $`T_c`$, which is now finite, both potentials predict $`g_c(T,ϵ)\stackrel{TT_c}{}(T_cT)`$ as in the lattice model (see eq.(31)), whereas for low $`T`$ the delta shows no re-entrance and the potential well displays a $`T^{1/2}`$ behaviour. We note once again that the low temperature behavior of such ‘quantum’ models is an artifact, when considered as polymer chain models, due to the ‘soft’ enforcing of the chain constraint. In this respect, discrete chain models in the continuum space are more realistic, either with constant bond length or with harmonic springs between consecutive beads. In the first case we have indeed proved (see eq. (35)) the existence of cold mechanical denaturation. In the second case it has been shown to occur by means of numerical simulations .

## V Conclusions

To conclude, we have studied simple models for DNA mechanical unzipping induced by a pulling force. By using analytical techniques, the relevance of forbidding the crossing of the two strands and the role played by bubble formation in the denaturation process has been discussed throughout the whole force-temperature plane. The scaling properties of the system along the transition line have also been determined. The existence of a cold mechanical denaturation for a finite range of forces at low enough temperatures has been exactly proved in a general case for both lattice and continuum space models. The fundamental role of chain discreteness has been emphasized in comparison with related models studied by means of quantum mapping techniques. Even if we neglected many important effects, such as base pairs heterogeneity, intrinsic helicity, wrong base pairing, and different stacking energy for the natured and denatured state, it would be interesting to experimentally test this prediction stemming from such a simple model.

Furthermore, on the basis of preliminary results (analytical and numerical) we anticipate that even in the presence of quenched randomness in the contact potential our models display a re-entrant transition, so in this sense it is really a “universal” feature.

Acknowledgements: We would like to thank F. Seno and S.M. Bhattacharjee for illuminating discussions. DM also acknowledges INFM funding. This work was supported by MURST grant.

## A Computation of $`B(z)`$ in models of DSAWs without crossing and of the phase diagram for models with crossing

We here first work out the details to obtain an explicit expression for $`B(z)`$, defined in eq.(6). Starting from equations (7), we can perform a Fourier and a discrete Laplace transform to obtain:

$`\stackrel{}{c}(z,\stackrel{}{q}){\displaystyle \underset{N=0}{\overset{\mathrm{}}{}}}z^N{\displaystyle \underset{\stackrel{}{x}}{}}\mathrm{exp}(i\stackrel{}{q}\stackrel{}{x})c_N(\stackrel{}{x})={\displaystyle \frac{1}{1zf(\stackrel{}{q})}}(\stackrel{}{c_0}(\stackrel{}{q})zf_0(\stackrel{}{q})c(z,\stackrel{}{0})),`$ (A1)
where we have made the following definitions:

$`f_0(\stackrel{}{q})2{\displaystyle \underset{i<j;i,j=1}{\overset{D}{}}}\mathrm{cos}\left(q_iq_j\right);f(\stackrel{}{q})D+f_0(\stackrel{}{q})`$ (A2)
$`c(z,\stackrel{}{x}){\displaystyle \underset{N=0}{\overset{\mathrm{}}{}}}z^Nc_N(\stackrel{}{x});\stackrel{}{c_0}(\stackrel{}{q}){\displaystyle \underset{\stackrel{}{x}}{}}\mathrm{exp}(i\stackrel{}{q}\stackrel{}{x})c_0(\stackrel{}{x})`$ (A3)
Using $`c(z,\stackrel{}{0})=_{[\pi ,\pi ]^D}\frac{d^D\stackrel{}{q}}{(2\pi )^D}\stackrel{}{c}(z,\stackrel{}{q})`$, equation (A1) is easily solved since, by permuting variables inside the resulting integrals, it is possible to see that they do not depend on the particular step which we forbid. We obtain:

$`\stackrel{}{c}(z,\stackrel{}{q})`$ $`=`$ $`{\displaystyle \frac{1}{1zf(\stackrel{}{q})}}{\displaystyle \frac{\stackrel{}{c_0}(\stackrel{}{q})zW_1(z)[A(D)f_0(\stackrel{}{q})\stackrel{}{c_0}(\stackrel{}{q})]}{1+zW_1(z)}},`$ (A4)
$`c(z,\stackrel{}{0})`$ $`=`$ $`A(D){\displaystyle \frac{W_1(z)}{1+zW_1(z)}},`$ (A5)
where $`A(D)=1\frac{1}{D(D1)}`$ is the reduction factor due to the crossing constraint. As a result, using $`B(z)=z^2_{ij;i,j=1}^Dc(z,\stackrel{}{e}_i\stackrel{}{e}_j)`$, we obtain

$$B\left(z\right)=A(D)\frac{z^2W_2(z)}{1+z_2W_1(z)},$$
(A6)
where

$$W_1(z)=_{[\pi ,\pi ]^D}\frac{d^D\stackrel{}{q}}{(2\pi )^D}\frac{f_0(\stackrel{}{q})}{1zf(\stackrel{}{q})},W_2(z)=_{[\pi ,\pi ]^D}\frac{d^D\stackrel{}{q}}{(2\pi )^D}\frac{f_0^2(\stackrel{}{q})}{1zf(\stackrel{}{q})}.$$
(A7)
The singularity $`z_1=1/D^2`$, leading to the usual random walk behavior, comes in when evaluating the above integrals; their denominator becomes negative as $`\stackrel{}{q}0`$ for $`z>z_1`$.

We give here also a brief outline of the calculation necessary to find out the phase diagram for the models of Section II of DSAWs when crossing is allowed. These models are simpler to solve than the corresponding models w.o.c., since it is not necessary to partition the molecule of DNA as done in the text in eq.(5), which allows us to write a more explicit formula for the grand partition function. Indeed we can simply start from the recursion relations (2) and, after Fourier and Laplace transform, find the expression:

$$\stackrel{}{p}(z,\stackrel{}{q},\beta ϵ)=\frac{1}{1zf(\stackrel{}{q})}\frac{1}{1\left(1\mathrm{exp}\left(\beta ϵ\right)\right)W_0(z)},$$
(A8)
for the quantity $`\stackrel{}{p}(z,\stackrel{}{q},\beta ϵ)`$, defined as:

$$\stackrel{}{p}(z,\stackrel{}{q},\beta ϵ)=\underset{N=0}{\overset{\mathrm{}}{}}z^N\underset{\stackrel{}{x}}{}\mathrm{exp}\left(i\stackrel{}{q}\stackrel{}{x}\right)p_N(\stackrel{}{x},\beta ϵ)$$
(A9)
with $`p_N(\stackrel{}{x},\beta ϵ)`$ defined as in eq.(1) and

$$W_0=_{[\pi ,\pi ]^D}\frac{d^D\stackrel{}{q}}{(2\pi )^D}\frac{1}{1zf(\stackrel{}{q})}.$$
(A10)
The final form of the grand partition function, for DSAWs w.c. in the presence of a pulling force, is:

$$Z(z,\beta ϵ,\beta g)=\frac{1}{1z/z_3(\beta g)}\frac{1}{1\left(1\mathrm{exp}\left(\beta ϵ\right)\right)W_0(z)},$$
(A11)
with $`z_3(\beta g)`$ defined as in eq.(16). Alternatively, one can proceed exactly as in Section II, but with initial conditions

$`c_0(\stackrel{}{x})`$ $`=`$ $`{\displaystyle \underset{ij;i,j=1}{\overset{D}{}}}\delta _{\stackrel{}{x},\stackrel{}{e}_i\stackrel{}{e}_j}.`$ (A12)
The only resulting difference would be to have $`A(D)=1`$ in eqs.(A4) and (A6); after a bit of algebra one can convince himself that the partition function found by recollecting the different factors in eq.(5) is the same as that found in eq.(A11). In doing this, one needs to compute the relations between the integrals $`W_0(z)`$, $`W_1(z)`$, $`W_2(z)`$, which may be done as, e.g., in . For example, the bubble generating function for both kind of models could also be simply expressed as

$$B(z)=A(D)\left(1Dz1/W_0(z)\right).$$
(A13)
Changing the dimension dependent factor $`A(D)`$ from $`A(D)=1`$, when crossing is allowed, to $`A(D)=1\frac{1}{D(D1)}`$, when crossing is forbidden, is enough for the melting transition temperature to become finite for $`D3`$.

## B Tauberian theorem and its applications

The discrete Tauberian theorem (see ) states that the relations

$$c\left(z\right)\underset{N=0}{\overset{\mathrm{}}{}}c_Nz^N\stackrel{zz^{}}{}\left(\frac{1}{1z/z^{}}\right)^\rho L\left(\frac{1}{1z/z^{}}\right),$$
(B1)
and

$$c_N\stackrel{N\mathrm{}}{}\frac{N^{\rho 1}}{\left(z^{}\right)^N\mathrm{\Gamma }(\rho )}L(N),$$
(B2)
are equivalent provided that $`\rho >0`$, $`\{c_N\}`$ is a positive monotonic sequence and $`L`$ is “slowly varying” in the sense specified in .

We apply this theorem here first to find the scaling laws of the unzipping transition to give an example of how the relations reported in section III can be recovered. When approaching the critical line at $`g=g_c`$ from the unzipped phase, the smallest singularity is $`z_3\left(\beta g\right)`$, as defined in eq. (16). The leading behavior as $`zz_3^{}`$ is:

$`Z(z,\beta ϵ,\beta g)`$ $`\stackrel{zz_3^{}}{}`$ $`{\displaystyle \frac{1}{\mathrm{exp}(\beta ϵ)Dz_3B(z_3)}}{\displaystyle \frac{1}{1z/z_3}},`$ (B3)
$`{\displaystyle \frac{Z(z,\beta ϵ,\beta g)}{(\beta g)}}`$ $`\stackrel{zz_3^{}}{}`$ $`{\displaystyle \frac{1}{\mathrm{exp}(\beta ϵ)Dz_3B(z_3)}}{\displaystyle \frac{1}{\left(1z/z_3\right)^2}}{\displaystyle \frac{z_3/(\beta g)}{z_3(\beta g)}},`$ (B4)
$`{\displaystyle \frac{Z(z,\beta ϵ,\beta g)}{(\beta ϵ)}}`$ $`\stackrel{zz_3^{}}{}`$ $`{\displaystyle \frac{\mathrm{exp}(\beta ϵ)}{\left(\mathrm{exp}(\beta ϵ)Dz_3B(z_3)\right)^2}}{\displaystyle \frac{1}{1z/z_3}},`$ (B5)
where we have neglected inessential factors when $`gg_c`$. Consequently, the application of the Tauberian theorem yields:

$`Z_N(\beta ϵ,\beta g)`$ $`\stackrel{N\mathrm{}}{}`$ $`{\displaystyle \frac{1}{\mathrm{exp}(\beta ϵ)Dz_3B(z_3)}}z_3^N,`$ (B6)
$`{\displaystyle \frac{Z_N(\beta ϵ,\beta g)}{(\beta g)}}`$ $`\stackrel{N\mathrm{}}{}`$ $`{\displaystyle \frac{N}{\mathrm{exp}(\beta ϵ)Dz_3B(z_3)}}z_3^N{\displaystyle \frac{z_3/(\beta g)}{z_3(\beta g)}},`$ (B7)
$`{\displaystyle \frac{Z_N(\beta ϵ,\beta g)}{(\beta ϵ)}}`$ $`\stackrel{N\mathrm{}}{}`$ $`{\displaystyle \frac{\mathrm{exp}(\beta ϵ)}{\left(\mathrm{exp}(\beta ϵ)Dz_3B(z_3)\right)^2}}z_3^N.`$ (B8)
The canonical partition function $`Z_N`$ can be expressed in terms of previously defined quantities as $`Z_N(\beta ϵ,\beta g)=_\stackrel{}{x}p_N(\stackrel{}{x},\beta ϵ)\mathrm{exp}(\beta \stackrel{}{g}\stackrel{}{x})`$ (see eqs.(1) and (4)).

Note that at the critical line $`z_2\left(\beta ϵ\right)=z_3\left(\beta g_c\right)`$, implying that the $`ϵ`$-dependent denominator in the above equations vanishes, as the critical line is approached (that is $`\gamma \left(gg_c\right)/g_c1`$ and $`\tau \left(TT_c\right)/T_c1`$). In that limit the order parameters obey

$`<n>{\displaystyle \frac{Z_N/(\beta ϵ)}{Z_N}}\stackrel{N\mathrm{}}{}{\displaystyle \frac{\mathrm{exp}(\beta ϵ)}{A_1\gamma +A_2\tau }},`$ (B9)
$`<x_g>{\displaystyle \frac{Z_N/(\beta g)}{Z_N}}\stackrel{N\mathrm{}}{}{\displaystyle \frac{z_3/(\beta g)}{z_3(\beta g_c)}}N.`$ (B10)
The coefficients $`A_1(T_c,g_c),A_2(T_c,g_c)`$ determine the normal direction in the ($`\mathrm{log}(T),\mathrm{log}(g)`$) plane, and are reported below:

$`A_1(T_c,g_c)=`$ $`\beta _cg_c\left(D+{\displaystyle \frac{B\left(z_3(\beta _cg_c)\right)}{z_3}}\right)z_3^2\left(\beta _cg_c\right)\left[4\mathrm{sinh}\left(2\beta _cg_c\right)+4\left(D2\right)\mathrm{sinh}\left(\beta _cg_c\right)\right],`$ (B11)
$`A_2(T_c,g_c)=`$ $`\beta _cg_c\left(D+{\displaystyle \frac{B\left(z_3(\beta _cg_c)\right)}{z_3}}\right)z_3^2\left(\beta _cg_c\right)\left[4\mathrm{sinh}\left(2\beta _cg_c\right)+4\left(D2\right)\mathrm{sinh}\left(\beta _cg_c\right)\right]+`$ (B13)
$`\beta _cϵ\mathrm{exp}\left(\beta _cϵ\right).`$
As expected, it may be easily checked that, whatever direction we choose to approach the critical curve, the physical quantities as (B9) only depend on the projection of this direction on the normal to the critical line. To show this, one only needs first to observe that criticality is achieved through the vanishing of the denominator in (B9), and then to apply the implicit function theorem to this denominator.

Recently, some authors () have suggested, as an alternative order parameter for the unzipping transition, the number $`m`$ of monomers from the last contact to the end. By using the same tools as before, we can find some quantities of interest such as the probability distribution of having $`m`$ monomers “liberated” (in the equations below we suppose 1$`<<m<<`$N). In general:

$$P(m)=\frac{_W\delta _{n_{u.z.},m}p_W}{_Wp_W}$$
(B14)
where W are the pairs of directed walks, $`p_W`$ is their Boltzmann weight and $`n_{u.z.}`$ is the number of “liberated” monomers of the configuration. The possible enforcing of the crossing constraint does not affect the following equations, simply shifting the critical melting temperature $`T_c`$.

At zero force and $`T<T_c`$, in the native state:

$$P(m)\mathrm{exp}\left[m\mathrm{log}(z_1/z_2)\right]m^{\alpha (d)},$$
(B15)
where $`m^{\alpha (d)}=m^{1/2}`$ in $`d=`$1, $`\left(\mathrm{log}(m)\right)^1`$ in $`d=`$2 and constant for $`d>`$2. We recall that $`z_1=1/D^2`$, with $`D=d+1`$. At criticality:

$$P(m)(Nm)^{\varphi _t1}m^{\alpha (d)},$$
(B16)
where in the first factor there are logarithmic corrections in $`d=2,4`$ to the power-law behavior controlled by the crossover thermal exponent $`\varphi _t`$ . Above $`T_c`$ we have:

$$P(m)(Nm)^{\beta (d)}m^{\alpha (d)}.$$
(B17)
with $`\beta (d)=|1d/2|1`$ if $`d`$2, in $`d=`$2 $`m^{\beta (d)}=1/(N\mathrm{log}^2(N))`$. At force different from zero, in the zipped phase, we get:

$$P(m)\mathrm{exp}\left(m\mathrm{log}(z_3/z_2)\right),$$
(B18)
while in the unzipped phase:

$$P(m)\mathrm{exp}\left((Nm)\mathrm{log}(z_2/z_3)\right),$$
(B19)
if $`T<T_c`$, and

$$P(m)\mathrm{exp}\left((Nm)\mathrm{log}(z_1/z_3)\right),$$
(B20)
if $`T>T_c`$. Just at criticality, the probability distribution is flat.

For the Y-shaped configuration, we get, at force different from zero:

$`P(m)`$ $``$ $`\mathrm{exp}\left((Nm)\mathrm{log}\left(z_y/z_3\right)\right)\mathrm{u}.\mathrm{z}.\mathrm{phase},`$ (B21)
$`P(m)`$ $``$ $`\mathrm{exp}\left(m\mathrm{log}\left(z_3/z_y\right)\right)\mathrm{z}.\mathrm{phase},`$ (B22)
$`P(m)`$ $``$ $`\mathrm{const}.\mathrm{at}\mathrm{criticality},`$ (B23)
where $`z_y=\mathrm{exp}\left(\beta ϵ\right)/\mu `$ is the singularity controlling the free energy per monomer in the Y-zipped phase, and $`\mu =D`$ the connective constant of a single directed walk. At zero force we get:

$`P(m)`$ $``$ $`\mathrm{exp}\left((Nm)\mathrm{log}\left(z_y/z_1\right)\right)m^{\alpha (d)}\mathrm{denatured}\mathrm{phase},`$ (B24)
$`P(m)`$ $``$ $`\mathrm{exp}\left(m\mathrm{log}\left(z_1/z_y\right)\right)m^{\alpha (d)}\mathrm{native}\mathrm{phase},`$ (B25)
$`P(m)`$ $``$ $`m^{\alpha (d)}\mathrm{at}\mathrm{criticality}.`$ (B26)
Let us give the details in a simple case, to get e.g. eqns.(B21,B22,B23). One has:

$$P(m)=\frac{Z_{Nm}^{PS}Z_m^{free}}{Z_N}$$
(B27)
where $`Z_{Nm}^{PS}`$ is the partition function of two $`\left(Nm\right)`$-monomer chains with the end points in common (this is the original Poland-Sheraga model), $`Z_m^{free}`$ is the partition of two $`m`$-monomer chains with no base pair in contact, and $`Z_N`$ is given by eq.(1). We have, for the Y at non zero force, $`Z_{Nm}^{PS}=z_y^{(Nm)}`$, $`Z_m^{free}z_3^m`$, and

$`Z_N`$ $``$ $`z_y^N\mathrm{native}\mathrm{phase},`$ (B28)
$`Z_N`$ $``$ $`z_3^N\mathrm{above}\mathrm{criticality}.`$ (B29)
## C Y configurations with SAWs

In this Appendix we extend the analysis of section IV to a Y configuration where the two chains are generic SAWs and not directed walks. We take here $`\stackrel{}{g}=g(1,0,\mathrm{},0)`$.

Let us call $`c_N(\stackrel{}{0},\stackrel{}{x})`$ the number of $`N`$-step SAWs ending in $`\stackrel{}{x}`$ and starting in $`\stackrel{}{0}`$, and $`c_{N,\stackrel{}{g}}=_\stackrel{}{x}c_N(\stackrel{}{0},\stackrel{}{x})\mathrm{exp}\left(\beta \stackrel{}{g}\stackrel{}{x}\right)`$ (so that $`c_Nc_{N,\stackrel{}{g}=\stackrel{}{0}}`$). As a first step, we need to prove the existence of a connective constant $`\mu _\stackrel{}{g}`$ for a single SAW in the presence of a pulling force i.e. $`c_{N,\stackrel{}{g}}\mu _\stackrel{}{g}^N`$ for large $`N`$. The existence of a connective constant then follows  from the inequality:

$`{\displaystyle \underset{\stackrel{}{x}}{}}c_N(\stackrel{}{0},\stackrel{}{x})\mathrm{exp}\left(\beta \stackrel{}{g}\stackrel{}{x}\right){\displaystyle \underset{\stackrel{}{x},\stackrel{}{x}_1}{}}c_{N_1}(\stackrel{}{0},\stackrel{}{x}_1)c_{NN_1}(\stackrel{}{x}_1,\stackrel{}{x})\mathrm{exp}\left(\beta \stackrel{}{g}\stackrel{}{x}_1\right)\mathrm{exp}\left(\beta \stackrel{}{g}(\stackrel{}{x}\stackrel{}{x}_1)\right),`$ (C1)
which implies the subadditivity of $`c_N`$ . From the connective constants of RWs and DSAWs, we can estabilish the following bounds (recall that here $`\stackrel{}{g}=g(1,0,\mathrm{},0)`$) for the connective constant of a SAW in $`d`$-dimensions:

$$d1+\mathrm{exp}\left(\beta |\stackrel{}{g}|\right)\mu _\stackrel{}{g}2(d1)+2\mathrm{cosh}\left(\beta |\stackrel{}{g}|\right).$$
(C2)
Let us now introduce the canonical partition function. It reads:

$$Z_N=\underset{m=0}{\overset{N1}{}}\mathrm{exp}\left(\beta (Nm)ϵ\right)c_{Nm1}\underset{\stackrel{}{x}_1,\stackrel{}{x}_2}{}c_{2m}^{}(\stackrel{}{x}_1,\stackrel{}{x}_2)\mathrm{exp}\left(\beta \stackrel{}{g}(\stackrel{}{x}_2\stackrel{}{x}_1)\right),$$
(C3)
where $`c_{2m}^{}(\stackrel{}{x}_1,\stackrel{}{x}_2)`$ counts the SAWs which do not cross the first $`Nm`$ attached monomers of the Y and $`\stackrel{}{x}_1,\stackrel{}{x}_2`$ are the end-points of the strands.

An upper bound is found if we allow the opened part of the Y to cross the first $`Nm`$ zipped monomers:

$$Z_N^{upper}=\underset{m=0}{\overset{N1}{}}\mathrm{exp}\left(\beta (Nm)ϵ\right)c_{Nm1}\underset{\stackrel{}{x}_1,\stackrel{}{x}_2}{}c_{2m}(\stackrel{}{x}_1,\stackrel{}{x}_2)\mathrm{exp}\left(\beta \stackrel{}{g}(\stackrel{}{x}_2\stackrel{}{x}_1)\right).$$
(C4)
The grand canonical partition function defined from (C4) displays singularities in $`zz_1=\frac{1}{\mu _\stackrel{}{g}^2}`$ and in $`zz_2=\frac{1}{\mu \mathrm{exp}\left(\beta ϵ\right)}`$, where $`\mu `$ is the connective constant of a standard SAW in $`d`$-dimensional space, i.e. $`\mu =\mu _{\stackrel{}{g}=\stackrel{}{0}}`$.

To find a lower bound for $`Z_N`$, we restrict the Y configurations to those with the joined part of the Y constrained to lay in the half-space $`\{\stackrel{}{x}|\stackrel{}{x}\stackrel{}{e}_1>0,\stackrel{}{e}_1\frac{\stackrel{}{g}}{|\stackrel{}{g}|}=0\}`$ (the origin is in the bifurcation point), and the unzipped part forced to stay in the other half-space.

$$Z_N^{lower}=\underset{m=0}{\overset{N1}{}}\mathrm{exp}\left(\beta (Nm)ϵ\right)c_{Nm1}^{left}\underset{\stackrel{}{x}_1,\stackrel{}{x}_2}{}c_{2m}^{right}(\stackrel{}{x}_1,\stackrel{}{x}_2)\mathrm{exp}\left(\beta \stackrel{}{g}(\stackrel{}{x}_2\stackrel{}{x}_1)\right),$$
(C5)
where we have indicated with $`c_n^{left(right)}`$ the number of $`n`$-step SAWs constrained to stay in the left(right) half space. As the connective constants can be rigorously proved to be the same even for walks constrained in such a way (see  for details, it is important here that walks generated from a reflection through the hyperplane $`\{\stackrel{}{x}|\stackrel{}{x}\stackrel{}{e}_1>0,\stackrel{}{e}_1\frac{\stackrel{}{g}}{|\stackrel{}{g}|}=0\}`$ leave the scalar product $`\stackrel{}{g}\stackrel{}{x}`$ unchanged), the singularities of the grand partition function obtained from (C5) are the $`\mathrm{𝑠𝑎𝑚𝑒}`$ as those of $`_{N=0}^{\mathrm{}}Z_N^{upper}z^N`$. Thus the grand partition function singularities associated to $`Z_N`$ in eq.(C3) are $`z=z_1,z_2`$ given above.

In summary, we have proved that a mechanical unzipping transition takes place when $`\mu _\stackrel{}{g}^2=\mu \mathrm{exp}\left(\beta ϵ\right)`$, similarly to the directed walk case. We suggest that a standard computation of $`\mu _\stackrel{}{g}`$ (for instance through exact enumeration) could give an accurate phase diagram for the “Y approximation” with SAWs.

The bounds in eq.(C2) give immediately two bounds within which the transition line $`g_c(T,ϵ)`$ lies. In the $`T0`$ limit eq.(33) follows since both bounds give the same asymptotic transition line. This should be the same also for generic SAWs, since in the $`T0`$ limit it is sufficient to consider only Y-shaped configurations. This demonstrates the existence of cold unzipping at sufficiently low temperature for the generic self-avoiding case.


format: text

text:
# A Search for Hi in E+A Galaxies

## 1 INTRODUCTION

In their seminal photometric study of rich clusters of galaxies at redshifts $``$ 0.5, Butcher and Oemler found a significant excess of blue, starforming galaxies in distant clusters relative to nearby clusters (Butcher & Oemler, 1978, 1984) and provided one of the clearest pieces of evidence for galaxy evolution. Although the cause for the relatively recent decrease in the mean star formation rates of cluster galaxies over the past few billion years is still not well understood, spectroscopic studies of these galaxies reveal a further puzzle. Sixty percent of these blue, supposedly starforming galaxies have strong Balmer absorption lines with little or no evidence for ongoing star formation from optical emission lines (Couch & Sharples, 1987). Because their spectra are a superposition of a young stellar population characterized by A stars and an older population characterized by K stars, typical of elliptical galaxies, these galaxies were termed E+A galaxies111We use the term E+A for historical reasons. Because the first detections of these galaxies were in distant clusters (Dressler & Gunn, 1983), “E” was used to represent the assumed morphology and “A” the dominant stellar lines of the spectrum. A better, purely spectroscopic designation would be K+A (c.f. Franx, 1993)..

E+A galaxies are believed to be post-starburst galaxies, in which star formation has been almost entirely quenched. Why star formation has ceased is unknown, although most speculate that the gas reservoir has been depleted either through highly efficient consumption or by external forces that removed the gas, such as ram pressure stripping. The purpose of this study is to detect and measure the gaseous content and morphology of E+A galaxies in an attempt to help us understand why star formation has ceased so suddenly in galaxies that were vigorously forming stars only $``$ 1 Gyr ago.

The population of E+A galaxies shows a dramatic change as a function of redshift. Dressler and Gunn (1990) found in seven clusters at 0.35 $`<z<`$ 0.45 that E+A galaxies constitute a significant fraction of the detected cluster members (23/172), while Belloni et al. (1995) found in the distant cluster Cl 0939+472 (z $`0.41`$) that 35 of the 169 member objects ($`21\pm 7`$ %) are E+A galaxies. These results have been strengthened by more recent surveys probing larger numbers of intermediate redshift clusters (e.g. Poggianti et al., 1999). On the other hand, at low redshift the fraction of E+A galaxies in clusters is close to zero. Fabricant et al. (1991) find that E+A galaxies form less than 1$`\%`$ of the population in nearby clusters.

So far, the E+A phenomenon has mostly been studied in cluster environments. The galaxy statistics of intermediate-redshift field surveys are much poorer than for clusters. The first environmentally unbiased search for E+A galaxies in the local universe was done by Zabludoff et al. (1996, hereafter Z96). From 11,113 galaxies at 0.05 $`<z<`$ 0.13 in the Las Campanas Redshift Survey (LCRS; Shectman et al., 1996), Z96 identified 21 well-defined nearby E+A galaxies, which have the strongest Balmer absorption lines (the average of the equivalent widths of $`H\beta \gamma \delta >`$ 5.5 Å ) and weakest \[O ii\] emission-line equivalent widths ($`<`$ 2.5 Å ) of all the galaxies in the survey. The fraction of nearby E+A galaxies in the LCRS sample of Z96 is less than $`0.2\%`$. However, the comparison to the higher redshift results is not straightforward because although the general spectral characteristics — strong Balmer absorption and no \[O ii\] emission — are similar for distant E+A’s and the Z96 sample, the latter sample had more stringent selection criteria. Interestingly 75$`\%`$ of the E+A’s in the Z96 sample lie in the field, well outside of clusters and rich groups of galaxies.

Although most scenarios explaining the E+A phenomenon have concentrated on star formation triggers specific to the cluster environment, the finding by Z96 that 75$`\%`$ of the local E+A sample are in the field indicates that either there are different triggers for cluster and field E+A’s or the trigger does not depend on global environment. Z96 argue that if only one mechanism is responsible for the E+A formation, then galaxy-galaxy interactions and mergers are that mechanism: optical images of at least five of the 21 galaxies in their sample show tidal features.

Finally, work by Smail et al. (1999) has cast some doubt on the post-starburst interpretation of E+A galaxies. From their radio continuum detections of E+A galaxies in a cluster at $`z=0.41`$, they suggest that these galaxies are not post-starburst galaxies, but actually host highly obscured starbursts.

Hi imaging can address a number of questions pertaining to the E+A phenomenon. If E+A galaxies are post-starburst galaxies, what happened to their gas — did the burst consume it all? A comparison of the gas content and surface density distribution between E+A’s and normal star forming galaxies might provide clues as to why E+A’s are now quiescent. If galaxy interactions are indeed responsible for triggering the original bursts, the Hi morphology and kinematics might present a tidal signature even in systems with undisturbed optical morphologies (e.g., the M81/M82/NGC3077 system; Yun, 1992). Could there be significant star formation in E+A’s that is obscured by dust (c.f. Smail et al., 1999)? Radio continuum measurements, obtained as a byproduct of the Hi observations, can set useful limits on such a scenario. Could the local environment of E+A’s play a secondary role, as a source for gas for subsequent star formation? As another byproduct of the observations, we obtain an estimate of the Hi content of neighboring galaxies and can place a limit on the total local gas reservoir for star formation.

In this paper, we present the results of a search for neutral hydrogen gas from a spectroscopically well-defined sample of nearby E+A galaxies. Five E+A galaxies, EA 1, EA 2, EA 3, EA 4, and EA 11, are drawn from the 21 E+A galaxies defined in the survey of Z96. Four of these galaxies (EA 1-4) are selected because 1) they have been observed already with the Hubble Space Telescope Wide Field/Planetary Camera 2 (Zabludoff et al., in prep.) and 2) they are easily accessible with the VLA. The fifth galaxy, EA 11, is included in our sample because, like EA 4, it lies in a known cluster environment (Z96) and thereby provides a comparison with the three field E+A’s. The numerical suffix in the galaxy name ranks the strength of the 4000 Å break, $`D_{4000}`$, in the galaxy, which is a measure of the composition of the stellar population of a galaxy. Galaxies, such as EA 1, with smaller $`D_{4000}`$ values are those most dominated by the A-stellar component and are thus bluer in color. The target galaxies lie at redshifts between 0.08 and 0.12. Our paper is organized as follows: Section 2 contains a summary of the observations and a brief description of the data reduction procedures. The Hi and the radio continuum properties of the galaxies are presented in Section 3. In Section 4, we discuss the implications of our results. We present our conclusions in Section 5. Throughout this paper we use $`H_0=100h`$ km sec-1 Mpc-1 and $`q_0=0.5`$.

## 2 OBSERVATIONS AND DATA REDUCTION

We observed the five E+A galaxies in our sample using NRAO’s Very Large Array (VLA) in January and March 1996 and June 1997. Table 1 lists the galaxy name, 1950 coordinates, heliocentric velocity ($`cz`$), B-band luminosity, luminosity distance, and the E+A environment. An E+A is identified as lying in a cluster if it is a member of a cluster in the LCRS group catalog (e.g., Tucker et al., 1998) or if there is a known cluster (from the NASA/IPAC Extragalactic Database (NED); Helou et al., 1991) at its position and redshift. Because the galaxies are at a relatively high redshift ($`z0.1`$) for an Hi detection experiment we opted for a conservative instrumental set up, using the 3 km (C) array as a compromise between optimizing surface brightness sensitivity and resolution (15 arcsec). An extended north arm was used in the first and third session to obtain a rounder beam for the sources at low elevation. We chose to probe a large velocity range (1400 km s-1), but with a rather poor velocity resolution of typically 50 km s-1; this was achieved by using a 6.25 MHz total bandwidth with 31 spectral channels and online Hanning smoothing. The observational parameters are summarized in Table 2.

We calibrated the data using standard AIPS procedures. For each field we made a continuum image covering the full primary beam centered on the E+A galaxy. We subtracted the continuum from the spectral line data by making a linear fit to the visibilities over the line-free channels, which were not affected by the edge effects of the band. For a first pass we used channels 5 to 27 to make the fit; if Hi was detected, a new fit was made excluding the channels with line emission. The search for Hi was done in the UV data and in the image plane. In the UV plane, we first vector-averaged over all baselines and over time, which is equivalent to making a natural-weighted synthesized beam at the field center. We examined the resulting spectrum for Hi line emission. By varying the weighting of data, we also obtain different spatial resolutions and, for example, can directly search for extended emission at the field center. In the image plane, we examined cubes by eye before and after continuum subtraction. We also smoothed the cubes spatially and in velocity to search for emission that was either spatially extended or covered many channels. All the cubes were made by using natural weighting to optimize sensitivity. The rms noise in the full resolution cubes is typically 0.4 mJy beam-1, which is equivalent to a few $`\times `$ 1019 cm-2 in column density sensitivity.

## 3 RESULTS

Of the five galaxies, only one was detected in emission — EA 1. No significant Hi emission was detected from any of the other galaxies, nor was radio continuum detected from any of the five galaxies. In the following sections, we present the Hi data for EA 1 along with the detection limits for the other systems without Hi detections. The general results are summarized in Table 3.

### 3.1 EA1

We detect a total Hi flux from EA 1 of 0.30 $`\pm `$ 0.02 Jy km s-1. We calculate that the mean column density is 1.3 $`\times `$ $`10^{20}`$ cm-2 and that the total Hi mass is 3.5 $`\pm `$ 0.2 $`\times `$ $`10^9`$ $`h^2`$ M. In the channel maps, shown in Fig. 1, we see emission associated with the galaxy over approximately 150 km s-1, from 22284 to 22428 km s-1. There is a north-south velocity gradient in the system with the north side approaching. Near the central velocity the emission broadens, consistent with the emission from a moderately inclined, asymmetric disk. Considering the very disturbed optical morphology (discussed below) it is also possible that the Hi arises in two kinematically distinct components that overlap in the central regions. A profile of the central Hi emission, shown in Fig. 2, also shows a slight asymmetry. At the extreme velocities Hi appears northwest and southeast of the galaxy at projected distances of about 50 $`h^1`$ kpc. Although this emission is weak, it almost certainly is real — it joins smoothly into the main body of emission in velocity structure and, at very low levels, in spatial distribution as well. Morphologically and kinematically, the extended emission appears very similar to the gas-rich tidal tails observed in nearby mergers (e.g., Hibbard & van Gorkom, 1996).

In Fig. 3, we show the Hi emission integrated over all velocities overlaid on an HST WFPC2 image (Zabludoff et al, in prep.). The overlay further supports the idea that the extended Hi is tidal debris, and the HST image shows that EA 1 looks very disturbed, much like an ongoing merger. The HST image also shows another somewhat disturbed galaxy at a small projected distance; unfortunately, we do not have a redshift for this galaxy. The only hint that it may be a true companion is that the Hi is extended in the direction toward the neighboring galaxy. In Fig. 4 we show a position-velocity profile along the kinematical major axis of the Hi, approximately the line connecting the tips of the tidal tails. The asymmetry in the central Hi distribution toward higher velocities is obvious. The extended Hi connects in velocity smoothly with the central gas, as would be expected from gas that has been tidally removed. Finally in Fig. 5 we show the intensity weighted velocity field of the Hi associated with EA 1.

EA1 had already been noted as having tidal features by Z96; however, it is the combination of the HST image and the VLA Hi results that show how truly strange this system is. The HST image suggests this is a very young, perhaps ongoing, merger. The two Hi tails suggests that two gas-rich disk galaxies are involved, yet the current star formation rate as derived from the \[O ii\] is remarkably low. From the HST image it is likely that the system is dusty, hence we may underestimate the amount of \[O ii\] emission due to extinction. However, a limit to the radio continuum emission, discussed in §4, confirms the very low upper limit on the current star formation rate. We can use the velocities of the Hi tails to give a rough estimate of the time elapsed since the interaction began (Schweizer, 1982). Taking $`t=0.8r_pv^1\mathrm{cot}\beta `$, where $`r_p`$ is the projected radius from the center to the tip of the tail, $`v`$ is the relative LOS velocity of that tip, and $`\beta `$ is the angle between the radius vector and the line of sight. The formula assumes that the tail motions are mostly radial, and the factor 0.8 corrects for minor decelerations and the non-central tail origin. Taking $`\beta =45^{}`$, $`r_p=52h^1`$ kpc, and $`v=70`$ km s-1, we get $`t=6\times 10^8h^1`$ yrs. While it is unclear if the Hi tails are related to the event that caused the starburst, the derived time scale is consistent with the notion that an interaction in the last Gyr triggered the original starburst that gave rise to the E+A spectrum.

### 3.2 The Hi Detection Limits

The smallest Hi mass that can be detected anywhere in a datacube is a 5$`\sigma `$ signal per beam per channel. In Table 4 we list the corresponding Hi mass limits, using M$`{}_{HI}{}^{}=2.35\times 10^5\times d_L^2\times \frac{{\scriptscriptstyle \mathrm{𝑆𝑑𝑉}}}{1+z}h^2`$ M, where $`d_L`$ is the luminosity distance in $`h^1`$ Mpc, $`\mathrm{𝑆𝑑𝑉}`$ is the integrated Hi emissivity in Jy km s-1, and $`z`$ is the redshift of the galaxy. These values are valid at the field center. At larger distances they have to be corrected for the primary beam response because, for example, at a distance of 15 from the field center, these limits are a factor two larger. Our limits are for detections at random positions in the cubes.

A signal at the position and velocity of an E+A galaxy would be believable at the 3$`\sigma `$ level. Because our only detection, EA1, is the most nearby of the E+A galaxies, we have also calculated whether Hi emission with the spatial and spectral characteristics of the gas in EA1 would be detectable at the distances of the other E+A galaxies. In the last two columns of Table 4 we list the rms noise after smoothing to the velocity extent of EA 1 and the expected signal per beam at the location of the other galaxies. EA 1 would have been detectable at the 3$`\sigma `$ level even at the greater distances of the other E+A galaxies with the possible exception of EA 11, for which the signal would only be 2.5$`\sigma `$.

### 3.3 The Continuum Detection Limits

In Table 5 we list the rms noise in the continuum images and the 3$`\sigma `$ detections limits at 1.3 GHz at the location of the E+A galaxies. The power per synthesized beam is $`P_{1.3}=1.2`$ $`\times `$ $`10^{17}`$ $`\times `$ $`S`$ $`\times `$ $`d_L^2`$, where $`P_{1.3}`$ is the power of the source at 1.3 GHz in $`h^2`$ W Hz-1, $`S`$ is the flux in mJy, and $`d_L`$ is the luminosity distance in $`h^1`$ Mpc. None of the five E+A’s are detected in continuum emission.

## 4 DISCUSSION

### 4.1 Hi Properties

Of our five spectroscopically selected E+A galaxies, we have detected Hi in one system, with an Hi mass of 3.5 $`\times `$ 109 $`h^2`$ M, while the other four non-detected systems have Hi upper limits of order 109 $`h^2`$ M. What does this small sample tell us about E+A galaxies as a class? The first important lesson is that not all E+A galaxies are very gas deficient — it is possible to stop star formation without exhausting all the gas. Secondly, even within our rigorously defined set of E+A galaxies there is a large spread in Hi properties. Although our one detected galaxy, EA 1, is the bluest galaxy in the entire sample of Z96, it is not unique in any sense. EA 1’s only distinction is its extremely disturbed optical morphology — two components are visible — but two galaxies with Hi non-detections, EA 2 and EA 3, are just slightly redder and are also tidally disturbed. Z96 note that even with their strict selection criteria, the spectra of the 21 E+A galaxies span a range of spectral characteristics that suggest that these objects are at different evolutionary stages, have experienced different evolutionary histories, and/or had morphologically different progenitors. Can we shed more light on the evolutionary status of these objects from their Hi properties?

The Z96 sample of E+A galaxies is far more extreme than any other E+A sample in terms of its spectroscopic properties, simply because the large number of galaxies ($``$11000) in the LCRS allows us to find extremely rare objects. Thus we can only compare the Hi properties of our sample with the little Hi data that exist for previously known, but less well-defined, E+A galaxies in the local universe, such as the merger remnants of Liu & Kennicutt (1995a, b) and the cluster galaxies of Caldwell and Rose (1997). The qualitative description of the E+A spectra, in general, is the strong Balmer absorption lines and the very weak emission from \[O ii\] and \[O iii\] forbidden lines. Merger remnants have strong Balmer absorption lines, but also significant emission from \[O ii\] lines (Liu & Kennicutt, 1995a, b; Schweizer, 1996; Kennicutt, 1992; Carter et al., 1988; Bica & Alloin, 1987). Several of these galaxies have been detected in Hi emission. For instance, two merger remnants, NGC 7252, which has clear tidal tails with a total Hi mass of 2 $`\times `$ $`10^9`$ $`h^2`$ M (Hibbard et al., 1994), and NGC 520, with a total Hi mass of 3.9 $`\times `$ $`10^9`$ $`h^2`$ M, exhibit strong Balmer absorption lines but have 6 Å \[Oii\] emission line equivalent widths (Liu & Kennicutt, 1995a). These galaxies have been called E+A galaxies in the literature, but clearly their spectroscopic properties are different from those of our E+A sample: contrary to our E+As these galaxies have low level ongoing star formation. Moreover, optically EA1 looks like a much younger object than, for example, the merger remnant NGC 7252, and more similar to the Antennae or NGC 6240, both of which are vigorously forming stars. This begs the question, why is EA 1’s current star formation rate so low? Could it be that the gas in EA 1, although significant in Hi mass, is too rarefied to form stars? Unfortunately, the relatively low spatial resolution of our Hi data makes it impossible to resolve the inner Hi mass distribution in EA 1 and address these questions. Furthermore, the upper limits to the Hi mass in EA 2 and EA 3 are not low enough to coin these systems as Hi deficient for merging or interacting systems. That EA 1 is bluer and apparently more gas-rich than EA 2 and EA 3 could point to a difference in gas properties of the progenitors of these systems or a different evolutionary history. On the other hand, EA 1 could represent an early stage of the same evolutionary sequence, before the gas disperses to the point where it is undetectable in EA 2 and EA 3.

We have measured only upper Hi mass limits for the two E+A galaxies, EA 4 and EA 11, in a cluster environment. A recent Hi study of 11 of the 15 E+A galaxies identified in the Coma cluster by Caldwell et al. (1993) and Caldwell and Rose (1997) placed upper limits to the Hi mass of these galaxies of a few $`\times `$ 107 $`h^2`$ M (Bravo-Alfaro et al., 2000). Because our Hi mass limits are larger and because the criteria for E+A galaxies differ between the samples, a direct comparison is not possible. Observations of larger, well defined, samples are necessary to determine whether cluster E+A galaxies have less Hi than their field counterparts.

The Hi upper limits also do not constrain the galaxy type. There is a large spread in Hi content for galaxies of all morphological types (Roberts & Haynes, 1994), and essentially no type can be ruled out based on our upper limits. Although the Hi content of individual galaxies varies enormously for a given Hubble type, the ratio of Hi mass and blue luminosity is more tightly constrained, log$`(M_{HI}/L_B)`$ varies smoothly from $`1`$ to just under $`0`$ from early (SOa) to late (Sm, Im) types (Roberts & Haynes, 1994). In Table 4 we list this ratio for the five galaxies. The blue luminosities are calculated from HST PC2 images in the F439W filter (Zabludoff et al., in prep.). For each galaxy, we determine the total magnitude within a 100 pixel radius, correcting for Galactic extinction and applying a $`k`$-correction based on the spectral energy distribution of an A star. The F439W band is approximately the Johnson B-band for an A stellar spectrum. EA 1 stands out with a rather high log$`(M_{HI}/L_B)`$ value of $`0.15`$; interestingly, this value is similar to that for tidal tails in gas rich mergers (Hibbard et al., 1994). The upper limits for the other E+A galaxies fall within the range of values for disk galaxies.

We also search a volume of $`40\mathrm{}\times 40\mathrm{}\times 1400`$ km s-1 around each of the galaxies. Taking into account the primary beam response, the total volume examined down to $`5\times 10^9h^2`$ M is $`340h^3`$ Mpc3. Using the Hi mass function derived by Zwaan et al. (1997) we expect to detect three galaxies with $`5\times 10^9`$ $`h^2`$ M of Hi in such a volume. We only detect one system. This result is somewhat lower than the expected three, but the difference is perhaps not statistically significant. However, we note that our fields are not randomly chosen, but are pointed at known galaxies, and galaxies tend to cluster. Early type galaxies detected in Hi are preferentially found in Hi rich groups (van Gorkom and Schiminovich, 1997), many late type (gas-rich) galaxies live in groups, and even cluster gas-rich galaxies tend to be subclustered (e.g., van Gorkom, 1996). The detection of only one large Hi mass in the entire volume suggests that the environment of the E+A galaxies is mildly Hi deficient.

### 4.2 The inferred star formation rate

The claim that the E+A galaxies have no ongoing star formation is based on the lack of significant optical emission lines. Our detection of Hi in EA 1 and our upper limits on Hi and radio continuum emission provide an independent way to constrain the on-going star formation rate. In this section we compare the rates estimated using the measured Hi column density, radio continuum emission, and \[O ii\] equivalent widths of the E+A galaxies. The results are listed in Table 6.

#### 4.2.1 Hi Limits

Although stars form out of molecular gas, Hi appears to trace the star formation rate (SFR) remarkably well (Kennicutt, 1989). Is the detection of a few $`\times `$ 109 $`h^2`$ M of Hi in EA 1 in contradiction with the claim that no on-going star formation occurs in E+A galaxies? We can make a crude estimate of the SFR in EA 1 using the Schmidt law, which relates the gas surface density to the star formation rate by a simple power law. According to Kennicutt (1998),

$`\mathrm{\Sigma }_{SFR}`$ = (2.5 $`\pm `$ 0.7) $`\times `$ $`10^4`$ $`\mathrm{\Sigma }_{gas}^{1.4\pm 0.15}`$ M yr-1 kpc-2,

where $`\mathrm{\Sigma }_{gas}`$ is the total (atomic + molecular) gas surface density averaged over the disks in M pc-2. The mean Hi surface density in the core of EA 1 is 2.2 $`\times `$ $`10^{20}`$ cm-2, or $``$ 1.8 M pc-2. If we simply adopt the Schmidt law to estimate the SFR, taking the mean Hi surface density of the central region of EA 1 as the total gas surface density and estimating the size of this region to be $``$ 940 $`h^2`$ kpc2, we get a SFR of $``$ 0.5 $`h^2`$ M yr-1. The upper limits on the Hi mass in the other E+A galaxies provide upper limits on their SFR of the same order of magnitude (Table 6).

Skillman (1986) points out, however, that there is a surface density threshold below which no star formation occurs. He finds that the observed local Hi column density threshold for star formation is 1 $`\times `$ $`10^{21}`$ cm-2, at a resolution of 500 pc. Although our resolution is of order 10 kpc, it is still interesting to note that the Hi column density of EA 1 and the column density upper limits for the other E+A galaxies are well below this threshold. If the gas is distributed smoothly over the beam area, then we might expect no on-going star formation in any of the E+A galaxies even though ample Hi exists.

#### 4.2.2 Radio Continuum Limits

Radio continuum emission associated with star-forming galaxies is thought to be due to synchrotron radiation produced by relativistic electrons in supernovae remnants. Therefore, the radio continuum emission provides a probe of massive (M $`>`$ 8 M) star formation in galaxies, independent of dust obscuration. Condon (1992) finds a relationship between 1.4 GHz luminosity and massive SFR of

SFR1.4(M $`>`$ 5 M) = $`\frac{L_{1.4}}{4.0\times 10^{21}WHz^1}`$ M $`yr^1`$.

assuming a Miller-Scalo IMF and using the Galactic radio luminosity and supernova rate. We list upper limits to the star formation rate using this formula in Table 6.

Smail et al. (1999) have cast some doubt on the post-starburst interpretation of the E+A phenomenon after detecting a significant fraction of cluster E+A galaxies at z=0.41 in radio continuum. Radio continuum emission can either be powered by active galactic nuclei (AGN) or by massive star formation. Observationally the two can easily be distinguished by morphology (jet/core morphology vs. extended disk emission). The radio power resulting from the two mechanisms is also quite different with a small range of overlap. The radio continuum luminosity from AGN is typically L $``$ $`10^{2326}`$ $`h^2`$ W Hz-1 at 1.4 GHz, while that due to star formation is L $``$ $`10^{2123}`$ $`h^2`$ W Hz-1 (Condon, 1989). The (unresolved) detections by Smail et al. (1999) in CL 0939+4713 at z$`=0.41`$ have radio luminosities of L $``$ $`10^{2223}`$ $`h^2`$ W Hz-1 and, based on their radio luminosity alone, could either be dust-enshrouded starbursts or AGN. However, these authors suggest that the fact that 4 out of 5 of those systems are disk galaxies favors the starburst interpretation. Our upper limits of a few $`\times `$ $`10^{21}`$ $`h^2`$ W Hz-1 make it unlikely that our sample of E+A galaxies host radio loud AGN. Although we cannot rule out ongoing star formation at a low level (Table 6), our upper limits are an order of magnitude smaller than the detections by Smail et al. (1999). We do not find evidence for hidden starbursts in the most extreme local E+A’s.

#### 4.2.3 $`[`$O ii$`]`$ Limits

The \[O ii\] emission line strength is another indicator of on-going star formation. The derived star formation rates from \[O ii\] are less accurate than those derived from H$`\alpha `$, but are comparable in accuracy to rates derived from the optical continuum flux and colors (Kennicutt, 1992b). Barbaro & Poggianti (1997) find that if the number of ionizing stars of an H ii region is small ($``$ 50), the \[O ii\] equivalent width itself is a good SFR indicator; however, the method fails if the number of ionizing stars is much larger than 50 (Stasinska & Leitherer, 1996). Because of their low current SFR’s, the E+A’s are probably in a low ionization regime and the \[O ii\] equivalent width should provide a rough estimate of the SFR in the E+A galaxies.

The relation between the \[O ii\] luminosity and equivalent width is given by (Kennicutt, 1992b)

L$`{}_{[\mathrm{O}\mathrm{II}]}{}^{}1.4\times 10^{29}\times L_{B,}\times \mathrm{EW}[\mathrm{O}\mathrm{II}]`$ erg s-1,

where $`L_{B,}`$ is the B band luminosity in $`L_{}`$, and EW stands for equivalent width. The $`L_{B,}`$ values were derived from our HST data and are listed in Table 1. Using the relationship between SFR and L$`_{[\mathrm{O}\mathrm{II}]}`$ (Barbaro & Poggianti, 1997):

SFR(M yr-1) $``$ 6.3 $`\times `$ $`10^{41}`$ L$`_{[\mathrm{O}\mathrm{II}]}`$

we calculate the star formation rates (Table 6) for all E+A galaxies. A 2$`\sigma `$ upper limit was used for the line intensity in the E+A galaxies that were not detected in \[O ii\].

The different methods used to derive limits on the current SFR in E+A galaxies each have their own limitations. The estimates based on the \[O ii\] equivalent width could be affected by extinction. Although there is a tight relation between radio continuum and other indicators of SFR (specifically CO emission and FIR flux), the reason for this correlation is not well understood. The SFR associated with Hi must also depend on the volume density of the gas, which is not known. Nonetheless, it is reassuring that the low star formation rates indicated by the \[O ii\] equivalent widths are not that much lower than the other values. These E+A galaxies are forming stars at a very modest rate, if at all.

## 5 CONCLUSIONS

This first search for gas in a rigorously defined set of the most extreme local E+A galaxies finds $`3.5\times 10^9h^2`$M of Hi in one galaxy and obtains upper limits for four more. Thus not all E+A galaxies are gas deficient; E+A formation does not always exhaust the entire gas reservoir.

In the one system detected in Hi, the morphology of the gas supports the hypothesis that the original starburst was induced via a galaxy interaction. The time elapsed since the interaction started, derived from the Hi kinematics, is roughly consistent with the time since the starburst must have stopped. However, the fact that the system contains a substantial amount of Hi, and that the merger appears not yet complete, raises questions of why the starburst has ceased.

None of the E+A galaxies in our sample were detected in radio continuum. The corresponding upper limits on the star formation rate rule out that these galaxies are highly obscured, extremely vigorous star forming systems (i.e. $`>`$ 10 M yr-1). Using spectral synthesis models (Bruzual & Charlot, 1993, and later models), one can show that the starbursts that produce E+A galaxies formed $`10\%`$ of their stellar mass during a burst of duration $`10^9`$ yrs; therefore, the corresponding SFR during the burst must be $``$ 10 M yr-1. From a comparison of this SFR and our observational limits, we conclude that the SFR in these systems must have dropped by at least an order of magnitude within the past 1 Gyr.

Although we can exclude large SFR’s, none of the available star formation indicators exclude star formation rates at a level comparable to or below that of a weak spiral galaxy. The \[O ii\] emission or lack thereof implies even lower star formation rates than inferred from the Hi or radio continuum luminosities, but the \[O ii\] estimate could be affected by dust.

We thank Tod Lauer for providing the B-band luminosities of the sample galaxies. This research was supported by an NSF grant to Columbia University. AIZ acknowledges support from NASA grant HF-01087.01-96A and HST grant GO-06835.01-95A. DZ acknowledges support from the David and Lucile Packard Foundation the Alfred P. Sloan Foundation, and a NASA LTSA grant (NAG 5-3501). JCM acknowledges support from an NSF CAREER award (NSF AST-9876143). The National Radio Astronomy Observatory is operated by Associated Universities, Inc., under cooperative agreement with the National Science Foundation.