text:
# Phase-field study for the splitting mechanism of coherent misfitting precipitates in anisotropic elastic media

## Abstract

Coherent misfitting precipitates in elastically stressed media such as $`\gamma `$ particles in nickel-based super-alloys show various splitting patterns such as doublets, quartets, or octets due to their misfit strain energy. While it is an interesting instability phenomenon defying conventional surface thermodynamics, its mechanism is not completely clear. Through a phase-field study upon the splitting behavior and morphological evolution of coherent precipitates, we show that an interface instability driven by elastic anisotropy and a diffusion field can generate elastically induced splitting during diffusional phase transition. Particle splitting is triggered by interface grooving which advances by penetrating grooves into the interior of the particle. The sequential evolution of shapes during the splitting process is in good agreement with previous experiments.

Elastic stresses arising during solid state phase transitions in some alloys can lead to microstructural features such as precipitate alignment, cuboidal precipitate shapes, and particle splitting. A wealth of experimental and theoretical works have reported how the mismatch strain between the precipitate and the surrounding matrix phase changes the characteristics of the microstructure Ardell66 ; Miyazaki79 ; Yoo95 ; Fratzl99 ; Miyazaki82 ; Khachaturyan88 ; Wang93 ; Zhang98 ; Lee ; Doi92 .

One of the primary reasons for considering the role of elasticity in microstructure formation is that it offers us the possibility of using elastic stresses to design desirable microstructures. In this context, one of the possibilities often cited is that of inverse coarsening. In inverse coarsening, the elastic fields promote the growth of a small precipitate at the expense of a larger neighbor - exactly the opposite to surface energy driven coarsening. In principle, elasticity could favor an array of equally sized precipitates, which would lead to the improved mechanical properties of the alloy.

Particle splitting is perhaps the most closely linked elastic effect to inverse coarsening. Particle splitting is a process in which one cuboidal particle splits into several particles, usually two (a doublet) or eight (an octet). Such splitting has been observed experimentally Miyazaki79 ; Yoo95  and has been theoretically identified as an elastically driven process Miyazaki82 ; Khachaturyan88 .

While splitting may be an important factor for inverse coarsening, its mechanism and the precise experimental and theoretical conditions under which it occurs are not completely clear. Some experimental pictures suggest that splitting occurs via a morphological instability in which the sides of a cuboidal precipitate become increasingly concave prior to splitting Yoo95 . Other pictures suggest a hollowing at the center of the particle Miyazaki82 . Many theoretical simulations of splitting have used diffuse interface techniques, and suggest that the instability occurs by the ’hollowing’ mechanism Wang93 ; Zhang98 . However, the hollowing mechanism cannot explain the formation of an octet and the fact that the splitting occurs mostly in the concave particle shape Miyazaki79 ; Yoo95  though it can describe splitting to a doublet. Recently, Lee, using a discrete atom method, observed splitting by an instability mechanism Lee .

In this letter we demonstrate, through phase field model calculations, that interface instability driven by elastic anisotropy and diffusion field can generate particle splitting. Splitting occurred only under the concave growth condition and the picture that a concentrated elastic field around interface groove promotes splitting is very similar to Asaro-Tiller-Grinfeld (ATG) instability atg . In addition, this splitting mechanism can explain the formation of both doublets and octets.

The model in this study assumes coherent interfaces (i.e. there are no dislocations) and takes into account the effect of anisotropic elasticity and elastic inhomogeneity which represents that the elastic constants of the matrix and the precipitate differ. We assume that the misfit between the matrix and the precipitate is purely dilatational with magnitude $`e^T`$, the elastic constants of both phases have cubic (four-fold) symmetry, and that the lattice parameter is not a function of composition. The model is formulated in terms of the concentration field $`c(\stackrel{}{r},t)`$, the phase field $`\varphi (\stackrel{}{r},t)`$ which is zero in the matrix phase and one in the precipitate, and the elastic strain tensor $`u_{ij}=1/2\left(u_i/x_j+u_j/x_i\right)`$, where $`\stackrel{}{u}(\stackrel{}{r},t)`$ is the atomic displacement field. The interfacial region is defined to be a mixture of matrix and precipitate with a different composition but with the same chemical potential skkim . Therefore, the concentration field of the system is defined as $`c=c^Ph(\varphi )+c^M[1h(\varphi )]`$, where $`c^A`$ is a concentration in the A-phase and $`h(\varphi )`$ is the smooth interpolation function, which is zero at $`\varphi =0`$ and one at $`\varphi =1`$ and changes smoothly from zero to one, i.e. $`0<\varphi <1`$. The superscript $`P`$ and $`M`$ represent the precipitate and matrix phases, respectively, which will be used later without description. From the definition of the interfacial region, the following constraint should be satisfied at any point in the system:

$$f_{c^P}^P(c^P)=f_{c^M}^M(c^M)$$
(1)
where $`f^A`$ is the chemical free energy of the A-phase and $`f_{c^A}^A`$ represents the derivative of $`f^A`$ with respect to $`c^A`$. The free energy functional is:

$$=_\stackrel{}{r}\left[f(\varphi ,c,u_{ij})+\frac{\epsilon ^2}{2}|\varphi |^2\right]𝑑\stackrel{}{r}$$
(2)
where $`f(\varphi ,c,u_{ij})`$ is the free energy density of the system including the coherent strain energy, defined as $`wg(\varphi )+f^P(c^P)h(\varphi )+f^M(c^M)[1h(\varphi )]+f^{el}(\varphi ,u_{ij})`$, where $`g(\varphi )`$ is the double well potential, which has minima at 1 and 0, the equilibrium values of the precipitates and matrix phases,respectively. In this study, $`g(\varphi )`$ was selected as $`\varphi ^2(1\varphi )^2`$. $`f^{el}`$ is the coherent elastic energy density caused by the lattice mismatch between two phases. In general, the elastic energy density is expressed as Leo :

$$f^{el}=\frac{1}{2}C_{ijkl}(\varphi )[u_{ij}e^T\delta _{ij}\psi (\varphi )][u_{kl}e^T\delta _{kl}\psi (\varphi )]$$
(3)
where $`C_{ijkl}`$ is the stiffness tensor that depends on the phase field due to elastic inhomogeneity and $`\psi (\varphi )`$ is a smooth function with the same characteristics as $`h(\varphi )`$. A summation convention over repeated indices is implicit. Considering the spatial inhomogeneity of the elastic stiffness, the stiffness tensor is given by:

$$C_{ijkl}(\varphi )=C_{ijkl}^0+C_{ijkl}^{}\xi (\varphi )$$
(4)
where $`\xi (\varphi )`$ is a smooth function which is zero at $`\varphi =1`$ and one at $`\varphi =0`$, $`C_{ijkl}^0=C_{ijkl}^P`$, and $`C_{ijkl}^{}=C_{ijkl}^MC_{ijkl}^P`$, where $`C_{ijkl}^A`$ is the elastic stiffness of the A-phase.

It is reasonable to suppose that the elastic field relaxes much faster than $`c`$ or $`\varphi `$. The elastic field can then be solved in terms of the phase field using the condition of local mechanical equilibrium:

$$\frac{\delta }{\delta u_i}=_j\sigma _{ij}=0$$
(5)
where $`\sigma _{ij}=\delta /\delta u_{ij}`$ is the stress tensor:

$$\sigma _{ij}=C_{ijkl}(\varphi )[u_{kl}e^T\delta _{kl}\psi (\varphi )]$$
(6)
In this work, in order to obtain accurate solutions of equation (5) with minimal calculation time, the high order approximation method similar to the method proposed by Chen and coworkers Chen  is used instead of the iteration method which gives very accurate solutions but requires considerable calculation time.

The solution of equation (5) to zeroth order in the elastic stiffness $`C_{ijkl}^{}`$ is

$$\widehat{u}_k^0=i\alpha _{ki}^1\sigma _{ij}^0q_j\widehat{\psi }(\widehat{\varphi })$$
(7)
where $`\widehat{A}`$ represents the value of $`A`$ in Fourier space, $`q_j`$ the *j*-direction component of the wave number vector, $`\sigma _{ij}^0=C_{ijkl}^0e^T\delta _{kl}`$, and $`\alpha _{ki}^1`$, indicates a component of the inverse tensor of $`\alpha _{ik}`$ which is defined as $`C_{ijkl}^0q_jq_l`$. Then, the displacement field to the *n*th order becomes

$`\widehat{u}_k^n=\widehat{u}_k^0`$ $`+`$ $`i\alpha _{ki}^1C_{ijlm}^{}e^T\delta _{lm}q_j\{\xi (\varphi )\psi (\varphi )\}_F`$ (8)
$``$ $`i\alpha _{ki}^1C_{ijlm}^{}q_j\{[\xi (\varphi )_m]u_l^{n1}\}_F`$
where $`u_k^{n1}`$ is the displacement field to the *n-1*th order and $`\{\}_F`$ represent the Fourier transform of the term in braces.

The elastic field can now be expressed in terms of the phase field $`\varphi `$. Substituting the solution for the strain field gives the free energy in terms of $`c`$ and $`\varphi `$. Assuming that the system evolves in time so that its total free energy decreases monotonically, the evolution equation for the phase field becomes:

$$\frac{\varphi }{t}=M\frac{\delta }{\delta \varphi }=M\left(\epsilon ^2^2\varphi f_\varphi \right),$$
(9)
with

$$f_\varphi =h^{}(\varphi )\{f^Pf^M+(c^Mc^P)f_{c^M}^M\}+f_\varphi ^{el}+wg^{}(\varphi )$$
(10)
while the concentration field should evolve with:

$`{\displaystyle \frac{c_k}{t}}`$ $`=`$ $`{\displaystyle \frac{D(\varphi )}{f_{cc}}}{\displaystyle \frac{\delta }{\delta c}}`$ (11)
$`=`$ $`D(\varphi )c+{\displaystyle \frac{D(\varphi )}{f_{cc}}}f_{c\varphi }\varphi `$
where $`M`$ is the phase-field mobility, $`D(\varphi )`$ is the diffusivity of solute atoms, and $`f_\varphi `$, $`f_\varphi ^{el}`$, $`f_{cc}`$, and $`f_{c\varphi }`$ are $`f/\varphi `$, $`f^{el}/\varphi `$, $`^2f/c^2`$, and $`^2f/c\varphi `$, respectively.

In the sharp interface limit, the phase field model presented here reduces to a modified Gibbs-Thompson equation with a coherent interface as follows:

$$\alpha \sigma v_n=\kappa _c\sigma +\mathrm{\Delta }\mu _C+G^{coh}+G^{trans}$$
(12)
Here, $`\alpha =D/M\epsilon ^2`$, $`v_n`$ is the interface velocity, $`\sigma `$ the interface energy, $`\kappa _c`$ the curvature of the interface, and $`\mathrm{\Delta }\mu _C`$ the chemical driving force of the phase transition defined as $`\mu _C^P\mu _C^M`$, where $`\mu _C^A`$ represents the chemical potential at the A-phase side of the interface. $`G^{coh}`$ means the elastic energy required to maintain coherence at the interface, defined by $`(u_{ij}^Mu_{ij}^P)\sigma _{ij}^M`$, where $`u_{ij}^A`$ and $`\sigma _{ij}^A`$ are elastic strain and stress tensors at the A-phase side of the interface, respectively LeoGT . $`G^{trans}`$ represents the difference in the elastic strain energies of the precipitate and the matrix, as defined by $`1/2\sigma _{ij}^P(u_{ij}^Pe^T\delta _{ij})1/2\sigma _{ij}^Mu_{ij}^M`$ LeoGT . The detailed derivation of Eq. (12) will emerge elsewhere further .

Numerical simulations on a discrete lattice were performed in two-dimensions. Euler’s method was used for the integration in time. For all simulations presented here, the mesh size $`\mathrm{}x=1.0\times 10^8m`$, the time step $`\mathrm{}t=0.0125`$, and $`\epsilon =3.52\times 10^5`$. This choice of $`\mathrm{}x`$ and $`\epsilon `$ guarantees that the interface is resolved by at least five points and the interfacial energy is 0.0243 $`J/m^2`$. We assume that the diffusivities of two phases are the same and $`1\times 10^{15}m^2/s`$. Quadratic functions were used for the free energy functions of both phases for simplicity, that is, $`f^P=\mathrm{\Delta }^P(c^Pc_e^P)^2/V_m`$ and $`f^M=\mathrm{\Delta }^M(c^Mc_e^M)^2/V_m`$, where $`\mathrm{\Delta }^P`$ and $`\mathrm{\Delta }^M`$ were 10000 and 2000, respectively, $`V_m`$ is molar volume and was set at $`7.0\times 10^6m^3/mol`$, and $`c_e^P`$ and $`c_e^M`$ are the equilibrium concentrations without elastic field and were set at 0.2 and 0.4, respectively. The elastic constants of two phases (in units of $`10^{11}erg/cm^3`$) are $`C_{11}^P=1.67`$, $`C_{12}^P=1.07`$, $`C_{44}^P=0.99`$, $`C_{11}^M=1.12`$, $`C_{12}^M=0.63`$, and $`C_{44}^M=0.57`$, respectively; for the Ni-based superalloy Miyazaki82 . The misfit strain between the precipitate and the matrix is dilatational and assumed to be $`0.8\%`$. The displacement field was calculated to the third order from Eq. (8) because the inhomogeneities of the elastic constants were very large. Periodic boundary conditions were employed in all directions.

Figure 1 shows the calculated splitting behavior of the precipitate and the corresponding elastic energy density distribution at various times. The morphological instability of the interface due to the diffusion field and anisotropy of the elastic stiffness is evident and this interface instability leads to particle splitting. Despite the simplicity of the two-dimensional calculation, the splitting pattern is in good agreement with the experimental results (see Fig. 1i and j). The scenario of the splitting is as follows: initially the precipitate grows in a circular shape due to the isotropic interface energy, and then transforms to cuboidal shape above a critical size, because the misfit strain energy proportional to its volume dominates the interfacial energy in proportion to the interfacial area. Subsequently, the growth rate in the $`<11>`$ direction is enhanced due to a high chemical driving force while the migration rate in the $`<10>`$ direction is suppressed due to enriched solute and a corresponding low chemical driving force. Therefore, the growing shape transforms to a concave shape. In the concave growth regime, the interface in the $`<10>`$ direction has a negative curvature, and so, the elastic strain energy concentrates around the concave edge . This concentrated elastic energy drives the dissolution of the precipitate, and then a groove forms and penetrates the precipitate (see Fig. 1c). This situation resembles Asaro-Tiller-Grinfeld instability, in which elastic media under the nonhydrostatic stress show the morphological instability atg . Through various numerical simulations, we confirmed that the splitting occurs only in the concave growth condition and that it is induced by the interface instability. This is in good agreement with the experimental observations Ardell66 ; Miyazaki79 ; Yoo95 .

Figure 2 shows the sequential variations of the $`<10>`$ length ($`a_1`$) and $`<11>`$ length ($`a_2`$) of the precipitate and the ratio of the two lengths which is called a shape parameter ($`\zeta `$) hereafter. There are four different growth modes (see Fig. 2 b). In the first growth mode the growth rates of $`a_1`$ and $`a_2`$ are proportional to $`t^{0.35\pm 0.005}`$ and $`t^{0.55\pm 0.01}`$, respectively and hence transition from circular through cuboidal to concave shapes occurs (for a perfect square $`\zeta =0.7`$). The second regime can be called the shape-conserving mode because the shape parameter remains constant in this period. The growth rates of $`a_1`$ and $`a_2`$ have the same dependence $`t^{0.31\pm 0.01}`$ on time. The growth of $`a_1`$ stops upon entering the third growth mode while $`a_2`$ continues to grow and the shape parameter starts to decrease again. In this regime the grooves are generated on the four interfaces of the $`<10>`$ direction (see Fig. 1c) and concavity around the grooves starts to relax because the elastic field around the grooves is relieved by local concentration of the stress field on the edges of the grooves. One can see this by comparing Fig. 1b with Fig. 1c. Once the elastic energy concentrated on the tips dominates the chemical driving force for precipitate growth which is reduced due to the buildup of solute around the grooves, the grooves penetrate into the precipitate and the fourth growth stage begins. Note that the particle in the non-splitting case also shows concave growth, but then returns to the cuboidal shape with time (see the variations of the filled triangles in Fig. 2b).

Fig. 3 shows the variations of the total free energy \[Eq. (2)\], elastic energy, and the total free energy excluding the elastic energy, with the evolution of the second particle. Here, the third free energy term mentioned above is called the chemical free energy although it includes the interfacial energy (i.e. gradient energy and double well potential energy). During the evolution of the particle, the total free energy and the chemical free energy decrease while the strain energy increases due to misfitting particle growth. When the splitting has been almost completed, the chemical free energy jumps due to the increased interfacial area, but the total strain energy reduces by interface grooving, which actually resembles the ATG instability atg . As the decrease of strain energy overwhelms the increased interfacial free energy, the total free energy decreases and the splitting process advances. As shown in Fig. 3b, the configuration of four split pieces is a more stable state than non-split single particle. Hence the non-split large particle is considered as a meta-stable phase and the activation barrier is the interfacial free energy corresponding to the increased interfacial area. Considering these points, the splitting phenomenon is kinetically rather than energetically driven.

In this work, it has been shown through phase field study that the interaction between the elastic field condensed at the concave region and solute enrichment causes the particle splitting phenomenon. The splitting is mediated by interface instability, which differs from previous studies that splitting occurs through the nucleation of the matrix phase at the center of the precipitate Miyazaki82 ; Wang93  or that it is caused by energy considerations Khachaturyan88 ; Doi92 . In particular, the sequential splitting behavior realized by our simulation shows an excellent agreement with previous experimental studies Miyazaki79 ; Yoo95 .


format: text

text:
# Twisted Klein curves modulo 2

## 1 Six inequivalent models modulo 2

The behaviour of the Klein curve is special in all the characteristics that divide 168. In this paper, we are particularly interested in the different models that exist in characteristic two. Theorem 1 shows that each model is characterised by a transformation $`PPSL(3,2)`$, such that the generic tangent of the curve at the point $`(x,y,z)`$ passes through $`(x^2,y^2,z^2)P^t`$ and $`(x^8,y^8,z^8)(P^t)^3`$. Two elements of $`PSL(3,2)`$ define equivalent curves if and only if they belong to the same conjugacy class. The conjugacy classes of $`PSL(3,2)`$ are of size 1 (identity), 24 (elements of order 7 with even trace), 24 (order 7, odd trace), 56 (order 3), 42 (order 4), and 21 (order 2). Up to equivalence we find six different models. For each of the six conjugacy classes, the corresponding curve has $`0`$ (identity), $`7`$ (order 7, even trace), $`0`$ (order 7, odd trace), $`3`$ (order 3), $`2`$ (order 4), or $`4`$ (order 2) rational points over the binary field.

The Klein curve modulo 2 represents a curve with $`3`$ rational points. In , Klein gives a different model

$`K^{}:{\displaystyle \frac{1}{21\sqrt[3]{7}}}(X^4+21X^2YZ147Y^2Z^2+49XY^3+49XZ^3)=0`$
which modulo 2 represents a curve with $`4`$ rational points (all points not on the line $`Y=Z`$). L.E. Dickson, in his paper “Quartic curves modulo 2” , classifies up to equivalence all the quartics that have $`0,7,6,`$ or $`5`$ rational points (for the last case, he finds twenty five equivalence classes and “Quartics with 1, 2, 3, or 4 real points have not been treated since they would probably not present sufficient novelty to compensate for the increased length of the investigation”). The paper provides us with three more inequivalent models,

$`\alpha :X^4+Y^4+Z^4+X^2Y^2+Y^2Z^2+Z^2X^2+X^2YZ+XY^2Z+XYZ^2=0`$
$`A:X^3Y+X^2Y^2+XZ^3+X^2Z^2+Y^3Z+YZ^3=0`$
$`\gamma _{1,0}:X^4+Y^4+Z^4+X^3Y+XZ^3+Y^3Z+XYZ^2=0`$
The model $`\alpha `$ is the unique model invariant under the full group $`PSL(3,2)`$. The model $`A`$ has the maximum number of rational points for a genus 3 curve. An equivalent model appears in , . We will denote it $`X_{N=7}`$.

$`X_{N=7}:XY(X+Y)(X+Z)+YZ(Z+Y)Y+ZX(X+Z)Z=0`$
The model $`\gamma _{1,0}`$, with no rational points, has trivial Jacobian. An equivalent model appears in the list by Madan and Queen  of the seven congruence function fields which have class number one and genus different from zero. We will denote it $`X_{h=1}`$.

$`X_{h=1}:Y^4+XY^3+(X^2+XZ)Y^2+(X^3+Z^3)Y+(X^4+XZ^3+Z^4)=0`$
The classes for $`A`$ and $`\gamma _{1,0}`$ are of same size and differ by a translation by $`\alpha `$. It happens that $`X_{N=7}=\gamma _{1,0}+\alpha `$, and $`X_{N=7}`$ and $`\gamma _{1,0}`$ have the same seven rational automorphisms. For the class with $`2`$ rational points, we give the model

$`X_{N=2}:X^4+Y^4+Z^4+X^3Y+XZ^3+Y^3Z+XYZ^2+`$
$`+X^2Y^2+X^2Z^2+Y^2Z^2+YZ^3=0`$
Comparison with $`\gamma _{1,0}`$ shows that the rational points are those not on $`X(Y+Z)=0`$. Properties of the different models are summarized in the Appendix.

## 2 Models in characteristic zero

The Klein curve arises in characteristic zero as a model for the elliptic modular curve $`X(7)`$. In particular this yields that the curve has the simple group $`PSL(2,7)`$ of order 168 as automorphism group. Many of the interesting properties of the Klein curve are directly related to this group. F. Klein derived the model

$$X^3Y+Y^3Z+Z^3X=0$$

as follows. A plane curve of degree four contains 24 flexpoints, that is points where the tangent intersects with multiplicity at least three. The group acts faithfully on the set of flexpoints and each flexpoint has stabilizer of order seven. There are eight subgroups of order seven, each fixing at most three points. It follows that each subgroup of order seven fixes three points, all of them flexpoints. The subgroup fixing a flexpoint also fixes the fourth point of intersection on its tangent, which must therefore be a flexpoint with the same stabilizer. Three flexpoints connected by their tangents form an oriented triangle, that Klein called a Wendedreiecke. The Klein quartic is now the curve with triangle of inflection tangents $`XYZ=0`$.

It is convenient to view the action of the automorphism group on the curve through its action on the eight triangles of inflection tangents. Each triangle has stabilizer of order 21. Let $`\zeta `$ be a primitive seventh root of unity. The stabilizer of the triangle $`XYZ=0`$ is generated by

$$\tau =\left(\begin{array}{ccc}\zeta & 0& 0\\ 0& \zeta ^4& 0\\ 0& 0& \zeta ^2\end{array}\right)$$

and the cyclic permutation $`\rho `$ of the coordinates $`X,Y,Z`$. It has 14 elements of order three. The full group has 56 elements of order three and each automorphism of order three appears in the stabilizer of two triangles. Say the other triangle stabilized by the cyclic permutation $`\rho `$ has a vertex $`(a:b:c)`$. Klein showed that

$$\sigma =\left(\begin{array}{ccc}a& b& c\\ b& c& a\\ c& a& b\end{array}\right)$$

acts on the curve, for

$$a=(\zeta ^5\zeta ^2)/\sqrt{7},b=(\zeta ^3\zeta ^4)/\sqrt{7},c=(\zeta ^6\zeta )/\sqrt{7},$$

with $`\sqrt{7}=\zeta +\zeta ^2\zeta ^3+\zeta ^4\zeta ^5\zeta ^6`$. The automorphism $`\sigma `$ is of order two and interchanges the two triangles. Indeed, $`\sigma `$ is of order two only if $`f=(xa)(xb)(xc)`$ is of the form $`x^3\pm x^2r`$, and it acts on the Klein curve only if $`r=0,\pm 1/7`$.

The full automorphism group $`PSL(2,7)`$ of the Klein curve is visible from its action on the eight triangles of inflection tangents. Let $`\mathrm{\Delta }_{\mathrm{}}`$ be the triangle $`XYZ=0`$, let $`\mathrm{\Delta }_0=\sigma \mathrm{\Delta }_{\mathrm{}}`$ be the other real triangle, and let $`\mathrm{\Delta }_i=\tau ^i\mathrm{\Delta }_0`$, for $`i=0,1,\mathrm{},6`$. The action of $`\sigma `$ and $`\tau `$ on the triangles is described by the permutations $`(\mathrm{1\hspace{0.17em}6})(\mathrm{2\hspace{0.17em}3})(\mathrm{4\hspace{0.17em}5})(0\mathrm{})`$ and $`(\mathrm{1\hspace{0.17em}2\hspace{0.17em}3\hspace{0.17em}4\hspace{0.17em}5\hspace{0.17em}6\hspace{0.17em}0})`$ respectively. And an isomorphism to $`PSL(2,7)`$ is given by

$$\tau \left(\begin{array}{cc}1& 1\\ 0& 1\end{array}\right)\sigma \left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)$$

Burnside pointed out that $`\tau `$ and $`\sigma `$ generate the automorphism group. Klein had originally included the cyclic permutation $`\rho `$ in the generating set.

The two triangles of inflection tangents that are left invariant by the cyclic permutation $`\rho `$ are the only two of the eight triangles with real vertices. If we change the orientation in the triangle $`XYZ=0`$, we obtain the different model

$$X^3Z+Y^3X+Z^3Y=0.$$

If we change the orientation in the other real triangle, we obtain the different model

$$\begin{array}{c}O_4:\frac{1}{7}((X^4+Y^4+Z^4)3(X^2Y^2+Y^2Z^2+Z^2X^2)\hfill \\ \hfill +3(X^2YZ+XY^2Z+XYZ^2)+6(X^3Z+Y^3X+Z^3Y))=0.\end{array}$$

This model was derived independently as Equation (1.22) in . After reduction modulo $`p=2`$, the model agrees with the invariant model $`\alpha `$ and has all automorphisms defined over the base field. Ciani  gives a model over $`Q(\sqrt{7})`$

$$C_4:X^4+Y^4+Z^4+3\omega (Y^2Z^2+Z^2X^2+X^2Y^2)=0,\omega ^2+\omega +2=0.$$

It coincides with Equation (1.11) in . The following model () has the $`S_3`$ symmetry of $`C_4`$ and the same binary reduction as $`O_4`$,

$$A_4:7s_1(s_1^3+s_3)(2s_1^2+s_2)^2=0.$$

All of the last three curves have their automorphisms defined over $`Q(\sqrt{7})`$. The six real flexpoints of the last curve have $`s_1=1,s_2=2,s_3=+1`$, i.e. have coordinates $`a,b,c`$ such that, for $`\theta =2\pi /7`$,

$$\begin{array}{c}(xa)(xb)(xc)=x^3+x^22x1=\hfill \\ \hfill =(x2\mathrm{cos}(\theta ))(x2\mathrm{cos}(2\theta ))(x2\mathrm{cos}(4\theta ))\end{array}$$

## 3 Bitangents in even characteristic

In characteristic zero, a non-singular quartic curve admits 28 bitangents. In even characteristic, the number of bitangents is at most $`7`$. In that case, the bitangents define a Fano plane. All non-singular quartic curves with seven bitangents are in the same isomorphism class. The unique curve with bitangents $`X,Y,Z,X+Y,Y+Z,X+Y+Z,X+Z`$ is the invariant curve

$$\alpha :X^4+Y^4+Z^4+Y^2Z^2+Z^2X^2+X^2Y^2+X^2YZ+Y^2ZX+Z^2XY=0$$

Uniqueness follows for example by considering that the bitangents must touch the curve in the fourteen points of $`P^2(F_4)`$ not in $`P^2(F_2)`$. Other curves with seven bitangents will have their bitangents defined over some extension of $`F_2`$. In those cases, the equation of the curve can be recovered from the bitangents by substitution in $`\alpha `$ of three non-concurrent bitangents.

The Klein curve has bitangents

$$K:a^2X+aY+a^4Z=0,a^8+a=0,a0.$$

It happens that the models chosen by L.E. Dickson have a similar description of the bitangents. The model $`A`$ has bitangents

$$A:b^2X+bY+b^4Z=0,b^8+b^2+b=0,b0.$$

The model $`\gamma _{1,0}`$ has bitangents

$$\gamma _{1,0}:c^2X+cY+c^4Z=0,c^8+c^4+c=0,c0.$$

If we choose different but equivalent models, the bitangents change but are still parametrized by the same values. The curve $`X_{N=7}`$, which is equivalent to $`A`$, has bitangents

$$X_{N=7}:bX+b^{16}Y+b^8Z=0,b^8+b^2+b=0,b0.$$

The curve $`X_{h=1}`$, which is equivalent to $`\gamma _{1,0}`$, has bitangents

$$X_{h=1}:cX+c^8Y+c^2Z=0,c^8+c^4+c=0,c0.$$

We chose the model $`X_{N=2}`$ such that it has bitangents

$$X_{N=2}:d^2X+dY+d^4Z=0,d^8+d^4+d^2+d=0,d0.$$

For the invariant curve, the bitangents have coefficients, in matrix form,

$$\alpha :\left[\begin{array}{ccccccc}1& 0& 0& 1& 0& 1& 1\\ \multicolumn{7}{c}{}\\ 0& 1& 0& 1& 1& 1& 0\\ \multicolumn{7}{c}{}\\ 0& 0& 1& 0& 1& 1& 1\end{array}\right]$$

So that the coefficients appearing in a given row are the zeros of $`(x^8+x^4)/x`$. After reduction, Klein’s model $`K^{}`$ has bitangents, with $`F_4=\{0,1,\omega ,\overline{\omega }\}`$,

$$K^{}:\left[\begin{array}{ccccccc}1& 1& 0& 0& 1& 0& 1\\ \multicolumn{7}{c}{}\\ 0& \omega & \overline{\omega }& \omega & 1& 1& \overline{\omega }\\ \multicolumn{7}{c}{}\\ 0& \overline{\omega }& \omega & \overline{\omega }& 1& 1& \omega \end{array}\right]$$

Substituion of $`Z=Z+X`$ in $`K^{}`$ gives the equivalent reduced model

$`X_{N=4}:X^2YZ+X^2Y^2+X^2Z^2+Y^2Z^2+XY^3+XZ^3+X^3Y+X^3Z=0`$
with bitangents,

$$X_{N=4}:\left[\begin{array}{ccccccc}1& \omega & \omega & \overline{\omega }& 0& 1& \overline{\omega }\\ \multicolumn{7}{c}{}\\ 0& \omega & \overline{\omega }& \omega & 1& 1& \overline{\omega }\\ \multicolumn{7}{c}{}\\ 0& \overline{\omega }& \omega & \overline{\omega }& 1& 1& \omega \end{array}\right]$$

The coefficients appearing in a given row are the zeros of $`(x^8+x^2)/x`$. The only non-trivial additive polynomial that can not be used to define coefficients of a Fano plane is $`x^8+x^4+x^2`$.

In all cases, the set of bitangents is defined over $`F_2`$. The Frobenius automorphism acts on the set of bitangents. For a careful choice of the coefficients of the bitangents (such that the linear expressions for the bitangents are closed under addition), the action is linear and can be represented by an element of $`PSL(3,2)`$.

###### Definition 1

To a quartic curve over $`F_2`$ with seven bitangents, we associate the unique element $`RPSL(3,2)`$ such that the action of $`R`$ on the bitangents agrees with the action of the Frobenius automorphism: $`(a^2,b^2,c^2)^T=R(a,b,c)^T`$, for each bitangent $`ax+by+cz=0`$.

In the next section, we give the inverse map, which gives for each $`RPSL(3,2)`$ the unique quartic curve such that the Frobenius action on the bitangents corresponds with the action of $`R`$.

## 4 Description of the 168 twisted curves

The group $`PSL(3,2)`$ is the full automorphism group of the curve $`\alpha `$. In , L.E. Dickson gives the modular invariants of the general linear group over a finite field. In particular, he shows that $`\alpha `$ is the smallest projective invariant of $`PSL(3,2)`$ acting on $`F_2[X,Y,Z]`$. The other generating invariants are of degree six and seven, and the three invariants are algebraicly independent , . Consider the additive polynomial in the variable $`T`$ with coefficients in $`F_2[X,Y,Z]`$ that cancels at $`X,Y,Z`$,

$$\underset{vX,Y,Z}{}(T+v)=T^8+I_4(X,Y,Z)T^4+I_6(X,Y,Z)T^2+I_7(X,Y,Z)T.$$

The algebra $`F_2[I_4,I_6,I_7]`$ is the ring of $`PSL(3,2)`$ invariants in $`F_2[X,Y,Z]`$. The coefficients follow by solving

$`\left(\begin{array}{ccc}X& X^2& X^4\\ Y& Y^2& Y^4\\ Z& Z^2& Z^4\end{array}\right)\left(\begin{array}{c}I_7\\ I_6\\ I_4\end{array}\right)=\left(\begin{array}{c}X^8\\ Y^8\\ Z^8\end{array}\right).`$
For $`I_4`$, we have

$`I_4(X,Y,Z)=det\left(\begin{array}{ccc}X& X^2& X^8\\ Y& Y^2& Y^8\\ Z& Z^2& Z^8\end{array}\right)/det\left(\begin{array}{ccc}X& X^2& X^4\\ Y& Y^2& Y^4\\ Z& Z^2& Z^4\end{array}\right).`$
Thus, $`I_4`$ is the unique non-trivial plane curve (i.e. different from a line), for which the line through a point $`(x:y:z)`$ and its Frobenius $`F(x:y:z)=(x^2:y^2:z^2)`$ also contains $`F^3(x:y:z)=(x^8:y^8:z^8)`$.

Let $`C`$ be a curve isomorphic to $`I_4`$. Then $`C`$ can be represented as a plane quartic. The canonical divisor of a quartic curve is given by the intersection of a line with the curve. Therefore the isomorphism preserves collinear points and can be represented by a linear transformation. Let $`(x,y,z)A(x,y,z)=(x^{},y^{},z^{})`$ be an isomorphism from $`C`$ to $`I_4`$. We have that $`(x,y,z)`$ is a point on $`C`$ if and only if $`(x^{},y^{},z^{})`$ is a point on $`I_4`$ if and only if

$$(x^{},y^{},z^{}),F(x^{},y^{},z^{}),F^3(x^{},y^{},z^{})$$

are collinear if and only if

$$A(x,y,z),F(A(x,y,z)),F^3(A(x,y,z))$$

are collinear if and only if

$$(x,y,z),(A^1FA)(x,y,z),(A^1FA)^3(x,y,z)$$

are collinear. We claim that $`(A^1FA)(x,y,z)=(PF)(x,y,z)`$, for some $`PPSL(3,2)`$. Let $`A^{(2)}`$ be the isomorphism $`A`$ with coefficients squared. We have

$$(A^1FA)(x,y,z)=(A^1A^{(2)}F)(x,y,z)$$

and $`A^{(2)}A^1\text{Aut}(I_4)=PSL(3,2)`$, say $`A^{(2)}A^1=Q`$. The element $`Q`$ is determined by the curve $`C`$ only up to conjugacy, for if we replace $`A`$ with $`SA`$, for $`S\text{Aut}(I_4)=PSL(3,2)`$, then $`(SA)^{(2)}(SA)^1=SQS^1`$. On the other hand, $`A^1QA=A^1A^{(2)}`$ does not change when $`A`$ is replaced with $`SA`$. The latter matrix has coefficients in $`F_2`$, for $`Q^2=A^{(2)}A^1Q=QA^{(2)}A^1`$, or

$$A^1QA=(A^{(2)})^1QA^{(2)}=(A^1QA)^{(2)}$$

and $`A^1A^{(2)}=A^1QA`$ belongs to $`PSL(3,2)`$.

###### Theorem 1

A curve is isomorphic to the invariant curve $`I_4`$ if and only if it is of the form

$$det(\left(\begin{array}{c}X\\ Y\\ Z\end{array}\right)P\left(\begin{array}{c}X^2\\ Y^2\\ Z^2\end{array}\right)P^3\left(\begin{array}{c}X^8\\ Y^8\\ Z^8\end{array}\right))/det(\left(\begin{array}{c}X\\ Y\\ Z\end{array}\right)P\left(\begin{array}{c}X^2\\ Y^2\\ Z^2\end{array}\right)P^2\left(\begin{array}{c}X^4\\ Y^4\\ Z^4\end{array}\right))$$

for some $`PPSL(3,2)`$. The rational automorphisms of such a curve are given by the elements of $`PSL(3,2)`$ that commute with $`P`$. The action of the Frobenius on the bitangents, as in Definition 1, is given by the matrix $`R=P^t`$. The tangent at the point $`(x,y,z)`$ has intersection divisor

$$2(x,y,z)+(x^2,y^2,z^2)P^t+(x^8,y^8,z^8)(P^t)^3$$

In particular, $`PPSL(3,2)`$ is determined by the generic tangent of the curve.

Proof. The curve is uniquely determined by the matrix $`P=A^1A^{(2)}`$, for any isomorphism $`A`$ from $`C`$ to $`I_4`$. Now $`SPSL(3,2)`$ defines an automorphism of $`C`$ if and only if $`AS`$ is an isomorphism from $`C`$ to $`I_4`$ if and only if $`(AS)^1(AS)^{(2)}=S^1PS=P`$. Let the columns of the matrix $`B`$ contain the coefficients of three non concurrent bitangents. Then $`B^{(2)}=RB`$. On the other hand, $`B^t`$ defines an isomorphism from $`C`$ to $`I_4`$. And $`P=(B^t)^1(B^t)^{(2)}=R^t`$. qed

The general construction of twisted models for a given object with non-trivial automorphisms is given by elementary galois cohomology. The case that we consider here seems to be particularly favourable for an explicit treatment. The invariant curve $`I_4`$ and all its automorphisms are defined over the base field, and the equation of the curve $`I_4`$ is a simple relation among Frobenius images of the generic point on the curve.

As illustration of the theorem, we may write for the Klein curve,

$$X^3Y+Y^3Z+Z^3X=det\left(\begin{array}{ccc}X& Y^2& X^8\\ \multicolumn{3}{c}{}\\ Y& Z^2& Y^8\\ \multicolumn{3}{c}{}\\ Z& X^2& Z^8\end{array}\right)/det\left(\begin{array}{ccc}X& Y^2& Z^4\\ \multicolumn{3}{c}{}\\ Y& Z^2& X^4\\ \multicolumn{3}{c}{}\\ Z& X^2& Y^4\end{array}\right).$$

## 5 Zeta functions for the twisted curves

The forms $`I_4`$ and $`I_6`$ define a pair of dual curves (up to an inseparable morphism of degree two), and are birationally equivalent. For both curves, all points on the tangent at $`P`$ are images of $`P`$ under the Frobenius automorphism $`F`$. And in both cases, the zeta function is determined by the intersection divisor of a generic tangent. For the tangent $`L_P`$ of $`I_4`$ at $`P`$,

$$L_PI_4=P+P+F(P)+F^3(P),$$

and the Frobenius eigenvalues are zeros of

$$(2+t+t^3)=(1+t)(2t+t^2).$$

Similarly, for the tangent $`L_P`$ of $`I_6`$ at $`P`$,

$$L_PI_6=P+P+P+P+F^2(P)+F^3(P),$$

and the Frobenius eigenvalues are zeros of

$$(4+t^2+t^3)=(2+t)(2t+t^2).$$

Thus, both $`I_4`$ and the normalization of $`I_6`$ have zeta polynomial

$$L(t)=(1t+2t^2)^3.$$

For the twisted models we see, either by elimination of $`P^t`$ in the generic tangent (given in Theorem 1) or by using properties of the zeta function, that the Frobenius eigenvalues are reciprocal zeros of

$$\underset{\zeta ^m=1}{}(1\zeta t+2\zeta ^2t^2)$$

where $`m`$ is the order of the matrix $`P`$. For the relevant values $`m=2,3,4,7,`$ we find the zeta polynomial immediately from the factorization of the above product into irreducibles and the known number of rational points. Let $`z^+=(1t+2t^2)`$ and $`z^{}=(1+t+2t^2)`$.

$$\begin{array}{ccccc}m=2:\hfill & & 1+3t^2+4t^4=z^+z^{},\hfill & & L(t)=z^+(z^{})^2.\hfill \\ m=4:\hfill & & 1t^4+16t^8=z^+z^{}z_2,\hfill & & L(t)=(z^+)^2z^{}.\hfill \\ m=3:\hfill & & 1+5t^3+8t^6=z^+z_3,\hfill & & L(t)=z^+z_3.\hfill \\ m=7:\hfill & & 1+13t^7+128t^{14}=z^+z_{7}^{}{}_{}{}^{+}z_{7}^{}{}_{}{}^{}.\hfill & & L(t)=z_{7}^{}{}_{}{}^{+},z_{7}^{}{}_{}{}^{}\hfill \end{array}$$

The two irreducible factors $`z_{7}^{}{}_{}{}^{+},z_{7}^{}{}_{}{}^{}`$ of degree six represent the two different conjugacy classes for $`P`$ of order $`7`$.

$$\begin{array}{ccc}\text{even trace:}\hfill & & L(t)=z_{7}^{}{}_{}{}^{+}=1+4t+9t^2+15t^3+18t^4+16t^5+8t^6\hfill \\ \text{odd trace:}\hfill & & L(t)=z_{7}^{}{}_{}{}^{}=13t+2t^2+t^3+4t^412t^5+8t^6\hfill \end{array}$$

The factor $`z^+`$ corresponds, as is well-known , to the elliptic curve with complex multiplication by $`7`$. In the cases $`m=2,3,4`$, it appears as the quotient of the twist by a rational automorphism. In the case $`m=7`$, the quotient gives the projective line. We give the elliptic curve as the quotient, in characteristic zero, of the Klein curve by the cyclic permutation of the coordinates. Let $`B`$ be the divisor sum of the two points in the intersection of the Klein curve and the invariant bitangent $`X+Y+Z=0`$. Recall that the Klein curve has two triangles of inflection tangents that are invariant under the cyclic permutation. Let $`\mathrm{\Delta }_{\mathrm{}}`$ and $`\mathrm{\Delta }_0`$, respectively, be the divisor sums of their vertices, respectively. Let $`x,y,zk[X,Y,Z]`$ intersect the Klein curve as follows

$$(x)=2\mathrm{\Delta }_{\mathrm{}}+2\mathrm{\Delta }_0,(y)=\mathrm{\Delta }_{\mathrm{}}+\mathrm{\Delta }_0+3B,(z)=4\mathrm{\Delta }_{\mathrm{}}.$$

We obtain

$`x=X^2Y+Y^2Z+Z^2X+XYZ=((X^3Y+Y^3Z+Z^3X)+s_s^2)/s_1,`$
$`y=(X+Y+Z)(XY+YZ+ZX),z=XYZ.`$
The morphism $`(x:y:z)`$ maps $`K`$ onto the elliptic curve

$$E:y^2+5xy=x^3x^2+7x.$$

In fact,

$$\begin{array}{c}y^2z5xyz+x^3+x^2z+7xz^2=\hfill \\ \hfill =(X^3Y+Y^3Z+Z^3X)(X^3Y^2+Y^3Z^2+Z^3X^2).\end{array}$$

Replacing $`x=((X^3Y+Y^3Z+Z^3X)+s_s^2)/s_1`$ by $`x=s_2^2/s_1`$ in the morphism, we find, as in ,

$$y^2z5xyz+x^3+x^2z+7xz^2=(X^3Y+Y^3Z+Z^3X)(X^3Z+Y^3X+Z^3Y)$$

The triangles $`\mathrm{\Delta }_{\mathrm{}}`$ and $`\mathrm{\Delta }_0`$ are mapped to the point at infinity and the origin respectively. The only ramification is at $`B`$.

## Appendix A Summary of curves

The curves appearing in this paper are summarized in the following table. It gives the number of equivalent curves in a class, the structure of the subgroup of rational automorphisms, and properties of the matrix $`P`$. The matrix $`P`$ for each of the curves follows (as determined by Theorem 1).

| Curve | $`\mathrm{\#}`$ | Rational | $`P`$ | |
| --- | --- | --- | --- | --- |
|  |  | automorphisms | order | trace |
| $`\alpha =I_4`$ | 1 | $`PSL(3,2)`$ | 1 | 1 |
| $`A`$, $`X_{N=7}`$ | 24 | $`Z/7Z`$ | 7 | 0 |
| $`\gamma _{1,0}`$, $`X_{h=1}`$ | 24 | $`Z/7Z`$ | 7 | 1 |
| $`K`$ | 56 | $`Z/3Z`$ | 3 | 0 |
| $`X_{N=2}`$ | 42 | $`Z/4Z`$ | 4 | 1 |
| $`K^{}`$, $`X_{N=4}`$ | 21 | $`D_4`$ | 2 | 1 |

$$\begin{array}{ccccccc}\alpha \hfill & \left[\begin{array}{ccc}1& 0& 0\\ & \multicolumn{2}{c}{}\\ 0& 1& 0\\ & \multicolumn{2}{c}{}\\ 0& 0& 1\end{array}\right]& A\hfill & \left[\begin{array}{ccc}0& 1& 1\\ & & \\ 0& 0& 1\\ & & \\ 1& 0& 0\end{array}\right]& X_{N=7}\hfill & \left[\begin{array}{ccc}1& 1& 0\\ & & \\ 0& 1& 1\\ & & \\ 1& 0& 0\end{array}\right]& \\ K\hfill & \left[\begin{array}{ccc}0& 1& 0\\ & \multicolumn{2}{c}{}\\ 0& 0& 1\\ & \multicolumn{2}{c}{}\\ 1& 0& 0\end{array}\right]& \gamma _{1,0}\hfill & \left[\begin{array}{ccc}0& 1& 0\\ & & \\ 0& 0& 1\\ & & \\ 1& 0& 1\end{array}\right]& X_{h=1}\hfill & \left[\begin{array}{ccc}0& 0& 1\\ & & \\ 0& 1& 1\\ & & \\ 1& 1& 0\end{array}\right]& \\ X_{N=2}\hfill & \left[\begin{array}{ccc}0& 1& 1\\ & \multicolumn{2}{c}{}\\ 0& 0& 1\\ & \multicolumn{2}{c}{}\\ 1& 0& 1\end{array}\right]& K^{}\hfill & \left[\begin{array}{ccc}1& 0& 0\\ & & \\ 1& 0& 1\\ & & \\ 1& 1& 0\end{array}\right]& X_{N=4}\hfill & \left[\begin{array}{ccc}1& 0& 0\\ & & \\ 1& 0& 1\\ & & \\ 1& 1& 0\end{array}\right]& \end{array}$$

Current address:
Department of Mathematics
University of Illinois at Urbana-Champaign
1409 W. Green Street
Urbana, IL 61801
E-mail: duursma@math.uiuc.edu


text:
# Poland-Scheraga models and the DNA denaturation transition

## 1 Introduction

When a solution of DNA is heated, the double stranded molecules denature into single strands. In this process, looping out of AT rich regions of the DNA segments first occurs, followed eventually by separation of the two strands as the paired segments denature. This denaturation process corresponds to a phase transition .

A simple model of the DNA denaturation transition was introduced in 1966 by Poland and Scheraga  (hereinafter referred to as PS) and refined by Fisher . The model consists of an alternating sequence (chain) of straight paths and loops, which idealize denaturing DNA, consisting of a sequence of double stranded and single stranded molecules. An attractive energy is associated with paths. Interactions between different parts of a chain and, more generally, all details regarding real DNA such as chemical composition, stiffness or torsion, are ignored. It was found that the phase transition is determined by the critical exponent $`c`$ of the underlying loop class. Due to the tractability of the problem of random loops, that version of the problem was initially studied by PS . The model displays a continuous phase transition in both two and three dimensions. It was argued  that replacing random loops by self-avoiding loops, suggested as a more realistic representation accounting for excluded volume effects within each loop, sharpens the transition, but does not change its order.

However, the sharp jumps observed in the UV absorption rate in DNA melting111 In recent years, a number of other properties of DNA molecules have been studied by refined techniques such as optical tweezers and atomic force microscopy, and theoretical descriptions of underlying effects such as unzipping  have been proposed. These will not be discussed here. experiments , which correspond to a sudden breaking of large numbers of base pairs, indicated that a first order phase transition would be the appropriate description. The question whether such an asymptotic description, which implies very long chains, is valid for relatively short DNA sequences, has been discussed recently . Nevertheless, a directed extension of the PS model, being essentially a one-dimensional Ising model with statistical weighting factors for internal loops, is widely used today and yields good coincidence of simulated melting curves with experimental curves for known DNA sequences . Another recent application of PS models analyses the role of mismatches in DNA denaturation . A numerical approach to DNA denaturation, which we will not discuss further, uses variants of the Peyrard-Bishop model , a Hamiltonian model of two harmonic chains coupled by a Morse potential.

With the advent of efficient computers, it has more recently been possible to simulate analytically intractable models extending the PS class, which are assumed to be more realistic representations of the biological problem. One of these is a model of two self-avoiding and mutually avoiding walks, with an attractive interaction between different walks at corresponding positions in each walk . The model exhibits a first order phase transition in $`d=2`$ and $`d=3`$. The critical properties of the model are described by an exponent $`c^{}`$ related to the loop length distribution , see also Fisher’s review article . This exponent is called $`c`$ again. Indeed, for PS models, it coincides with the loop class exponent $`c`$ if $`1<c<2`$, see below. Within a refined model, where different binding energies for base pairs and stiffness are taken into account, the exponent $`c^{}`$ seems to be largely independent of the specific DNA sequence and of the stiffness of paired walk segments corresponding to double stranded DNA parts . There are, however, no simulations of melting curves for known DNA sequences which are compared to experimental curves for this model.

An approximate analytic derivation of the exponent related to this new model was given by Kafri et al.  using the theory of polymer networks. They estimated the excluded volume effect arising from the interaction between a single loop and two attached walks. This approach (refined recently ) yields an approximation to the loop length distribution exponent $`c^{}`$, which agrees well with simulation results of interacting self-avoiding walk pairs . There is a recent debate about the relevance of this approximation to real DNA . The polymer network approach, as initiated previously , led to a number of related applications .

In this article, we reconsider PS models for three principal reasons. Firstly, the older articles are short in motivating the use of particular loop classes, which may be misleading in drawing conclusions about the thermodynamic effects of different loop classes. In fact, the loop classes discussed in the early approaches  are classes of rooted loops and lead to chains which are not self-avoiding. This seems unsatisfactory from a biological point of view, since real DNA is self-avoiding. Secondly, the common view holds that PS models with self-avoiding loops cannot display a first order transition in two or three dimensions. In fact, this view led to extending the PS class  in order to find a model with a first order transition. However this view is incorrect, as we demonstrate by a self-avoiding PS model with self-avoiding loops. Thirdly, the two exponents $`c`$ and $`c^{}`$, extracted from different expressions as described above, are used in the literature without distinction, although there are subtle differences, which we will point out.

This paper is organized as follows. In the next section, we give a rigorous discussion of PS models and their phase transitions. We derive the scaling functions which describe the behavior below the critical temperature, as the critical point is approached. We then derive the loop statistics, thereby analyzing the occurrence of the loop class exponent in the loop length distribution. The third section reviews the prevailing PS models, with emphasis on motivation for the underlying loop classes. We then introduce a self-avoiding PS model displaying a first order phase transition in two and three dimensions. In Section 5, exactly solvable directed PS models are discussed. We will give explicit expressions for generating functions. The critical behavior of some of these models has been analyzed previously by different methods . This is followed by a discussion of models extending the PS class. We conclude with a discussion of some open questions.

## 2 Poland-Scheraga models: general formalism

We define PS models and consider analytical properties of their free energies, thereby refining previous expositions  in a rigorous manner. We analyze phase transitions, extract the asymptotic behavior of chains at temperatures near the phase transition by means of scaling functions, and consider the loop length distribution of chains.

### 2.1 Definition of the model

We consider a discrete model, defined on the hypercubic lattice $`^d`$. Double stranded DNA segments are modeled by paths, and single stranded DNA segments are modeled by loops on the edges of the lattice. Each loop is assumed to have two marked vertices to indicate where paths are attached. (Like PS, we do not a priori assume the “DNA condition” that the marked vertices divide the loop into parts of equal lengths. There is a recent discussion about the effect of mismatches .) Any alternating sequence of paths and (marked) loops is called a chain. A PS model consists of all chains obtained by concatenation of paths and loops from a given path class and a given loop class, where the initial segment and the final segment of a chain are both paths. (The case of an open end affects the behavior above the critical temperature  and will not be considered here ). Note that, in general, such chains are not self-avoiding, in contrast to real DNA. Self-avoidance may be violated by paths, by loops, or by the way segments are concatenated. A chain is called segment-avoiding if there are no overlaps, i.e., every two non-neighboring segments have no vertex in common, and every two neighboring segments have exactly the marked vertex in common. We call a PS model self-avoiding if paths and loops are self-avoiding and if all chains of the model are segment-avoiding. We will discuss below several examples of self-avoiding PS models, which arise from self-avoiding walks.

The requirement of self-avoidance restricts the admissible path classes and loop classes. A simple subclass of self-avoiding PS models are directed PS models: We call a chain directed if there is a preferred direction such that the order of the chain segments induces the same order on the vertex coordinates (w.r.t. the preferred direction), for each pair of vertices taken from two different chain segments. Such chains are then segment-avoiding. We call a PS model directed if paths, loops and chains are directed. We will discuss below several examples of directed PS models which arise from classes of directed walks.

For a given PS model, let $`z_{m,n}`$ denote the number of chain configurations with $`m`$ contacts and length $`n`$. The generating function is defined by

$$Z(x,w)=\underset{m,n}{}z_{m,n}w^mx^n,$$
(2.1)
where the activity $`x`$ is conjugate to the chain length $`n`$. The Boltzmann factor $`w=e^{E/kT}`$ takes into account the attractive interaction (achieved by setting the energy $`E<0`$) between bonds. $`T`$ is the temperature, and $`k`$ is Boltzmann’s constant. Over the relevant temperature range $`0<T<\mathrm{}`$, we have $`\mathrm{}>w>1`$. Note that there is no interaction between different segments in a chain.

The generating function $`Z(x,w)`$ can be expressed in terms of the generating functions for paths $`V(x)`$ and loops $`U(x)`$. These are

$$V(x)=\underset{n=0}{\overset{\mathrm{}}{}}b_nx^n,U(x)=\underset{n=1}{\overset{\mathrm{}}{}}p_{2n}x^n,$$
(2.2)
where $`b_n`$ is the number of paths of length $`n`$, and $`p_{2n}`$ is the number of loops222 In simulations of melting curves like MELTSIM , a cooperative parameter $`\sigma 10^5`$ is used. In our context, this amounts to replacing the loop generating function $`U(x)`$ by $`\sigma U(x)`$.  of length $`2n`$. Due to the chain structure, we get a geometric series in $`V(wx)U(x)`$,

$$Z(x,w)=\frac{V(wx)}{1U(x)V(wx)}=\underset{n=1}{\overset{\mathrm{}}{}}Z_n(w)x^n.$$
(2.3)
Since we want to analyze phase transitions of the model, which can only occur in the infinite system, we define the free energy of the model as

$$f(w)=\underset{n\mathrm{}}{lim}\frac{1}{n}\mathrm{log}Z_n(w)=\mathrm{log}x_c(w),$$
(2.4)
where, for fixed $`w`$, $`x_c(w)`$ is the radius of convergence of $`Z(x,w)`$. Concatenation arguments and supermultiplicative inequalities can be used to show that the free energy exists . We will alternatively investigate properties of the free energy in terms of the generating functions for paths and loops333 The reader may find it illuminating in following this very general discussion that now follows to refer to a specific model, discussed in Section 5.1.. Throughout the paper, we employ the following assumption.

###### Assumption 1.

Assume that the generating functions $`U(x)`$ and $`V(x)`$, defined in (2.2), have radius of convergence $`x_U`$ and $`x_V`$, respectively, where $`0<x_U<x_V1`$. At the critical point $`x=x_V`$, assume that $`V(x_V^{}):=lim_{xx_V^{}}V(x)=\mathrm{}`$.

Remark. $`U(x)`$ is a generating function (i.e., a series with non-negative coefficients), such that $`U(x)`$ and the derivative $`U^{}(x)`$ are strictly positive for $`0<x<x_U`$. A corresponding statement holds for $`V(x)`$. The assumption that $`0<x_U<1`$ reflects the requirement that the number of configurations $`p_{2n}`$ grows exponentially in length. The models commonly discussed  have $`x_V=1`$ and $`V(0)=1`$. For loops and paths of the same type, we have $`x_U=x_V^2<x_V<1`$, see also the examples discussed below. The assumption $`V(x_V^{})=\mathrm{}`$ is satisfied for typical classes of paths. Note that typically $`U(x_U^{}):=lim_{xx_U^{}}U(x)<\mathrm{}`$, but there are loop classes where $`U(x_U^{})=\mathrm{}`$ such as convex polygons . We have $`U(0)=0`$ by definition.

The radius of convergence $`x_c(w)`$ of $`Z(x,w)`$ is the minimum of the radius of convergence of $`U(x)V(wx)`$ and the point $`x_1(w)`$ where the denominator in (2.3) vanishes. Define $`F(x,w):=U(x)V(wx)`$. As a function of argument $`x`$, $`F(x,w)`$ is continuous and monotonically increasing for $`0<x<\mathrm{min}(x_U,x_V/w)`$. Note that $`F(0,w)=0`$. If $`wx_V/x_U`$, then $`lim_{xx_V/w}F(x,w)=\mathrm{}`$, such that there exists a unique solution $`x_1(w)<x_U`$ with $`F(x_1(w),w)=1`$. Assume that $`1w<x_V/x_U`$. If $`U(x_U^{})=\mathrm{}`$, there exists a unique solution $`x_1(w)<x_U`$. If $`U(x_U^{})<\mathrm{}`$ and $`U(x_U^{})V(x_U)1`$, we have $`F(x_U,w)1`$ for all $`w1`$, i.e., there exists a unique solution $`x_1(w)<x_U`$. If $`U(x_U^{})<\mathrm{}`$ and $`U(x_U^{})V(x_U)<1`$, we define $`w_c>1`$ by the condition $`U(x_U^{})V(w_1x_U)=1`$. Then, for $`ww_c`$, we have a unique solution $`x_1(w)<x_U`$. If $`w<w_c`$, there is no such solution, and the dominant singularity of $`Z(x,w)`$ occurs at $`x=x_U`$. We thus proved the following theorem.

###### Theorem 1.

If Assumption 1 is satisfied, the PS model (2.3) has free energy

$$f(w)=\{\begin{array}{cc}\mathrm{log}x_c(w)\hfill & (w_cw<\mathrm{})\hfill \\ \mathrm{log}x_U\hfill & (1<w<w_c),\hfill \end{array}$$
(2.5)
where the radius of convergence $`x_c(w)`$ of $`Z(x,w)`$ is, for $`w_cw<\mathrm{}`$, the unique positive solution of

$$U(x_c(w))V(wx_c(w))=1.$$
(2.6)
If $`U(x_U^{}):=lim_{xx_U^{}}U(x)=\mathrm{}`$ or $`U(x_U^{})V(x_U)1`$, we have $`w_c=1`$. Otherwise, $`w_c>1`$ is implicitly given by $`U(x_U^{})V(w_cx_U)=1`$. ∎

Remark. If $`w_c>1`$, this point is a critical point, see the following section.

### 2.2 Phase transitions

For $`w>w_c`$, there is no phase transition, since the radius of convergence $`x_c(w)>0`$ is analytic for $`w_c<w<\mathrm{}`$, as can be inferred from (2.6). A necessary condition for a phase transition at a finite temperature $`w_c>1`$ is $`U(x_U^{})V(x_U)<1`$. The nature of the transition at $`w=w_c>1`$ is exhibited by the fraction of shared bonds $`\theta (w)=w\frac{df(w)}{dw}`$, which is defined for $`ww_c`$. Note that $`\theta (w)=0`$ for $`w<w_c`$. For $`w>w_c`$, we have

$$\theta (w)=w\frac{x_c^{}(w)}{x_c(w)}=w\left(\frac{U^{}}{U}\frac{V}{V^{}}+w\right)^1>0.$$
(2.7)
Consider the limit $`ww_c^+`$. We have $`0<V<\mathrm{}`$, since $`UV=1`$ and $`U`$ is finite. Since $`V(x)`$ is a generating function, this implies $`0<V^{}<\mathrm{}`$. If $`U^{}(x_U^{})=\mathrm{}`$, it follows that $`\theta (w)0`$, such that the phase transition is continuous. If $`U^{}(x_U^{})<\mathrm{}`$, $`\theta (w)\theta _c>0`$, such that the phase transition is of first order. This leads to the following statement.

###### Theorem 2.

Let Assumption 1 be satisfied. If $`U(x_U^{})=\mathrm{}`$ or $`U(x_U^{})V(x_U)1`$, the PS model (2.3) has no phase transition at finite temperature. Otherwise, if $`U(x_U^{})<\mathrm{}`$ and $`U^{}(x_U^{})=\mathrm{}`$, a continuous phase transition will occur at $`w=w_c>1`$ defined in Theorem 1. If $`U(x_U^{})<\mathrm{}`$ and $`U^{}(x_U^{})<\mathrm{}`$, a first order phase transition will occur at $`w=w_c>1`$ defined in Theorem 1. ∎

Remark. The phase transition condition $`U(x_U^{})V(x_U)<1`$ is typically satisfied in more realistic models, where $`U(x)`$ is multiplied by the cooperativity parameter $`\sigma 10^5`$.

We conclude that the nature of the transition is determined by the singularity of the loop generating function $`U(x)`$ at $`x=x_U`$. It can be more directly related to the asymptotic properties of loops, which are typically of the form444 $`A_nB_n`$ for $`n\mathrm{}`$ means that $`lim_n\mathrm{}A_n/B_n=1`$. Similarly, $`f(x)g(x)`$ for $`xx_c`$ means that $`lim_{xx_c}f(x)/g(x)=1`$.

$$p_{2n}Ax_U^nn^c(n\mathrm{}),$$
(2.8)
for some constants $`A>0`$ and $`c`$. The exponent $`c`$ determines the singularity at $`x=x_U`$, which is, to leading order and for $`c`$, algebraic. This specializes Theorem 2.

###### Proposition 1.

Let Assumption 1 be satisfied and assume that, for $`c`$, the singular part of $`U(x)`$ is, as $`xx_U^{}`$, asymptotically given by

$$U^{(sing)}(x)U_0(x_Ux)^{c1}(xx_U^{})$$
(2.9)
for some constant $`U_00`$. Then, the PS model (2.3) has no phase transition if $`0<c<1`$ or if $`U(x_U^{})V(x_U)1`$. Otherwise, a continuous phase transition occurs for $`1<c<2`$, and a first order transition occurs if $`c>2`$, at $`w=w_c>1`$ defined in Theorem 1. ∎

Remark. If $`c=1`$ in (2.8), we get a logarithmic singularity in $`U(x)`$, hence no phase transition. If $`c=2`$ in (2.8), we get a logarithmic singularity in the derivative $`U^{}(x)`$, resulting in a continuous phase transition. See also Fisher .

### 2.3 Scaling functions

In the vicinity of a phase transition point, critical behavior of the form

$$Z(x,w)(x_cx)^\theta F((ww_c)/(x_cx)^\varphi )(x,w)(x_c^{},w_c^+),$$
(2.10)
uniformly in $`w`$, is expected with critical exponents $`\theta `$ and $`\varphi `$. Here $`F(s)`$ is a scaling function which only depends on the combined argument $`s=(ww_c)/(x_cx)^\varphi `$. The scaling function is extracted from the generating function $`Z(x,w)`$ by replacing $`w`$ using the variable $`s`$ of combined argument and expanding to leading order in $`x_cx`$. To this end, let us assume as in (2.9) that $`U(x)`$ behaves as

$$U(x)=U(x_U^{})+U_0(x_Ux)^{c1}+o\left((x_Ux)^{c1}\right)(xx_U^{}),$$
(2.11)
and $`1<c`$. In this case, the constant $`U_0`$ is negative. (If $`c<1`$, no phase transition occurs at positive temperature according to the considerations above). We have the following result.

###### Theorem 3.

Let Assumption 1 be satisfied. If the leading singularity of the loop generating function $`U(x)`$ is of the form (2.11) with an exponent $`c>1`$ and if $`U(x_U^{})V(x_U)1`$, the PS model (2.3) has critical exponents and scaling function (2.10)

$$\begin{array}{cc}\hfill \theta =\varphi =c1,& F(s)=\frac{1}{|U_0|U^2(x_U^{})V^{}(w_cx_U)x_Us}(1<c<2),\hfill \\ \hfill \theta =\varphi =1,& F(s)=\frac{1}{U^2(x_U^{})V^{}(w_cx_U)[w_cx_Us]}(c>2),\hfill \end{array}$$
(2.12)
where $`w_c1`$ is implicitly given by $`U(x_U^{})V(w_cx_U)=1`$.∎

Remark. If $`c>2`$, the critical exponents are independent of $`c`$. In both cases, the scaling function is proportional to $`F_a(s)=1/(1as)`$ with a positive constant $`a`$.

It can be shown  that under mild assumptions such a scaling function implies a certain asymptotic behavior of the function $`Z_n(w)`$ for large $`n`$. That is to say,

$$Z_n(w)x_c^{n\theta }n^{\theta 1}h\left(\left(\frac{n}{x_c}\right)^\varphi (ww_c)\right)(n\mathrm{},ww_c^+),$$
(2.13)
uniformly in $`w`$, where $`h(x)`$ is the finite-size scaling function $`h(x)=_{k=0}^{\mathrm{}}f_kx^k/\mathrm{\Gamma }(k\varphi +\theta )`$, where $`\mathrm{\Gamma }(z)`$ denotes the Gamma function, and the coefficients $`f_k`$ appear in the Taylor expansion of the scaling function $`F(s)=_{k=0}^{\mathrm{}}f_ks^k.`$ This can be used to derive the critical behavior of the fraction of shared bonds above the phase transition  for $`1<c<2`$,

$$\theta (w)=A(ww_c)^{\frac{2c}{c1}}(ww_c^+),$$
(2.14)
for some constant $`A>0`$. In our case $`F_a(s)=1/(1as)`$, we have $`h_a(x)=e^{ax}`$, if $`\theta =1`$. For $`0<\theta =c1<1`$ and $`\theta `$ rational, the finite size scaling function can be expressed in terms of hypergeometric functions. For the case $`c=3/2=\theta +1`$, which will be relevant below, we get

$$h_a(x)=\frac{1}{\sqrt{\pi }}+axe^{(ax)^2}\left(1+\text{erf}(ax)\right),$$
(2.15)
where $`\text{erf}(x)=2/\sqrt{\pi }_0^xe^{t^2}dt`$ denotes the error function.

### 2.4 Loop statistics

Within a given PS model, consider the set of all chains of length $`n`$ with $`m`$ contacts. Denote the number of loops of length $`2l`$ in this set by $`g_{m,n,l}`$. Define the generating function

$$L(x,w,y)=\underset{m,n,l}{}g_{m,n,l}w^mx^ny^l.$$
(2.16)
Again, due to the directedness of the model, this generating function can be expressed in terms of the loop generating function $`U=U(x)`$, $`\overline{U}=U(xy)`$ and the path generating function $`V=V(xw)`$. An argument as in Section 2.1 yields

$$\begin{array}{cc}\hfill L(x,w,y)& =V\overline{U}V+(V\overline{U}VUV+VUV\overline{U}V)+\mathrm{}\hfill \\ & =V\overline{U}V(1+2UV+3(UV)^2+\mathrm{})\hfill \\ & =V\overline{U}V\frac{1}{(1UV)^2}=U(yx)Z^2(x,w).\hfill \end{array}$$
(2.17)
The generating function for the total number of loops is then given by setting $`y=1`$,

$$T(x,w)=L(x,w,1)=U(x)Z^2(x,w),$$
(2.18)
and the generating function for the sum of all loop lengths is given by

$$S(x,w)=y\frac{d}{dy}L(x,w,y)|_{y=1}=xU^{}(x)Z^2(x,w).$$
(2.19)
The finite size behavior of these quantities about the phase transition can be computed from the scaling function results of Section 2.3. For the total number of loops, we get

$$\underset{ww_c^+}{lim}[x^n]T(x,w)Ax_U^{n2\theta }n^{2\theta 1}(n\mathrm{}),$$
(2.20)
where $`A>0`$ is some amplitude555 In order to simplify notation, we will denote all following amplitudes by the letter $`A`$, with the convention that their values may be different., and for the number of loops of length $`2l`$, we get

$$\begin{array}{cc}\hfill \underset{ww_c^+}{lim}[x^ny^l]L(x,w,y)& p_{2l}x_U^{ln2\theta }n^{2\theta 1}(n\mathrm{})\hfill \\ & Al^cx_U^{n2\theta }n^{2\theta 1}(1ln,n\mathrm{}).\hfill \end{array}$$
(2.21)
Note that we obtained (2.21) under the assumption that $`l`$ is asymptotically large but that $`ln`$. This is due to the estimate $`[y^l]U(yx)=p_{2l}x^lp_{2l}x_U^l`$ for $`xx_U^{}`$ and $`l`$ finite, but $`p_{2l}Ax_U^ll^c`$ for $`l\mathrm{}`$, where we assumed (2.8) to be valid.

For the sum of the loop lengths, we distinguish between the cases $`1<c<2`$ and $`c>2`$, where the derivative of the loop generating function $`U^{}(x_U^{})`$ is finite or infinite at the critical point. We find

$$\underset{ww_c^+}{lim}[x^n]S(x,w)\{\begin{array}{cc}Ax_U^{nc}n^{c1}\hfill & (1<c<2)\hfill \\ Ax_U^{n2}n\hfill & (c>2)\hfill \end{array}(n\mathrm{}).$$
(2.22)
Let us compute the probability of a loop of length $`2l`$ within chains of length $`n`$. We get, in the limit $`ww_c^+`$, for the loop length distribution

$$P(l,n)=\underset{ww_c^+}{lim}\frac{[x^ny^l]L(x,w,y)}{[x^n]T(x,w)}Al^c(1ln,n\mathrm{}).$$
(2.23)
Furthermore, for the mean loop length $`l_n`$ we get, in the limit $`ww_c^+`$,

$$l_n=\underset{l=0}{\overset{n}{}}lP(l,n)=\underset{ww_c^+}{lim}\frac{[x^n]S(x,w)}{[x^n]T(x,w)}\{\begin{array}{cc}An^{2c}\hfill & (1<c<2)\hfill \\ An^0\hfill & (c>2)\hfill \end{array}(n\mathrm{}).$$
(2.24)
The behavior of the mean loop length reflects the nature of the phase transition: If $`1<c<2`$, the mean loop length diverges, such that the transition is continuous. If $`c>2`$, the mean loop length stays finite, indicating a first order transition.

## 3 Two prominent loop classes

Since the critical behavior of PS models is essentially determined by the properties of loops, PS , and later Fisher , were led to consider various loop classes (together with straight paths for the double stranded segments). Whereas PS analyzed loop classes derived from random walks, Fisher considered loop classes derived from self-avoiding walks.

### 3.1 Loops and walks

An oriented, rooted loop of length $`n`$ is a walk of length $`n1`$, whose starting point and end point are lattice nearest neighbors. As usual \[28, Sec. 3.2\], we identify such loops if they have the same shape, i.e., if they are equal up to a translation, possibly followed by a change of orientation. These objects we call unrooted, unoriented loops, or simply loops. Each loop of length $`n`$ has at most $`2n`$ corresponding walks. If the walks are self-avoiding, each loop has exactly $`2n`$ corresponding walks. The number of loops of length $`n`$ is denoted by $`p_n`$. For example, for self-avoiding loops on $`^2`$ we have $`p_4=1`$ and $`p_6=2`$. For a given class of walks, the above description defines the corresponding (unmarked) loop class. Within a chain structure, two paths are attached to each loop. Different choices for attachment positions increase the number of (marked) loop configurations to $`\stackrel{~}{p}_np_n`$. If we assume the DNA condition that the two paths attached to a loop bisect it into pieces of equal length, then the number of possible attachments of two paths to a loop of length $`n`$ is less than or equal to $`2n`$. (We distinguish the two strands of a marked loop). PS and Fisher consider classes of oriented rooted loops. Self-avoiding oriented rooted loops can be interpreted as loops with $`2n`$ possible attachment positions of paths to a loop of length $`n`$. A similar interpretation for oriented rooted random loops is not obvious. We stress that these loop classes result in chains which are not segment-avoiding, as paths will intersect the loops. Both models cannot therefore represent real (self-avoiding) DNA.

### 3.2 Oriented rooted random loops

The first simple example of loops, which are discussed by PS , are oriented rooted random loops derived from random walks. Random walks on $`^d`$ have the generating function $`V_d(x)=1/(12dx)`$. The asymptotic behavior of the number of oriented rooted random loops of length $`2n`$ is given \[28, App. A\] by

$$\stackrel{~}{p}_{2n}A(2d)^{2n}(2n)^{d/2}(n\mathrm{}).$$
(3.1)
This implies no phase transition in $`d=2`$ and a continuous phase transition in $`d=3`$, since the phase transition condition $`U(x_U^{})V(x_U)<1`$ is satisfied in $`d=3`$: We have

$$\stackrel{~}{p}_{2n}=\underset{k+l+m=n}{}\frac{(2n)!}{(k!)^2(l!)^2(m!)^2}=\frac{((2n)!)^2}{(n!)^4}{}_{3}{}^{}F_{2}^{}(n,n,n;1,n+1/2;1/4),$$
(3.2)
where $`{}_{3}{}^{}F_{2}^{}(a_1,a_2,a_3;b_1,b_2;z)`$ is a hypergeometric function. With $`x_U=1/36`$, we extracted the amplitude $`U(x_U^{})`$ numerically  using the values $`\stackrel{~}{p}_{2n}`$ where $`n30`$, getting $`U(x_U^{})=0.51638461326(7)`$. The result follows for straight paths since $`1/(1x_U)=36/35`$ and also for random walks since $`V_3(x_U)=6/5`$.

### 3.3 Oriented rooted self-avoiding loops

The PS results led to the question  whether accounting for excluded volume effects within a loop increases the loop class exponent $`c`$, which might change the order of the phase transition. This led to considering self-avoiding loops, which are loops derived from self-avoiding walks. (A self-avoiding walk \[28, Sec. 1.1\] is a random walk which never visits a vertex twice.) By definition, the self-avoiding loop class fully accounts for excluded volume interactions within a loop. In $`d=3`$, self-avoiding loops of length $`n24`$ may be knotted. Self-avoiding walks and loops are well studied objects . Fisher considered oriented rooted loops $`\stackrel{~}{p}_{2n}=4np_{2n}`$. Their loop class exponent $`c`$, $`\stackrel{~}{p}_{2n}B\mu _d^nn^c`$, is related to the mean square displacement exponent $`\nu `$ of self-avoiding walks by the hyperscaling relation \[28, Sec. 2.1\] $`c=d\nu `$, where $`\nu =1/2`$ for $`d4`$, $`\nu =0.5877(6)`$ for $`d=3`$ and $`\nu =3/4`$ for $`d=2`$. At present, there is no proof for the values of $`\nu `$ in dimensions $`d4`$. Explicitly, we have

$$c=\{\begin{array}{cc}d/2\hfill & d4\hfill \\ 1.7631(18)\hfill & d=3\hfill \\ 3/2\hfill & d=2\hfill \end{array}$$
(3.3)
For unknotted self-avoiding loops, which is the preferrable model from a biological point of view, it has been proved  that the exponential growth constant is strictly less than that of all self-avoiding loops, while the exponent (if it exists) is expected to coincide with that of all self-avoiding loops.

Fisher concluded that the above values of the loop class exponent $`c`$ imply a continuous phase transition in $`d=2`$ and $`d=3`$. We have to check the phase transition condition: In $`d=2`$, a numerical analysis of the oriented rooted SAP series data  gives $`x_U=0.1436806285(8)`$ and $`U(x_U^{})=0.6523866(2)`$. The value of the generating function for straight paths is $`1/(1x_U)<1.678`$. In $`d=3`$, a corresponding analysis with unknotted oriented rooted SAP data and SAW data \[28, App. C\] gives $`x_U=0.045578(3)`$ and $`U(x_U^{})=0.10(1)`$. The path generating function value is bounded by that of random walks $`1/(16x_U)<1.38`$. Thus, the phase transition condition is satisfied in both $`d=2`$ and in $`d=3`$ for straight paths. Note that in $`d=2`$, self-avoiding walks as paths will result in no phase transition, since $`V(x_U)>2`$ for such paths, as follows from SAW data analysis.

## 4 Self-Avoiding Poland-Scheraga models

As discussed above, the previous PS models are not self-avoiding. Given a particular class of walks, for example, SAW, a PS model with segment-avoiding chains may be defined as follows. As paths, take only those walks with extremal first and last vertex: If $`v(0),v(1),\mathrm{},v(n)`$ are the vertices of an $`n`$-step walk $`v`$, this walk is taken as a path iff $`v_x(0)<v_x(i)<v_x(n)`$ for all $`1<i<n`$. Such walks are bridges \[28, Sec. 1.2\], whose last step is in $`x`$-direction. For (unmarked) loops, take loops derived from the walks, as explained in Section 3.1. In our example these will be SAW loops. Marking of the loops may be achieved in different ways. To this end, consider for a given loop the sets $`A_l`$ ($`A_r`$) of vertices of smallest (largest) $`x`$-coordinate. We distinguish four different types of marking: complete marking (with DNA constraint), where we mark a loop at all vertex pairs from $`A_l`$ and $`A_r`$ (whenever the DNA condition is satisfied), and unique marking (with DNA constraint), where we only mark a loop at a single vertex pair, for example the bottom vertex and the top vertex in a lexicographic ordering (if they satisfy the DNA condition).

If the walk class is self-avoiding walks, this will result in self-avoiding PS models. Unique marking would then imply $`\stackrel{~}{p}_n=2p_n`$ (we distinguish the two strands of a marked loop) and hence increase the previous exponents by one. Hence, such a PS model displays a first order transition with $`c=5/2`$ in $`d=2`$ and with $`c=2.7631(18)`$ in $`d=3`$. (The phase transition condition is satisfied, as follows from the values $`U(x_U^{})`$ given in Section 3.3, and from the estimate $`V(x_U)<1/(1dx_U)1`$ for bridges with last step in $`x`$-direction.) Unique marking with DNA constraint results in $`\stackrel{~}{p}_n2p_n`$. If we assume that the exponential growth constant for marked self-avoiding loops is still given by $`\mu _d`$ (which we expect to be true, compare Section 5.2), this implies a critical exponent $`c`$ greater or equal to the model with unique marking, i.e., a first order phase transition in $`d=2`$ and $`d=3`$. For complete marking, we have $`2p_n\stackrel{~}{p}_n2np_n`$, which rules out a decrease of $`c`$ by more than one. We however expect that the number of possible markings is of order $`1`$ as $`n\mathrm{}`$, such that $`c`$ remains unchanged, and the model again displays a first order phase transition. Similar considerations apply for the case of complete marking with DNA constraint.

## 5 Directed Poland-Scheraga models

We present two classes of directed PS models which are exactly solvable and, by definition, take into account excluded volume interactions between all parts of the structure. The first class is derived from fully directed walks and is solvable in arbitrary dimension . We will give explicit expressions for the generating function in $`d=2`$ and in $`d=3`$. This extends , where the critical behavior has been derived. The model displays a first order phase transition only for $`d>5`$. The second class, derived from partially directed walks, is considered in $`d=2`$ only. Different variants of the second model display both a first order phase transition and a continuous phase transition in $`d=2`$.

### 5.1 Fully directed walks and loops

This directed PS model consists of fully directed walks for the paths in the chain. These only take steps in positive directions. The corresponding loops are staircase polygons, which consist of two fully directed walks, which do not intersect or touch, but have a common starting point and end point. Paths are attached to these points. We distinguish the two strands of a loop. This model satisfies all the assumptions discussed in Section 2.1, and also satisfies the DNA condition that the two segments of the loops are the same length.

In $`d=2`$, the generating functions for paths and marked loops are

$$V(x)=\frac{1}{12x},U(x)=12x\sqrt{14x}.$$
(5.1)
$`U(x)`$ is twice the generating function of staircase polygons. The loop class exponent is $`c=3/2`$. We have $`x_V=1/2`$, $`x_U=1/4`$ and $`w_c=1`$, i.e., a phase transition at $`T=\mathrm{}`$. (If the empty path would not be allowed, a phase transition occurred at a finite temperature.) The free energy $`f(w)`$ is given by

$$f(w)=\mathrm{log}\left(\frac{2(w1)^2}{\sqrt{1+(w1)^2}1}\right)(1w<\mathrm{}).$$
(5.2)
The fraction of shared bonds follows as

$$\theta (w)=\frac{2w}{w1}\frac{\sqrt{1+(w1)^2}1(w1)^2/2}{1+(w1)^2\sqrt{1+(w1)^2}}(1w<\mathrm{}),$$
(5.3)
which approaches zero linearly in $`w1`$. The asymptotic behavior of $`Z_n(w)`$ about $`w=1`$ is given by

$$Z_n(w)\frac{4^{n1}}{n^{1/2}}h_{\sqrt{n}/8}(w1)(n\mathrm{},w1^+),$$
(5.4)
uniformly in $`w`$, where $`h_a(x)`$ is given by (2.15).

The PS model of fully directed walks and loops is exactly solvable in arbitrary dimension . Consider fully directed walks on $`^d`$. If there are $`k_i`$ steps in direction $`i`$, the number of distinct walks, starting from the origin, is given by the multinomial coefficient $`\left(\genfrac{}{}{0pt}{}{k_1+k_2+\mathrm{}+k_d}{k_1,k_2,\mathrm{},k_d}\right)`$. This results in $`V_d(x)=1/(1dx)`$. The generating function of a pair of such walks is given by

$$\stackrel{~}{Z}_d(x)=\underset{k_1,k_2,\mathrm{},k_d=0}{\overset{\mathrm{}}{}}\left(\genfrac{}{}{0pt}{}{k_1+k_2+\mathrm{}+k_d}{k_1,k_2,\mathrm{},k_d}\right)^2x^{k_1+\mathrm{}+k_d}.$$
(5.5)
The generating function $`\stackrel{~}{Z}_d(x)`$ can be interpreted as a chain, where each link consists of a staircase polygon or a double bond. Thus, it is related to the generating function $`U_d(x)`$, where staircase polygons are counted twice, by

$$\stackrel{~}{Z}_d(x)=\frac{1}{1(dx+U_d(x))}.$$
(5.6)
The functions $`\stackrel{~}{Z}_d(x)`$ satisfy Fuchsian differential equations of order $`d1`$, from which the singular behavior of $`\stackrel{~}{Z}_d(x)`$, and hence of $`U_d(x)`$ can be derived . In dimensions $`2<d<5`$, $`U_d(x)`$ can be expressed in terms of Heun functions. We get in $`d=3`$

$$\stackrel{~}{Z}_3^2(x)=\left(\frac{2}{\pi }\right)^2(19x)^1(1x)^1K(k_+)K(k_{}),$$
(5.7)
where $`K(k)`$ is the complete elliptic integral of the first kind, and

$$\begin{array}{cc}\hfill k_\pm ^2& =\frac{1}{2}\pm \frac{x_1}{4}(4x_1)^{\frac{1}{2}}\frac{1}{4}(2x_1)(1x_1)^{\frac{1}{2}}\hfill \\ \hfill x_1& =\frac{16x}{9x^210x+1}\hfill \end{array}$$
(5.8)
We have $`x_U=1/9`$. In arbitrary dimension $`d3`$, it has been shown  that the models have a critical point $`x_U=1/d^2`$ with exponent $`c=(d1)/2`$, with logarithmic corrections in $`d=3`$. Noting that $`U_d(x_U^{})=11/d1/\stackrel{~}{Z}_d(x_U^{})`$ and $`V_d(x_U)=d/(d1)`$, we conclude that the phase transition condition $`U_d(x_U^{})V_d(x_U)<1`$ is satisfied iff $`\stackrel{~}{Z}_d(x_U^{})<\mathrm{}`$. This is the case  in $`d4`$. For the corresponding PS model, the phase transition is thus first order for $`d6`$, and we have a continuous phase transition in dimensions $`2d5`$, being at finite temperature in $`d4`$ only.

### 5.2 Partially directed walks and loops

We consider a directed PS model, where paths are partially directed walks on the square lattice $`^2`$. These walks are self-avoiding and not allowed to take steps in the negative $`x`$ direction. The number of paths of length $`n`$ can be obtained  by considering the ways a walk of length $`n`$ can be obtained from a walk of length $`n1`$. This yields the path generating function

$$V(x)=\frac{1+x}{12xx^2},$$
(5.9)
whose dominant singularity is a simple pole at $`x_V=1\sqrt{2}`$. The associated loop class is column-convex polygons, which has been analyzed in , see Fig. 1 for an example.

We assume unique marking as explained in Section 4. The (unmarked) loop generating function $`U(x)`$ satisfies the algebraic equation of order four \[7, Eqn. 5.2\]

$$\begin{array}{cc}\hfill 0=& (x1)^4x^2+(x^37x^2+3x1)(x1)^3U(x)\hfill \\ & +2(2x^311x^2+10x4)(x1)^2U(x)^2\hfill \\ & +(5x^335x^2+47x21)(x1)U(x)^3\hfill \\ & +(2x^323x^2+38x18)U(x)^4.\hfill \end{array}$$
(5.10)
The dominant singularity at $`x_U=x_V^2=32\sqrt{2}`$ is a square-root singularity with an exponent $`c=3/2`$. The critical amplitude is given by $`U(x_U^{})=0.086983(1)`$. Thus $`2U(x_U^{})V(x_U)<1`$. This implies a continuous phase transition for the PS model at $`w_c=1.922817(1)`$.

However, the loop class defined above does not obey the DNA constraint, which might increase the exponent $`c`$. Using standard techniques , it is possible to derive a functional equation for the generating function of column-convex polygons counted by perimeter and, in addition, by their upper and lower walk length. The corresponding expression is, however, much more difficult to handle than (5.10), and there seems to be no practical way to extract the loop class exponent for the model with DNA constraint. On the other hand, it is possible to bound the loop class exponent from above by considering a refined (artificial) model. As we will argue below, this refined model has a loop class exponent $`c=5/2`$. For the former model, we thus expect a value of $`3/2c5/2`$, with a value of $`c=2`$ seeming plausible. This PS model would then be at the border between a first order and a continuous phase transition.

For the refined model, let us further assume that the upper and lower walk in a loop are not allowed to touch the $`x`$-axis. The loop class of the refined model can be described in terms of bar-graph polygons . These are column-convex polygons, one of whose walks is a straight line. Loops are composed of two bar-graph polygons, one of them in the half-plane $`y0`$, the other in the half-plane $`y0`$. Both polygons have the same horizontal coordinates and are constrained to have equal walk lengths. Paths of the model are partially directed walks, which are attached to the two loop vertices with vertical coordinate zero.

Let $`g_{m,n}`$ denote the number of bar-graph polygons with $`2m`$ horizontal and $`2n`$ vertical steps. The anisotropic perimeter generating function $`G(x,y)=g_{m,n}x^my^n`$ satisfies an algebraic equation of degree two ,

$$G(x,y)=xy+(y+x(1+y))G(x,y)+xG^2(x,y).$$
(5.11)
A generalization of the Lagrange inversion formula  can be used to obtain a closed formula for its series coefficients. It is

$$[x^my^n]G(x,y)=g_{m,n}=\underset{k=0}{\overset{n}{}}\frac{1}{m+k}\left(\genfrac{}{}{0pt}{}{m+k}{m}\right)\left(\genfrac{}{}{0pt}{}{m}{n1k}\right)\left(\genfrac{}{}{0pt}{}{m}{nk}\right).$$
(5.12)
The number $`p_{2n}`$ of all combinations of bar-graph polygons subject to equal walk lengths is then

$$p_{2n}=\underset{k=1}{\overset{(n1)/2}{}}g_{n2k,k}^2.$$
(5.13)
The first few numbers $`p_{2n}`$ are $`1,1,2,10,38,126,483,126,483,\mathrm{}`$ for $`n=3,4,5,\mathrm{}`$. We used standard methods of numerical series analysis  to estimate the critical point and critical exponent of $`U(x)`$. An analysis with first order differential approximants, using the coefficients $`p_{2n}`$ for $`70n80`$, yields the estimates

$$x_U=0.17157287(1),c=2.5000(1)$$
(5.14)
This is, to numerical accuracy, indistinguishable from $`x_U=32\sqrt{2}=0.171572875253\mathrm{}`$, as expected, and $`c=5/2`$. Note that the phase transition condition $`2U(x_U^{})V(x_U)<1`$ is satisfied, since the loop generating function is bounded by that of column-convex polygons, which is known to satisfy the above inequality.

## 6 Models extending the Poland-Scheraga class

The PS model with oriented rooted self-avoiding loops , reviewed in Section 3.3, displays a continuous phase transition in $`d=2`$ and $`d=3`$. This led to the question whether introducing excluded volume effects between different segments of the chain can change the nature of the transition . In Section 4, we introduced self-avoiding PS models in order to answer that question in the affirmative, staying within the PS class.

Another approach, which extends the PS class, consisted in simulations of self-avoiding walk pairs . We emphasize that these models are not PS models as defined in this paper, and that loop class exponents as defined in (2.8) bear no meaning for these models. The appropriate generalization is the exponent $`c^{}`$ of the loop length distribution $`P(l,n)`$ of the chain, which is assumed to behave in the limit $`ww_c^+`$ like

$$P(l,n)l^c^{}g(l/n)(1ln,n\mathrm{}),$$
(6.1)
where $`g(x)`$ is a scaling function, assumed to be constant in some analyses . Then, the same conclusions as in Section 2.2 about the nature of the phase transition determined by the loop length distribution exponent $`c^{}`$ hold. This follows from the behavior of the mean loop length $`l_n`$ in the limit $`ww_c^+`$, which is obtained from (6.1) by integration,

$$l_n=\underset{l=0}{\overset{n}{}}lP(l,n)An^{2c^{}}(n\mathrm{}).$$
(6.2)
If $`1<c^{}<2`$, the mean loop length diverges, indicating a continuous transition. If $`c^{}>2`$, the mean loop length stays finite, indicating a first order transition . Note that assumption (6.1) is not fully consistent with PS models: The behavior (2.23) is consistent with (6.1), which justifies the use of the same name for the two exponents. We found (6.2) to be satisfied for PS models if $`1<c<2`$ in Section 2.4. If $`c>2`$, however, the mean loop length exponent is independent of $`c`$ for PS models—see (2.24). This implies that a scaling form like (6.1) cannot hold for PS models with loop class exponent $`c>2`$.

The simulation results are $`c^{}=2.44(6)`$ in $`d=2`$  and $`c^{}=2.14(4)`$ in $`d=3`$ , indicating a first order phase transition in $`d=2`$ and in $`d=3`$. This led to the conclusion that introducing excluded volume effects between different segments to the PS model of rooted self-avoiding loops drives the transition from continuous to first order . The perturbative analytic arguments using the theory of polymer networks  yield a very good approximation to the value of the loop length distribution exponent found in simulations. This suggests that the main mechanism responsible for the change of critical behavior in the PS model of oriented rooted self-avoiding loops is due to “local” excluded volume effects arising from forbidden attachments of paths to loops.

## 7 Conclusion

We discussed PS models with emphasis on the order of their phase transitions. We re-analyzed the old model of oriented rooted self-avoiding loops  and found that the resulting PS model is not self-avoiding. Hence the conclusions about the order of the phase transition, which are most relevant to the recent discussions , rely on a model where most chain configurations are not self-avoiding and thus cannot represent real DNA. Our self-avoiding PS model (unique marking, with self-avoiding bridges and (unrooted) self-avoiding loops as defined in Section 4) yields a first order phase transition in both $`d=2`$ and $`d=3`$. It is expected that the other variants considered in Section 4 also yield a first order phase transition, but a detailed (numerical) analysis is needed to answer this question, as in the discussion in Section 4.

For the presumably more realistic model of pairs of interacting self-avoiding walks , our results suggest an interpretation of excluded volume effects, which complements the common one . The self-avoiding PS models of Section 4 correctly account for excluded volume effects within a loop, but overestimate excluded volume effects between different segments of a chain, due to their directed chain structure. Since this leads to a first order phase transition in $`d=2`$ and $`d=3`$, we conclude that the relaxation of excluded volume effects between different segments of the chain does not change the nature of the transition.

We also discussed several directed examples. One of these models displays a first order phase transition in $`d=2`$. Despite being exactly solvable, these models seem to be of limited relevance to the biological problem due to their directed structure and other limitations.

The key question, as to which effects are responsible for the observed behavior in real DNA, is only partially answered by these results. With respect to directed models, which also account for different base pair sequences, it seems surprising however that simulations of melting curves  agree so well with experimental curves. In these simulations, a heuristic critical exponent accounting for the statistical weight of internal loops has to be inserted. It has recently been suggested  that a value more realistic than the commonly used Fisher loop class exponent  $`c1.763`$ might be the loop length distribution exponent of the polymer network approximation  $`c^{}2.115`$. This suggestion implicitly assumes that the loop length distribution exponent of a non-directed model can be interpreted as an “effective” loop class exponent within a directed model. Although appearing plausible, it seems difficult to give this assumption a quantitative meaning. On the other hand, the MELTSIM simulation approach  is robust against a change of parameters. Indeed, melting curves can be simulated with satisfactory coincidence for both exponents, if the cooperativity parameter for the loops is adjusted accordingly . It would certainly be illuminating to obtain such simulated melting curves from a self-avoiding walk pair modeling . We also mention the recent debate about the relevance of the self-avoiding walk pair model to real DNA .

The stiffness of the double stranded segments does not seem to qualitatively alter the critical properties of chains , as can be seen by comparing the PS model result of Section 2.1, showing that critical properties are largely independent of the path generating function. The corresponding question of introducing energy costs for bending and torsion to single stranded segments has, to our knowledge, not been discussed. Even within the PS class, it is not clear how this affects the loop class exponents, since the changes just discussed might lead to loop classes with different exponential growth constants, but there is no argument as to how the exponents might change.

In conclusion, the question of the mechanisms applying in real DNA which are responsible for the denaturation process and which explain multistep behavior as observed in melting curves, are still far from being satisfactorily answered in our opinion.

## Acknowledgments

We would like to acknowledge useful discussions with Stu Whittington. CR would like to acknowledge financial support by the German Research Council (DFG) and like to thank the Erwin Schrödinger International Institute for Mathematical Physics for support during a stay in December 2002, where part of this work was done. AJG would like to acknowledge financial support from the Australian Research Council, and helpful discussions with Marianne Frommer. We acknowledge clarifying comments on an earlier version of the manuscript by Marco Baiesi, Enrico Carlon, Andreas Hanke and Ralf Metzler.


text:
# Pressure Balance between Thermal and Non-Thermal Plasmas in the 3C129 Cluster

## 1. Goals

3C129 and 3C129.1 are two radio galaxies in a nearby cluster of galaxies (fig. 1) which lies close to the galactic plane. We proposed Chandra X-ray observations in order to find signatures of hydrodynamic interactions between the hot intracluster gas and the radio lobes and tails and to estimate the non-thermal pressures within the radio structures by determining the external gas pressure and assuming pressure balance. It seemed likely that given the long tail of 3C129 ($``$ 0.5 Mpc projected length), that most of the source was close to the plane of the sky and that projection effects would thus be minimal. Since tailed radio galaxies are found only in clusters of galaxies, it also seemed likely that the projected distance between the cluster center and 3C129 was not significantly less than the actual distance: 3C129 was already known to lie close to the edge of the detected gas distribution (Leahy & Yin, 2000). For these reasons, we believed 3C129 was a good target for detecting the sort of depression in X-ray surface brightness caused by cavities in the hot gas coinciding with radio lobes (e.g. Cygnus A, Carilli et al. 1994; and Hydra A, McNamara et al. 2000).

The reason we wanted to use pressure balance between thermal and non-thermal plasmas is that this procedure provides an estimate of the non-thermal pressure which can be compared to that obtained from the usual synchrotron formulae assuming equipartition between the energies in relativistic particles and the magnetic field. It is then possible to evaluate the likelihood for various values of the critical synchrotron parameters. To do this confidently, we require that unmeasurable ram pressures be negligible and that by defining various areas of the radio source, we can infer the emitting volumes.

We use the redshift of 3C129.1 (z=0.0208, Spinrad, 1975) as the distance indicator of the cluster, and with Ho=65 km s-1 Mpc-1, one arcsec corresponds to 450 pc.

## 2. Cluster Gas

The X-ray emission from the cluster gas is shown in fig. 1. It is relatively smooth unlike the clumpy distribution seen for the gas around the radio galaxies Hydra A (McNamara et al. 2000) and M84 (Finoguenov & Jones, 2001). We performed a spectral deprojection analysis (Krawczynski 2002) and $`\beta `$ model fits to the radial distributions in two pie sections are shown in fig. 2. From these data we can obtain the thermal pressure as a function of distance from the cluster center which was judged to lie $`1^{}`$ SW of 3C129.1. It seems likely that the cluster has suffered a recent merger because the gas distribution is significantly elliptical (figs. 1 and 2); the cluster contains a radio galaxy with a long tail; and we find no evidence for a cooling flow (Krawczynski, 2002);

## 3. Non-thermal Pressures

The calculation of minimum non-thermal pressure involves 4 components: the assumption of equipartition between the relativistic particles and the magnetic field; the assumption that the filling factor is 1; the assumption that protons do not contribute significantly to the particle energy; and the assumption that most of the relativistic electrons have been ’counted’ when the synchrotron luminosity is integrated over some frequency range. If any of these assumptions are violated, the total non-thermal pressure will be greater than the minimum value.

While the usual invocation of equipartition leads to the classical $`\frac{4}{7}`$ power for the pressure’s dependence on the filling factor and the proton energy density, we should not forget that if we were to abandon equipartition and argue that we knew the average magnetic field strength from some other method (e.g. detection of inverse Compton emission from a known photon distribution), then the dependence is linear. Compare the basic equation for the total pressure with the expression for the minimum non-thermal pressure.

$$P=nkT+\frac{1}{3}\frac{c_{12}(1+k)L}{\varphi VB^{\frac{3}{2}}}+\frac{B^2}{8\pi }$$

For the field which minimizes the total pressure, BminP

$$P_{NT}=0.265[\frac{c_{12}(1+k)L}{\varphi V}]^{\frac{4}{7}}$$

where

$`c_{12}=1.06\times 10^{12}\frac{\nu _1^{(12\alpha )/2}\nu _2^{(12\alpha )/2}}{\nu _1^{1\alpha }\nu _2^{1\alpha }}`$

$`L=4\pi D_L^2(1+z)^{\alpha 1}k_s\nu ^\alpha 𝑑\nu `$

$`\varphi `$ is the filling factor for the emitting volume,V; k is the ratio of particle energy densities (protons to electrons), and $`k_s`$ is the amplitude of the radiation power law: $`S_\nu =k_s\nu ^\alpha `$.

It is also the case that various conditions can be chosen instead of the classical equipartition (particle energy density equals magnetic field energy density). We have chosen to use the magnetic field strength which minimizes the total pressure, BminP. Other choices which can change the value of the field by factors of up to 1.48 include assuming the field is smooth or tangled and equalizing pressures instead of minimizing the total pressure, or equalizing the energy densities instead of minimizing the total energy (see the appendix of Harris et al. 1995, for further details).

We have chosen regions in 3C129.1 and 3C129 on the basis of minimizing uncertainties in converting rectangular areas to cylindrical volumes and circular areas to spherical volumes (i.e. we chose regions for which the assumption that the depth dimension can be found from transverse dimensions is most likely valid). These regions are shown in fig. 3.

Flux densities for these regions were measured at 0.33, 1.4, 5, and 8 GHz. Spectral indices were determined for single or broken power laws from these flux densities as well as from spectral index maps produced by scaled arrays (5 and 8 GHz, Taylor et al. 2001) and from the similar uv coverage of our 0.33 and 1.4 GHz data.. Details of this analysis can be found in Krawczynski et al. (2003).

The results are shown in fig. 4 and given in table 1 and it is immediately clear that the minimum non-thermal pressures are almost always less than the thermal pressures, and that the difference for most regions is a factor of order 0.5.

## 4. Implications of Pressure Balance

Given the uncertainties in calculating various pressures, we find it remarkable that the thermal and non-thermal pressures are so close to each other for most of the regions selected. Extending the integration of the synchrotron spectrum down to 1 MHz (and thereby including electrons with Lorentz factors in the range $`\gamma `$ = 100 to 500) increases the minimum non-thermal pressures to values only a factor of two less than the thermal estimates for most regions. This factor of two is easy to obtain for example by invoking a value of the filling factor, $`\varphi =\frac{1}{3}`$. Taken at face value, the small factors between the thermal and non-thermal pressures indicate that we cannot accommodate large values of k (expected if relativistic protons are present); significant departures from equipartition (as hypothesized in a different context to facilitate explanations of excess euv emission as IC emission, Bowyer, this volume); anomalous numbers of low energy electrons (i.e. an excess over the extrapolation to low energies from the electron spectra inferred from the radio data); and/or values of the filling factor $`<<`$ 0.1.

## 5. Features and Conclusions

Although not all cavities will cause a depression in the X-ray surface brightness (Clarke et al. 1997), we interpret the absence of evidence for cavities associated with the radio structure of 3C129 to mean that the filling factor is likely to be something of order 1/3 or 1/10 (e.g. filaments of magnetic field and relativistic particles embedded within the ambient thermal plasma). For 3C129.1, the X-ray surface brightness is produced by a long integration of emissivities along the line of sight through the cluster center and the relatively small volume occupied by the radio structures would not be expected to produce an observable change in the X-ray surface brightness.

Whilst analyzing the morphology, we discovered X-ray emission coincident with the first few arcsec of the northern radio jet of 3C129 (Harris, Krawczynski, and Taylor, 2002) and we find hard emission from the region next to the nucleus of 3C129.1, which might be associated with a radio jet in that source (fig. 5).

We have also constructed hardness ratio maps (H-S)/H+S) with H = 2 to 5 keV and S = 1 to 2 keV. For both radio galaxies, there is a small region of significantly softer emission trailing off to the NW from the host galaxies (figs. 6 and 7), and in the case of 3C129, there is a harder region to the SE of the galaxy core (fig. 8). We will examine the statistical significance of these features in Krawczynski et al. (2003), but the obvious interpretation would be a stand-off bow shock heating the gas ahead of 3C129, and a short trail of cooler ISM being swept out of each galaxy. Although there is little debate about the direction of the relative velocity between 3C129 and the ambient gas because of the strong bending of both radio jets, previous evidence for relative motion for 3C129.1 has been marginal at best. The fact that both soft trails are in the same quadrant, is circumstantial evidence in favor of a large scale mass motion of the ICM instead of, or in addition to, the classical explanation for the long tail of 3C129 as being caused by a large velocity of the galaxy relative to a stationary cluster gas.

### Acknowledgments.

Our collaborators W. Lane, N. Kassim, and G. Taylor provided and/or reduced the radio data used in this analysis and A.G. Willis obtained new HI data from the DRAO which was used in the spectral deprojection. The work at the CfA was supported by NASA grant GO1-2135A and contract NAS8-39073; and at Yale (HK) by NASA grant GO 0-1169X.

References

Clarke, D. A., Harris, D. E., and Carilli, C. L. 1997, MNRAS 284, 981 \[Cyg A cavities\]

Carilli, C.L., Perley, R.A., & Harris, D.E. 1994, MNRAS 270, 173-177. \[Cyg A cavities\]

Finoguenov, A. & Jones, C. 2001 ApJ 547, L107 \[M84\]

Harris, D. E., Willis, A.G., Dewdney, P.E., & Batty, J. 1995, MNRAS 273, 785-789 \[appendix deals with different sorts of equipartition\]

Krawczynski, H. 2002, ApJ 569, L27 \[3C129 - the cluster deprojection\]

Krawczynski, H., Harris, D. E., Grossman, R., Lane, W., Kassim, N., & Willis, A.G. 2003, MNRAS (submitted) \[3C129-pressure; http://xxx.lanl.gov/abs/astro-ph/0302027\].

Lane,, W. M., Kassim, N. E., Ensslin, T. A., Harris, D. E., and Perley, R. A. 2002, AJ 123, 2985-2989

Leahy, D. A. & Yin, D. 2000 MNRAS 313, 617 \[ROSAT analysis\]

McNamara, B. R. et al. 2000, ApJ 534, L135 \[Hydra A cavities\]

Spinrad, H. 1975 ApJ 199, L1

Taylor, G. B., Govoni, F., Allen, S. A., & Fabian, A. C. 2001 MNRAS, 326, 2