text:
# Optimizing the fast Rydberg quantum gate

## Abstract

The fast phase gate scheme, in which the qubits are atoms confined in sites of an optical lattice, and gate operations are mediated by excitation of Rydberg states, was proposed by Jaksch et al. Phys. Rev. Lett. 85, 2208 (2000). A potential source of decoherence in this system derives from motional heating, which occurs if the ground and Rydberg states of the atom move in different optical lattice potentials. We propose to minimize this effect by choosing the lattice photon frequency $`\omega `$ so that the ground and Rydberg states have the same frequency-dependent polarizability $`\alpha (\omega )`$. The results are presented for the case of Rb.

Recently, a number of schemes for quantum computation with neutral atoms have been proposed qp1 ; qp2 ; qp3 ; qpn1 ; qp4 ; qp5 ; qpn2 ; 1 ; qpn3 . In those schemes, qubits are realized as internal states of neutral atoms trapped in optical lattices or magnetic microtraps. The two-qubit quantum gates are realized using controlled cold collisions qp1 ; qp2 ; qp3 ; qpn1 , controlled dipole-dipole interactions qp4 ; qp5 ; qpn2 , or by conditional excitations of atoms into the Rydberg states by a series of laser pulses 1 ; qpn3 . This approach to quantum computation has many advantages, such as scalability, possible massive parallelism, long decoherence times of the internal states of the atoms, flexibility in controlling atomic interactions, and well-developed experimental techniques. As extensive experimental studies of the feasibility of using neutral atoms for quantum computation are under way, more detailed theoretical investigation of the various schemes is needed. In this work, we investigate optimization of the quantum gate scheme realized by excitations to Rydberg states 1 . We refer to this scheme below as the Rydberg gate. The choice of this particular scheme results from its potential for fast (sub-microsecond) gate operations.

In the Rydberg gate scheme, the basic qubit is based on two ground hyperfine states of neutral atoms confined in an optical lattice. A two-qubit phase gate may be realized by conditionally exciting two atoms to low-lying Rydberg states. Different versions of the scheme have been proposed 1 . In those schemes, either one or both atoms may occupy the Rydberg state for much of the duration of the gate operation. However, an atom in a Rydberg state will, in general, move in different optical lattice potential than that experienced by the ground state. Therefore, the vibrational state of the atom in the lattice may change after the gate operation is completed, leading to decoherence due to motional heating. The optical potential for a given state depends on its ac polarizability, so we can seek to minimize this motional heating effect by the choice of a particular Rydberg state or of the lattice photon frequency, $`\omega `$. In this paper, we describe a method for accomplishing this by matching the frequency-dependent polarizabilities, $`\alpha (\omega )`$, of the atomic ground state and Rydberg state. The results are presented for case of Rb; however, the approach used here is applicable for logic gates with other alkali-metal atoms.

In this work, we calculate polarizabilities of the Rb atom in its ground state, $`\alpha _{5s}(\omega )`$, and in various Rydberg $`ns`$ states $`\alpha _{ns}(\omega )`$. We demonstrate that there exists a lattice frequency $`\omega `$ such that

$$\alpha _{5s}(\omega )\alpha _{ns}(\omega ).$$
(1)
This finding relies upon the fact Claude  that the frequency-dependent polarizabilities of alkali Rydberg states are close to that of a free electron, $`\alpha _{\mathrm{free}}(\omega )=e^2/m_\mathrm{e}\omega ^2`$, where $`e`$ is the elementary charge in statcoulombs, and $`m_\mathrm{e}`$ is the electron mass. For the optical frequencies treated in this paper, $`\alpha _{nl}(\omega )`$ for Rydberg states is comparable in magnitude to $`\alpha _{5s}(0)`$, but it has a negative sign. As a function of $`\omega `$, $`\alpha _{5s}(\omega )`$ increases as $`\omega `$ increases from zero, but changes sign when $`\omega `$ exceeds the first resonance frequency. We show that it will always be possible to find a frequency for which the Rydberg and ground-state polarizabilities are equal. An optical lattice constructed with light of that frequency will have well depths characteristic of far-off-resonance traps.

An alternative approach to polarizability matching of ground and Rydberg states is to find Rydberg states with resonant frequencies close to those of the ground state. In such cases the polarizabilities can be matched at near-resonant frequencies, resulting in a tighter trap for a given laser intensity. Such possibilities can be identified by a straightforward search of tabulated atomic energy levels, and they turn out to be rare. The two cases of interest in Rb are discussed at the end of this paper.

We begin with a description of our calculation of the polarizabilities. It is convenient for this purpose to use the system of atomic units, a.u., in which $`e,m_\mathrm{e}`$, $`4\pi ϵ_0`$ and the reduced Planck constant $`\mathrm{}`$ have the numerical value 1. Polarizability in a.u. has the dimensions of volume, and its numerical values presented here are thus measured in units of $`a_0^3`$, where $`a_00.052918`$ nm is Bohr radius. The atomic units for $`\alpha `$ can be be converted to SI units via $`\alpha /h`$ \[Hz/(V/m)2\]=2.48832$`\times 10^8\alpha `$ \[a.u.\], where the conversion coefficient is $`4\pi ϵ_0a_0^3/h`$ and Planck constant $`h`$ is factored out. The atomic unit of frequency $`\omega `$ is $`E_h/\mathrm{}4.1341\times 10^{16}`$ Hz, where $`E_h`$ is Hartree energy.

The valence contribution to the dynamic polarizability for an alkali atom in an $`ns`$ state can be calculated using the formula (in a.u.) relsd

$`\alpha _{ns}(\omega )`$ $`=`$ $`{\displaystyle \frac{1}{3}}{\displaystyle \underset{n^{}}{}}({\displaystyle \frac{(E_{n^{}p_{1/2}}E_{ns})n^{}p_{1/2}Dns^2}{(E_{n^{}p_{1/2}}E_{ns})^2\omega ^2}}`$ (2)
$`+`$ $`{\displaystyle \frac{(E_{n^{}p_{3/2}}E_{ns})n^{}p_{3/2}Dns^2}{(E_{n^{}p_{3/2}}E_{ns})^2\omega ^2}}),`$
where $`D`$ is the electric dipole operator and $`E_i`$ is the energy of the state $`i`$. In this formula, $`\omega `$ is assumed to be at least several linewidths off resonance with the corresponding transition. This condition is satisfied for the frequencies considered in this work, as the lattice has to be significantly detuned from the resonance to minimize spontaneous emission. The core contribution to the polarizability, calculated in the Hartree-Fock (HF) approximation, is found to be small (9.3 $`a_0^3`$) and is weakly dependent on $`\omega `$ in the frequency range considered here. The static value for the polarizability of Rb+ calculated in the random phase approximation (RPA) der ; jkh  is 9.1 $`a_0^3`$, close to the value of 9.0 $`a_0^3`$ obtained from by Johansson joh  from analysis of the observed term values of nonpenetrating Rydberg states. The accuracy of the core contribution is estimated to be 5 % in Ref. relsd . We use the RPA value for the core polarizability of Rb+ as a baseline, and adjust it to account for valence electron and the frequency dependence by using HF calculations. The RPA and HF values differ by only 2 %.

In calculating the valence contribution, the term with $`n=5`$ dominates the the sum over states, and the contributions from higher values of $`n`$ converge rapidly for all values of $`\omega `$ considered in this work. The $`5s5p`$ matrix elements were measured in Ref. 2  and the $`5snp`$ matrix elements with $`n=6,7,8`$ were calculated using single-double all-order method in Refs. relsd ; 4 . We have conducted the calculation of the ground ($`5s`$) state $`\alpha _{5s}(\omega )`$ in Rb, including terms with $`n=5,6,7`$, and 8 using matrix elements from Refs. 2 ; 4  and experimental energies from exp . The contributions to the ground state polarizability from states with $`n^{}=8\mathrm{}`$, including the continuum, and from terms with $`n=2,3,4`$ are estimated to be very small (0.2 $`a_0^3`$ and -0.3 $`a_0^3`$, respectively). The results for the values of $`\omega `$ near first two resonances, $`5s5p_{1/2}`$ and $`5s5p_{3/2}`$ are listed in Table 1. As expected, the value of the frequency-dependent polarizability varies very rapidly in the vicinity of the $`5s5p_{1/2}`$ and $`5s5p_{3/2}`$ resonances. The second set of values of $`\alpha _{5s}(\omega )`$ listed in column four illustrates the existence of the $`\alpha _{5s}(\omega )=0`$ point between the $`5s5p_{1/2}`$ and $`5s5p_{3/2}`$ resonances. The next resonance occurs at $`\omega =E_{5s}E_{6p_{1/2}}`$.

Molof et al. Molof  measured $`\alpha _{5s}(0)`$ to be $`47.3\pm 0.9`$Å3, equivalent to $`319\pm 6a_0^3`$. Our calculated value agrees with this result to within the experimental uncertainty. The static polarizability of the ground state of Rb, $`\alpha _{5s}(0)`$, was previously calculated using the techniques described above in Ref. der . The general behavior of the frequency-dependent polarizability of the ground state of Rb has been investigated previously in 1a . A result for one particular value of $`\omega =0.04298`$ a.u., corresponding to $`\lambda =1.06\mu `$m wavelength\_comment , has been calculated in the same work, using a model potential method. We compare our value of the ground state dynamic polarizability for this frequency with the result of Marinescu et al. 1a . Our result $`\alpha =693.5(9)a_0^3`$, listed in Table 1, is in good agreement with the result $`\alpha =711.4a_0^3`$ from Ref. 1a . It is just outside the range of uncertainty of the value $`\alpha =769\pm 61a_0^3`$ inferred by Bonin and Kadar-Kallen from an atomic deflection experiment Bonin . A subsequent publication tbs  will provide further details of the polarizability calculations, as well as calculation of the radiative widths of $`np`$ levels to quantify the behavior of the ground state polarizability near the resonances.

Second, we calculate the dynamic polarizabilities of the Rydberg states for a specific value of the lattice photon frequency $`\omega `$ and then investigate the dependence of polarizability on $`\omega `$ for several Rydberg states. The most likely realization of the Rydberg gate scheme involves two-photon transitions from the ground state to either $`ns`$ or $`nd`$ states. For clarity, we calculate the polarizabilities of the $`ns`$ states in this work. To make an estimate of the $`ns`$ Rydberg state polarizability $`\alpha _{ns}(\omega )`$ we first calculate relevant matrix elements and energies in Hartree-Fock approximation. The HF calculations are done on a non-linear grid, and the question of the numerical stability of the calculation of the properties of the Rydberg states was investigated by conducting the calculation with several different grids to insure the stability of the results. We use the resulting matrix elements to calculate frequency-dependent polarizabilities of the $`ns`$ states in HF approximation. The results for the $`ns`$ state polarizabilities with $`\omega =0.0576645`$ a.u. are listed in Table 2 in the columns labeled $`\alpha _{\mathrm{HF}}`$. This frequency corresponds to the value of the $`\alpha _{5s}(\omega )=290a_0^3`$ and to a detuning approximately 1/3 of the distance between the two resonances; $`\alpha _{\mathrm{free}}(\omega )=301a_0^3`$ at this frequency.

We observe that there is no significant change of the polarizability values after $`n8`$. The summation over $`n`$ in Eq. (2) is truncated at $`n=23`$. The contribution of the states with $`n>23`$ and continuum is evaluated by carrying out a calculation of the polarizability of the $`10s`$ state with B-spline basis set wrj . The summation over the entire basis set yields the result $`292`$ a.u. which differs from the value in Table 2 by only 1 %. To evaluate the uncertainty of the HF approximation further we repeat the calculations using HF matrix elements and experimental energies from Ref. exp  in Eq. (2). We find a substantial, from 10 % to 30 %, difference between these two approximations despite only a few percent differences between Hartree-Fock and experimental $`[E_{ns}E_{np}]`$ energies. Such a large discrepancy is explained by severe cancellations of the different terms in the sum of Eq. (2). As an illustration, we list the contributions to the polarizability of the $`15s`$ state from several dominant terms together with the values for the corresponding dipole matrix elements $`npD15s`$, energy differences $`\delta E=E_{np}E_{15s}`$, and denominators of Eq. (2) in Table 3. We find that the dominant contributions come from $`n=14`$ and $`n=15`$ terms which have different sign owing to the sign change in the energy differences. We also note that the denominator $`(\delta E)^2\omega ^2`$ is completely dominated by the $`\omega ^2`$ term, and, therefore, nearly identical for all of the contributions. We also note the cancellations of the smaller pairs of terms with $`n=13,16`$ and $`n=12,17`$.

The cancellations of dominant contributions have also been observed in the calculations of the vector transition polarizability $`\beta `$ and parity-nonconserving amplitude in Cs and Fr ( see, for example relsd ; 4  and references therein) when calculated using similar “direct-summation” method. Generally, the values obtained using sets of data of the consistent accuracy are more reliable in such cases, so we use HF data below. The accuracy of the values can be increased by estimating correlation correction contributions to the values of the dipole matrix elements. We conduct such calculation for $`\alpha _{8s}(\omega )`$ for the same frequency as data in Table 2 using experimental energies and high-accuracy all-order matrix elements for the dominant terms with $`n^{}=7,8`$. The result (-295 $`a_0^3`$) is listed in Table 2 and agrees very well with Hartree-Fock value. We note that substituting HF energy by experimental values and leaving matrix elements unchanged yields substantially different (-356 $`a_0^3`$) value, thus confirming our conclusion that using HF values for both matrix elements and energies produces more accurate values then replacing HF energies by experimental results. We note that current accuracy of $`\alpha _{ns}(\omega )`$ (estimated at 10 %) is sufficient for the purpose of present paper since the polarizability of the ground state varies vary rapidly with $`\omega `$ as illustrated in Table 1 and the variation of the value $`\omega `$ within the uncertainty of the Rydberg state polarizability is small. Consequently, the exact point where $`\alpha _{5s}(\omega )=\alpha _{ns}(\omega )`$ can be determined experimentally by detuning of the lattice frequency near the matching point. Further investigations of the Rydberg state polarizabilities are beyond the scope of this paper and will be considered elsewhere tbs .

Finally, we investigate the dependence of the Rydberg state polarizabilities on the value of $`\omega `$ in the vicinity of the $`5s5p_{1/2}`$ and $`5s5p_{3/2}`$ resonances. As noted above, the polarizability of $`ns`$ state with $`n>8`$ does not change significantly with $`n`$. The results of the calculation of valence dynamic polarizability for the $`8s`$, $`15s`$, and $`20s`$ states with $`\omega =0.0560.060`$ a.u. are illustrated in Fig. 1. We find that the dynamic polarizabilities of the Rydberg $`ns`$ states vary very weakly with $`\omega `$.

Summarizing the results above we find that the values of the ground state and Rydberg $`ns`$ state polarizabilities can be matched near the point between $`5s5p_{1/2}`$ and $`5s5p_{3/2}`$ resonances where the ground state polarizability changes sign as demonstrated in Fig. 2. The horizontal dashed line corresponds to polarizability of the Rydberg $`ns`$ level. The resonances are shown by the vertical dashed lines for clarity as we assumed $`\omega `$ to be a few linewidths from the resonances. We predict that value of $`\omega `$ where these polarizabilities match equals to $`\lambda _{vac}=790.14(2)`$ nm for $`15s`$ state. This value includes 10 % uncertainty in the value of the 15s state polarizability and 6 % uncertainty in the value of the ground state polarizability. The next point at which the polarizabilities can be matched is close to the $`5s6p_{1/2}`$ resonance which requires much larger detuning. The above discussion is valid for higher Rydberg states, which may be chosen for the gate implementation because of the longer lifetimes M , since we found Rydberg state ac polarizability to vary weakly with the principal quantum number.

Next, we consider an alternative approach to the issue of the matching polarizabilities of the ground and arbitrary Rydberg state. This discussion is not limited to $`ns`$ Rydberg states. It follows from the expression for the dynamic polarizability that if the energy difference $`E_{ns}E_{n^{}p}`$ or $`E_{nd}E_{n^{}p}`$ is accidentally close to the value of either $`5s5p_{1/2}`$ or $`5s5p_{3/2}`$ transition energy then the polarizability of the corresponding $`s`$, $`p`$, or $`d`$ state can be made large by detuning to the appropriate frequency. We have investigated the spectrum of Rb to locate energy differences which are close to the resonance $`5s5p_{1/2}`$ or $`5s5p_{3/2}`$ transition energies. The closest matches are $`6s15p`$ and $`4d11p`$ transition energies which differ from the $`\delta E=E_{5s}E_{5p_{3/2}}`$ by $`20`$ cm-1 and $`60`$ cm-1, respectively. The corresponding energy level scheme is illustrated in Fig. 3. However, in both of these cases the terms with such denominators contribute to the polarizabilities of the $`p`$ states. The required detuning in this case is smaller but only $`11p`$ and $`15p`$ levels may be used in this method. We note that there are no such Na levels for the situation where the lattice light is near Na(3p) levels. The advantage of this scheme is high value of the ac polarizabilities at the matching point; in the previous scheme the value of the ac polarizability at the matching point is relatively small, leading to the higher laser power requirement and subsequent higher scattering rate (on the order of 10 $`s^1`$). The disadvantage of using $`np`$ vs. $`ns`$ states in a Rydberg gate scheme is that excitation of transitions from the ground state requires either ultraviolet radiation or a three-photon process. Another difficulty of this approach is the necessity of the very fast gate operation times to achieve high gate fidelity because of the short, 4 $`\mu s`$ T , lifetime of the $`15p`$ level.

We note that by choosing the appropriate Rydberg level and a longer wavelength trapping laser (e.d. CO2 laser) it should often be possible to find a matching scheme. In these situations, because of the number of nearby resonances for the Rydberg state, one probably would have to consider contributions from several intermediate states.

In summary, we have found two ways to match frequency-dependent polarizabilities of the Rydberg state and the ground state for optical lattices tuned near the Rb(5p) states. First, the value of the ground state polarizability at $`\lambda _{vac}=790.032(8)`$ nm crosses zero and changes sign. Therefore, it can be matched with polarizabilities of various Rydberg levels near that point as demonstrated in Fig. 2. For the case of $`15s`$ state, polarizabilities match at $`\lambda _{vac}=790.14(2)`$ nm. Second, the value of $`\omega `$ can be detuned to nearly match the values of the $`15p6s`$ and $`11p4d`$ transition energies which are accidentally close to the $`\delta E=E(5p_{3/2})E(5s)`$ energy. Then, the dynamic polarizability of the $`15p`$ or $`11p`$ state can be made large enough to match the polarizability of the ground state near $`5p_{3/2}5s`$ resonance. The matching of the frequency-dependent polarizabilities of the atom in its ground and Rydberg states results in the matching of the optical potential seen by the atom during the gate operation and provides the optimal scheme for the Rydberg gate operation with respect to the motional decoherence.

Another way of eliminating the differential in the trapping potential between the ground and Rydberg states is to switch the trap off during the time of the gate action, and turn it back on after the action is completed. This procedure also induces heating of the atomic center of mass motion, the effect of which can be characterized as follows. Suppose that, at time $`t=0`$, an atom of mass $`M`$ is in the ground state of a trap site, where for convenience we assume the potential to be that of an isotropic oscillator of frequency $`\omega _0`$, with ground state energy $`E_0=3\mathrm{}\omega _0/2`$. If the trap is turned off suddenly, the wavefunction $`\mathrm{\Psi }\{𝐫,t\}`$ for the center of mass coordinate $`𝐫`$ evolves during time t as

$$\mathrm{\Psi }\{𝐫,t\}=\frac{\mathrm{exp}(r^2/2d_0^2(1+i\omega _0t))}{\pi ^{3/4}d_0^{3/2}(1+i\omega _0t)^{3/2}},$$
(3)
where $`d_0=\sqrt{\mathrm{}/M\omega _0}`$. In free expansion of the atomic wavepacket from the released trap, its mean kinetic energy, $`<T>=E_0/2`$ remains constant, while its mean-square radius grows as $`<r^2(t)>=3d_0^2(1+\omega _0^2t^2)/2`$. If the trap is suddenly turned back on again at time $`\tau `$, the mean energy of the wavepacket in the restored trapping potential is $`E=E_0/2+M\omega _0^2<r^2(\tau )>/2=E_0(1+\omega _0^2\tau ^2/2)`$. This corresponds to a heating of $`k_\mathrm{B}T=(\mathrm{}\omega _0)\omega _0^2\tau ^2/4`$ per cycle of trap release and restoration. Such heating can be avoided by use of a matched potential scheme of the type proposed here. For a trap frequency of a 1 MHz and a 1 $`\mu `$s gate time one obtains $`k_BT0.006(\mathrm{}\omega _0)`$. The gate operation time in principle can be optimized by appropriate choice of the Rydberg level, the applied dc field, and laser power.

We gratefully acknowledge discussions with Steven Rolston, Wendell Hill, Walter Johnson, Philippe Grangier, Mark Saffman, and Thad Walker. This work was partially supported by the Advanced Research Development Activity, the National Security Agency, and NIST Advanced Technology Program.


format: text

text:
# Division of Labor as the Result of Phase Transition

## 1 Introduction

The economic system is no doubt a many-particle system-it can be viewed as a collection of numerous interaction agents. So it is possible that methods and concepts developed in the study of strongly fluctuation systems might yield new results in this area. In fact, in the past decades the approaches from statistical physics have been applied in economics and a lot of interesting results, including empirical laws and theoretical models, have been achieved (see  and  as a review of Econophysics). Among all these studies, a great deal of researchers is on the agent-based modelling and related non-trivial self-organizing phenomena. In the economic system, the agents learn from each other, and their activities may be influenced by others’ actions. These interactions between agents may be simple and local, but they may have important consequence related with the emergence of global structure. To understand the mechanism behind these innovation phenomena, the methods and concepts in phase transitions and critical phenomena are helpful. For instance, in the study of majority and minority game, opinion formation, and computational ecosystems.

In this paper, we focus on the formation of labor division. Roughly speaking, an economic organizational pattern is said to involve the division of labor if it allocates the labor of different individuals to different activities. Hence the specialization of individuals and the number of professional activities are two sides of division of labor. It is a common functional organization observed in many complex systems, and it is a fundamental way to improve efficiency and utilization so as to get global optimization for the system. In order to investigate the mechanism behind the formation of labor division, we have constructed a simple model with many interacting agents. Every agent has only two kind of tasks, namely $`A`$ and $`B.`$ we describe the level of specialization of agent by his working-time share spent on producing $`A`$ or $`B`$. Each agent make their decisions for working-time in different tasks, and receive payoffs according to their and other agents’ choices. The agents can adapt by evaluating the performance of their strategies from past experience so as to get maximum returns. The returns for any agent is determined by its production with endogenous technical progress-through the mechanism of learning by doing, and its cooperation with other agents. Just like the Hamiltonian in statistical physics, the payoff function determines the behavior of the agent in economic system. Because of the bounded rationality and incompleteness information in the system, we have introduced a parameter, named social temperature $`T`$ (in the model, we have absorbed the $`\beta =k_BT`$ into the other parameters. Such an approach is traditionally used in Statistical Physics, for example, let $`H^{^{}}=\beta H`$, so the new $`J`$ in Hamiltonian of Ising model means $`\beta J`$ actually.), to describe the degree of randomness in decision-making. Then we assume that the system should obey the canonical ensemble distribution, that is the probability $`P\left(\stackrel{}{x}\right)`$ of a microstate $`\stackrel{}{x}`$ is proportional to its ‘Boltzmann factor’ determined by the total returns of the microstate $`\stackrel{}{x}`$:

$$P\left(\stackrel{}{x}\right)=\frac{e^{+\beta E\left(\stackrel{}{x}\right)}}{e^{+\beta E\left(\stackrel{}{x}\right)}𝑑\stackrel{}{x}}$$
(1)
Then using the Metropolis simulation method, we can get a Master equation to investigate the evolution of the system. With the continual change of system parameters, we have found a so called ‘social phase transitions’ phenomenon related to the emergence of labor division.

The model is defined in Section 2. The economies of specialization is introduced by increasing returns from learning by doing. And the economies of complementarity is described by an additional payoffs from the combination of two products. Section 3 gives the numerical results. The effects of parameters on critical point are discussed in detail. It is revealed that although the technical progress is a key factor that determine whether or not the labor division will happen, the competitive cooperation among agents has very important effects on the critical point. Our results are summarized in Section 4.

## 2 The Model

Let’s consider a system consists of $`N`$ homogeneous agents. In a given time period $`T`$, any agent has two different kind of tasks, namely $`A`$ and $`B`$. In any time unit, the fraction of working time on task $`A`$ for $`ith`$ agent is denoted by $`p^i\left(t\right)[0,1]`$. we are interested in the long-term ($`T\mathrm{}`$) evolution of the system, especially the emerging of labor division.

$`p^i`$ is a real number in $`[0,1]`$, which describes the working pattern for the $`ith`$ agent. If $`p^i`$ equals $`1`$ or $`0`$, the agent is full specialized in $`A`$ or $`B`$ and we call the system is in complete division of labor. If $`p^i`$ equals $`0.5`$, it means that the agent spends same time interval on tasks $`A`$ and $`B`$, and we call it as the time-dividing working mode in the following discussion. Usually $`p^i`$ could be any real number between $`0`$ and $`1`$, in this case the agent does not have any preference on job $`A`$ or $`B.`$We focus on the global behavior of the system. Based on the above description of the agent, we introduce the following two order parameters to describe the behavior of the system on macroscopic level:

$$\lambda _1\left(t\right)=\frac{1}{N}\underset{i}{}\left|2p^i\left(t\right)1\right|$$
(2)
$$\lambda _2\left(t\right)=\frac{1}{N}\underset{i}{}\left(2p^i\left(t\right)1\right)$$
(3)
Where $`\lambda _1`$ describes the intensity of labor division and cooperation in the system. It has three special values, $`\{0,0.5,1\}`$, represent time-dividing working mode, every $`p^i`$ equals $`0.5`$, no-preference working mode, $`p^i`$ distributed in (0,1) randomly, and full specialization, every $`p^i`$ equals $`0`$ or $`1`$, respectively. $`\lambda _2`$ gives the agents allocation on tasks $`A`$ and $`B`$ on global (in average on macroscopic level).

In the following discussion, we try to specify a suitable evolution rule for $`p^i\left(t\right)`$ based on the similar approach in statistical physics. So we can then determine the dynamical equations for $`p^i\left(t\right)`$, and get the final steady states for $`\lambda _1\left(t\right)`$ and $`\lambda _2\left(t\right)`$.

Analogous with Ising model in statistical physics, $`p^i`$ could be treated as the spin in $`ith`$ point of the lattice, and then $`\lambda _2`$ is the magnetization. So if we can give the Hamiltonian of the system, the evolution of the system could be determined. In the economic system, the agent switches his behavior according to his evaluations on the returns from production. Every agent try to get maximum return just as any practical tends to stay in the state with lowest energy. So the payoff function from the agent’s production should have the same effect as the Hamiltonian. The working mode for each agent should be determined by its returns.

### 2.1 Production function

According to the previous literatures on specialization and economic organization, we know that there are several factors that related to the division of labor, such as the increasing returns to specialization and transaction costs. In fact, limited ability on learning and incomplete information on the technology will lead to the labor division directly. But we don’t take them into account in this model. In our model, every agent know the technology for producing $`A`$ and $`B`$. And for any agent there is not any comparative advantages for producing $`A`$ or $`B`$ in the initial. The technical progress is achieved by the mechanism of learning by doing without any cost. The payoff function for agent $`i`$ is given by the following formulas:

$$\epsilon ^i=\epsilon _{tech}^i+\epsilon _{self}^i+\epsilon _{co}^i$$
(4)
where

$$\epsilon _{tech}^i=p^i\times \mathrm{\Gamma }_A^i+\alpha \left(1p^i\right)\times \mathrm{\Gamma }_B^i$$
(5)
$$\epsilon _{self}^i=\beta \mathrm{min}(p^i\times \mathrm{\Gamma }_A^i,\left(1p^i\right)\times \mathrm{\Gamma }_B^i)$$
(6)
$$\epsilon _{co}^i=\frac{1}{2}\beta \left(\mathrm{\Gamma }_B^i\left(1p^i\right)\mathrm{\Gamma }_A^ip^i\right)\underset{ji}{}\left(\mathrm{\Gamma }_A^jp^j\mathrm{\Gamma }_B^j\left(1p^j\right)\right)$$
(7)
In the functions above, $`\mathrm{\Gamma }_A^i`$ and $`\mathrm{\Gamma }_B^i`$ are the technology for producing $`A`$ and $`B`$ respectively. For simplicity and without losing any generality, we assume that unit production $`A`$ will get unit return. And $`\alpha (1,\mathrm{})`$ is the return of unit production $`B`$. $`\epsilon _{self}^i`$ gives the additional benefits for composite $`A`$ and $`B`$ into a final product. So it is related to the less one between $`A`$ and $`B`$. If agent $`i`$ has some more single product after the composition himself, he can also get another return, named $`\epsilon _{co}^i`$, from the cooperation with other agents. But the return got from the composition with the product from other agent is factored by $`\frac{1}{2}\beta `$. Because of the incomplete information in the system, any agent couldn’t know the situation of production on global. So in order to get the corresponding returns $`\epsilon _{co}^i`$ for every agent, we used the average-field approach similar with that in statistical physics. That means, for every agent, he will match his product with all the other agents and then the result is the average. This assumption indicate that there are already a public market in economics. In the function for $`\epsilon _{co}^i`$, $`\left(\mathrm{\Gamma }_B^i\left(1p^i\right)\mathrm{\Gamma }_A^ip^i\right)`$ gives the surplus product $`B`$ for agent $`i`$, and the sum of $`\left(\mathrm{\Gamma }_A^jp^j\mathrm{\Gamma }_B^j\left(1p^j\right)\right)`$ gives the surplus product $`A`$ of the all other agents. When these two terms have different signs, their product is negative and $`\epsilon _{co}^i`$, the corresponding payoff from the cooperation, will be also negative. In this case, the surplus product of the agent is the same kind as the final surplus of all the other agents. Because of the transaction cost and the diseconomies of incomplementarity, it is rationale that the agent would get negative payoff.

### 2.2 Two descriptions for learning by doing

The way for technical progress is learning by doing. This mechanism is related to the accumulation of agent’s historical production behavior. That means the development of comparative advantages is determined by related working time. Let $`x_{AorB}^i`$ denote the integral working time on $`A`$ or $`B`$ for agent $`i`$:

$$\begin{array}{c}x_A^i=_0^t𝑑\tau p^i\left(\tau \right)\\ x_B^i=_0^t𝑑\tau (1p^i\left(\tau \right))\end{array}$$
(8)
We have introduced two mechanisms for technical progress, named $`\gamma `$-mechanism and $`\mu `$-mechanism. They are given by following functions:

$$\mathrm{\Gamma }_{A/B}^i=\gamma ^{x_{A/B}^i}$$
(9)
$$\mathrm{\Gamma }_{A/B}^i=1+\frac{\mu x_{A/B}^{i}{}_{}{}^{2}}{\nu +x_{A/B}^{i}{}_{}{}^{2}}$$
(10)
These two mechanisms will give different results. Equation(9) is an exponential function which gives unlimited growth on technology, while there is an upper limit in equation(10). Figure(1) shows three functions for $`\mathrm{\Gamma }_{A/B}^i`$ in our computer simulations.

### 2.3 The Master equation for the system evolution

In the production function discussed above, term $`\epsilon _{co}^i`$ is the return from the agents’ cooperation. It introduces interactions among all agents (the same as the interactions among spins). So if we discuss the dynamics in $`\mu `$-space, we should construct $`N`$ evolutionary equations for $`p^i`$. But they are difficult for both theoretical analysis and computer simulation. In the following discussion, we try to describe the evolution of the system in $`\mathrm{\Gamma }`$-space. That is to discuss the evolutionary behavior of the joint density function $`P(\left\{p^i\right\},t)`$ in the space sponsored by $`\{p^i\}`$. Then the total production for the system is:

$$E=\underset{i=1}{\overset{N}{}}\epsilon _i$$
(11)
Where $`\epsilon _i`$ is the production of agent $`i`$ given by eq(4). The same as the approach of critical dynamics in Ising model, the Master equation for the joint density function $`P(\left\{p^i\right\},t)`$ in $`\mathrm{\Gamma }`$ space is

$$\frac{d}{dt}P(p^1,\mathrm{},p^i,\mathrm{},p^N;t)=\underset{i}{}_0^1𝑑x^i\begin{array}{c}[\omega ^i(x^ip^i)P(p^1,\mathrm{},x^i,\mathrm{},p^N;t)\\ \omega ^i(p^ix^i)P(p^1,\mathrm{},p^i,\mathrm{},p^N;t)]\end{array}$$
(12)
where $`\omega ^i`$ is the transition probability which is determined by the Boltzman factors as the following:

$$\omega ^i\left(p^ix^i\right)=\frac{e^{\mathrm{\Delta }E\left(p^ix^i\right)}}{_0^1𝑑x^ie^{\mathrm{\Delta }E\left(p^ix^i\right)}}$$
(13)
In the above transition probabilities, $`\mathrm{\Delta }E\left(p^ix^i\right)`$ is the change of returns related to the variation of working mode for $`ith`$ agent. From eq.(11) we can get

$$\mathrm{\Delta }E\left(p^ix^i\right)=\begin{array}{c}[(\mathrm{\Gamma }_A^i\alpha \mathrm{\Gamma }_B^i+\frac{1}{2}\beta (\mathrm{\Gamma }_A^i\mathrm{\Gamma }_B^i))(x^ip^i)\\ \frac{\beta }{2}\left(\left|x^i\mathrm{\Gamma }_A^i\left(1x^i\right)\mathrm{\Gamma }_B^i\right|\left|p^i\mathrm{\Gamma }_A^i\left(1p^i\right)\mathrm{\Gamma }_B^i\right|\right)\\ +2\epsilon _{co}^i(x^i)2\epsilon _{co}^i(p^i)]\end{array}$$
(14)
With the master equation(12) and certain initial conditions, we can get evolving behaviors of $`\lambda _1`$ and $`\lambda _2`$. Then we can discuss the phenomena on labor division. It is difficult to have some theoretical results. So in the next section we will show some numerical simulation results by Monte Carlo simulation method. Based on the Master equation(12) and corresponding transition probability(13), the simulation is proceed on the following Metropolis algorithm:

1. For any given state with $`p^i`$ for agent $`i`$, a new state $`x^i`$ is randomly selected;
2. If $`\mathrm{\Delta }E\left(p^ix^i\right)>0`$, then the transition from $`p^i`$ to $`x^i`$ is proceeded;
3. If $`\mathrm{\Delta }E\left(p^ix^i\right)0`$, then a random number $`\xi (0,1)`$ is selected. If $`\xi e^{\beta \mathrm{\Delta }E}`$, the transition from $`p^i`$ to $`x^i`$ is proceeded or else the agent $`i`$ keeps the original state.
4. For another agent $`j`$, goto step 1.

For any given initial state, the system will achieve a certain steady state after some transient process. Figure(2) gives a typical evolution behavior under $`\mu `$ mechanism.

## 3 The Results of Monte Carlo Simulation

There are several parameters related to the mechanism of labor division in the model. We will show the simple or maybe trivial results of $`\gamma `$-mechanism first and then emphasize our discussion on the results of $`\mu `$-mechanism. And we let $`\alpha =1`$ in all the simulations before we discuss the effects of parameters $`\alpha `$ and $`N`$ on the system evolution in the end.

### 3.1 $`\gamma `$-mechanism

As shown in(3) because of $`\gamma `$-mechanism is described by an exponential function, the system will achieve the state of labor division for any given initial state and $`\gamma >1`$ no matter how small is $`\gamma 1`$. That means if the intensity of technical progress is big enough, the agent would definitely specialized in a kind of task. Other factors such as competitive cooperation described by parameter $`\beta `$ have no effects on the long term evolution.

### 3.2 $`\mu `$-mechanism

But for the model with $`\mu `$-mechanism, $`\beta `$ has some important effects on the system evolution. As shown in Figure(4), there are phase transitions results in labor division and $`\beta `$ has strong effects on the critical point $`\mu _c`$ for parameter $`\mu `$. $`\mu _c`$ decreases as the $`\beta `$ increasing. The $`\mu `$-mechanism for learning by doing gives a logistic growth for technology. That is much more realistic than the $`\gamma `$-mechanism. The results reveal that the competitive cooperation among agents is very important for the emergence of labor division. However, as discussed in the end of this section, when $`\mu =0`$, that is there is no mechanism of increasing returns, the system has no phase transition to labor division. This result is rationale for the model and it is consistent with the theory on specialization. The comparative advantages in production can only be introduced by a positive feedback caused by increasing returns. So the mechanism of learning by doing is a dominant factor for labor division.

In order to investigate the effects of all terms in production function in detail, we have simulated some special cases. The first case is $`\beta =0`$. Then $`\epsilon _{self}^i`$ and $`\epsilon _{co}^i`$ in equation (4) are all zero and technical progress is the only one factor that affects the final results. As shown in figure(5), there is a phase transition in the system. Under the given conditions, the critical point for $`\mu `$ ($`\mu _c`$) is around 50.

Then we have simulated the model with $`\epsilon _{self}^i`$ but no $`\epsilon _{co}^i`$. The results are shown in (6). When $`\beta `$ is greater than a critical value $`\beta _c`$, the system will stabilized in the state with $`\lambda _1=0`$. That is $`p^i=0.5`$ for any agent $`i`$. Every agent does the tasks $`A`$ and $`B`$ himself and spends the same time on both tasks. The system is in the time-dividing working mode. Only when $`\beta `$ is smaller than $`\beta _c`$, does the system have phase transition to labor division. But the critical value $`\mu _c`$ is much larger than that of in Figure (4). It is even larger than the $`\mu _c`$ in Figure (5). So the term $`\epsilon _{self}^i`$ has the effect of anti-specialization. From the simulation results in Figure (4) and just as we have indicated in the beginning of this section, $`\epsilon _{co}^i`$ is helpful to labor division. Actually the term $`\epsilon _{co}^i`$ is the benefit for agent $`i`$ get from the cooperation with all the other agents. It should has the same effect as the decreasing of transaction costs. So it is not surprising that this term will enforce specialization.

Another interesting case is the system without technical progress, that is $`\mathrm{\Gamma }_A^i`$ and $`\mathrm{\Gamma }_B^i`$ are all equal 1. In this case, the system is exactly the same as Ising model. Let $`S^i=2p^i1`$ and $`\alpha =1`$, then the total returns for the system is

$$E=\beta \underset{i,j,ij}{}S^iS^j+\frac{1}{2}\underset{i}{}\beta \left(1\left|S^i\right|\right)$$
(15)
Where $`S^i(1,1)`$. The distribution of the state for the system is determined by $`e^E`$. Comprising with Ising model, let the corresponding Hamiltonian is

$$H=\frac{1}{2}\underset{i}{}\beta \left(1\left|S^i\right|\right)+\beta \underset{ij}{}S^iS^j$$
(16)
and let $`k_\beta T=1`$, then the system could be described by Canonical Ensemble with $`e^{\frac{H}{k_\beta T}}`$. That is an anti-ferromagnetic Ising model with global interaction. The similar approach shows that this system either does’t have phase transition related to labor division. The result is shown in Figure(7).

### 3.3 The effects of parameter $`\alpha `$ and $`N`$

Parameter $`\lambda _2`$ reflects the preference for tasks $`A`$ and $`B`$. So it is expected that parameter $`\alpha `$, the ratio of returns for $`A`$ and $`B`$, will affect final results on $`\lambda _2`$. The diagram in phase space for $`\lambda _2`$ when $`\alpha =1`$ is shown in Figure (8), in which we can see $`\lambda _2`$ fluctuates around zero, and no phase transition with symmetry breaking emerges. And as shown in Figure(9), $`\alpha `$ indeed has some effects on $`\lambda _2`$. But because of the combination mechanism described by parameter $`\beta `$, $`\alpha `$ has only a little effect on $`\lambda _2`$. It also has some detailed effects on labor division described by $`\lambda _1`$, as shown in Figure (10). But $`\alpha `$ could not change the qualitative behavior of phase transition.

Here are other two phenomena that should be pay attention to, first, the model has no scale effect for limited $`N`$, as shown in Figure (11). The phase curves for different $`N`$ are almost consistent. It has only some effects on the accurate of simulation. Second, When $`\gamma <1.0`$ in $`\gamma `$-mechanism for technical progress, the productivity declines with the time evolution. Then the system could not reach the state of specialization. This result is also meaningful for real economic system.

## 4 Concluding Remarks

A distinctive feature of the organization of a human society is the division of labor. From the classical economic theory, the division of labor comes from the development of endogenous comparative advantages. So there is the intrinsic relationship between technical progress and the evolution of the division of labor. But how and why it emerges from the system consists of identical individuals? We have studied the formation of labor division by the approach of statistical physics. The results reveal that there is a phase transition with this pattern formation. Although the progress of productivity dominated the phase transition occurs or not, the competitive cooperation among the agents has important effects on the critical point. So the market formation and labor division are usually reinforced each other. All the above results give us deep understanding to the evolution of labor division.

Studying the economy as an evolving complex system, we can avoid the standard economic assumptions of equilibrium based on rational behavior by agents. And the concepts and methods developed in statistical physics is helpful to uncover fundamental principles governing the evolution of complex adaptive systems. So the approaches presented here have potential applications in a variety of economical, biological and financial problems.

## 5 Acknowledgement

This work was supported by the National Science Foundation under Grant No.10175008 and National 973 Program under Grant No.G2000077307. We also thanks to all members of Professor Yang’s Group and faculties from Department of System Science for their warm discussion.


text:
# A model for the force stretching double-stranded chain molecules

## I Introduction

In recent years, diverse mechanical manipulation techniques of single-molecule, such as optical tweezers, atomic force microscopy, and soft microneedles have been applied in probing basic mechanical, physical and chemical properties of single biological macromolecule, such as proteins, nucleic acids, and molecular motors. The basic principle is to exert force in the $`pN`$ range on the studied molecules, and then force-extension (FECs) or extension-force curves (EFCs) are recordedsmith ; kishino ; riefp ; visscher . Many useful insights about biomolecules were obtained by analyzing these curves, going from the detailed elastic propertiessmith  to complex structure

experimentsessevaz ; bockelmann ; maier ; rief ; liphardt . In experiments of unzipping of double-stranded DNA (dsDNA), FECs show gross features, which reflect simply the local GC versus AT content along the sequenceessevaz ; bockelmann . While in stretching single-stranded DNA (ssDNA), very different EFCs have been observed: the EFCs recorded by Maier et al. showed extremely smoothmaier ; in contrast, Rief et al. found that ssDNA end-to-end distance (EED) suddenly changed during narrow force range. In recent unzipping shorter RNA experiment, except EED increasing abruptly, EFCs with intermediates were observedliphardt .

On the theoretical side, compared to the great effort in understanding the unzipping of dsDNAthompso ; cocco ; lubensky ; marenduzzo , little attention has been paid to the RNA or ssDNA stretching problemmontanari ; muller225 ; zhou1 ; zhou2 ; zhang ; muller348 ; gerland . Montanari and M$`\stackrel{´}{e}`$zard discussed the EFCs of a homogeneous ssDNA sequence. Their model exhibited a second order phase transition from a collapsed globular state to an extended necklace-like phasemontanari . Zhou and Zhang pointed out the important roles of electrostatic and base-pair stacking interactions played in stretching ssDNAzhou2 ; zhang . M$`\ddot{u}`$ller et. al. considered how the random disorder sequence modifies EFCs characteristicsmuller348 . Gerland et al. explored quantitatively how nuclei acids structure determines a broad range of FECs from very rugged to completely featureless. They found that apart from the entropy elasticity of the unpaired single strandthompso , there are two additional mechanism, the compensation effect and the contribution of suboptimal secondary structures smoothing the features of FECsgerland . Though their results were discussed on the constant extension ensemble, it may be true for constant force ensemblelubensky .

There are several reasons to attract increasing attention on the force stretching RNA or ssDNA. First, studies of stretching RNA may generate new insights into the RNA folding problemtinoco , including folding pathways and the energy landscapechen00 . Unlike the simpler hairpin structure, such as dsDNA, RNA may fold into complicated branched structures, and the native state has to compete with a large number of suboptimal states. Exerting force on the molecule will drive the native state into metastable states, which might be seen from the FECs or EFCs. Second, force stretching provides a direct approach to measure RNA thermodynamic parameters, which are key in predicting RNA secondary structuretinoco ; liphardt . In addition, the force stretching RNA is also an intriguing physical problem, e.g., how RNA secondary structures determine the outcoming EFCsgerland , and how to describe the RNA properties far from equilibrium induced by external force is still an open questionlubensky . Because the structure complexity of RNA is intermediate between dsDNA and protein, the exploration may also be valuable in understanding properties of dsDNA and proteinriefp .

Although the theoretical models about force stretching RNA moleculesmontanari ; zhou1 ; zhang ; gerland  have been proposed, more realistic and detailed models are neededtinoco . We note that with a few exceptionsgerland , almost all models were restricted the study of homogeneous or completely random RNA sequences. Whereas stretching experiments showed that sequences affect EFCs dramaticallymaier ; rief ; liphardt . In very recent experimentdessinges , Dessinges et. al  found that excluded-volume interactions play an essential role as stretching ssDNA, which were assumed to be negligible before. To take into account the two factors, in this paper we will modify and extend the statistical model of double-stranded chain conformations developed by Chen and Dill to constant force stretching ensemblechen95 ; chen98 . We think that their model is suitable for our aim. First, the model retains a relatively high degree of realism, which the sequence dependence, nonlocal intrachain contacts, and the excluded volume interactions are accounted for explicitly. Second, the model is more “microscopic”, i.e., only one parameter, base-pair contact energy is required. The entropy is calculated from the number of conformations directly. In addition, because the model can completely characterize the full energy landscape of secondary structure having chain lengths up to hundreds baseschen00 , the exploration of how force changes the folding pathways or energy landscape of RNA is valuable.

The organization of the paper is as follows. We first, in Sec II, simply overview the statistical model of double-stranded chain molecules developed by Chen and Dill. In Sec. III, we extend the models about two classes of chain molecules, hairpin and RNA-like secondary structures conformations to constant force stretching cases, respectively. Their corresponding EFCs, contact probability, population distribution and re-entering phenomena, and the difference are discussed in Sec. IV. Section. V is our conclusion. Some numerical process are relegated to appendices. In the Appendix A we illustrate the calculation of EED with a special case in hairpin conformations, namely when the closed graph is of type 1. In the Appendix B, we simply describe how to get the EED distribution of “open” self-avoiding walk (OSAW) (self-avoiding walks involving no neighboring contacts) along one coordinate through extrapolating approach.

## II The double-stranded chain model

The details of the statistical model of double-stranded chain molecules can be found in Refs. chen95  and chen98 . Here we give a brief review.

### II.1 Polymer graph theory

The model is based on polymer graphs, diagrammatic representations of the self-contacts made by different chain conformations. Fig. 1 shows a hairpin conformation and corresponding polymer graph: vertices represent the chain monomers, straight line links symbolize the covalent bonds, and curved links stand for spatial contact between monomers.

A given polymer graph represents an ensemble of chain conformations that are consistent with given contacts. Conformations having contacts other than those specified by the polymer graph belong to other graphs. In general, fewer curved links exist in a graph, a more larger number of chain conformations are consist with this graph. Any two pairs of curved links in a polymer graph must bear one of three relationships: nested, unrelated and crossing linked. Graphs of the double-stranded chain conformations involve no crossing links, examples of which include the simplest hairpin structure, such as dsDNA, and mainly secondary structures among nucleic acids and antiparallel $`\beta `$-sheets in proteins. In terms of polymer graph, the partition function of a $`(N+1)`$ monomers ((N+1)-mer) chain molecule is given as a sum over all possible polymer graphs,

$`𝒬_N(T)={\displaystyle \underset{\mathrm{\Gamma }}{}}\mathrm{\Omega }(\mathrm{\Gamma })e^{E(\mathrm{\Gamma })/k_BT},`$ (1)
where $`k_B`$ is Boltzmann constant, $`T`$ is temperature, $`\mathrm{\Gamma }`$ is an index of a possible polymer graph, $`E(\mathrm{\Gamma })`$ and $`\mathrm{\Omega }(\mathrm{\Gamma })`$ are the energy and the number of conformations of the given polymer graph $`\mathrm{\Gamma }`$, respectively. To calculate $`\mathrm{\Omega }(\mathrm{\Gamma })`$, Chen and Dill developed a matrix multiplication methodchen95 ; chen98 .

### II.2 The matrix method for a given polymer graph

A complex polymer graph can be divided into a series of faces consecutively, in which face is a region in graph that is bounded by curved and straight lines and contains no other edges; see the shadowed region in Fig. 1. Faces are classified into five types: left (L), middle (M), right (R), left-right (LR) and isolated (I), according to the arrangement of the nested curved links that bound the face. The calculation of the full partition function $`\mathrm{\Omega }(\mathrm{\Gamma })`$ for a given graph $`\mathrm{\Gamma }`$ is correspondingly separated into two steps: to count all conformations for each face, as if they are isolated and independent of each other, then to assemble these conformations into $`\mathrm{\Omega }(\mathrm{\Gamma })`$. To avoid that conformations of different faces bump into each other (excluded-volume interaction), more detailed information about the conformations of faces is needed. On two dimension (2D) lattice, it is realized by exact enumeration. First, conformations of each face are classified into sixteen types according to the shape of ports (inlet and outlet) through which it is connected to other faces; see Fig. 2. Then the compatibility between neighboring faces can be checked exactly through the spatial compatibility between the outlet of one face and the inlet of next face. It is convenient to introduce two matrices: the Face Count Matrix (FCM) $`𝐒_t`$, which matrix element $`(𝐒_t)_{ij}`$ is the number of conformations having an inlet conformation of type $`i`$ ($`1i4`$) and an outlet conformation of type $`j`$ ($`1j4`$) for a type $`t`$ face, and viability matrix $`𝐘_{t_1t_2}`$, which element $`(𝐘_{t_1t_2})_{ij}`$ is $`1`$ or $`0`$ if connection between type $`i`$ outlet of a type $`t_1`$ face and type $`j`$ inlet of a type $`t_2`$ face is viable or not viable. For a hairpin graph $`\mathrm{\Gamma }`$ having $`M`$ faces, $`\mathrm{\Omega }(\mathrm{\Gamma })`$ is derived as a product of matrices:

$`\mathrm{\Omega }(\mathrm{\Gamma })=𝐔𝐒_{t_M}𝐘_{𝐭_𝐌𝐭_{𝐌\mathrm{𝟏}}}𝐒_{t_{M1}}\mathrm{}𝐒_{t_1}𝐔^t,`$ (2)
where $`𝐔=\{1,1,1,1\}`$, $`𝐔^t`$ is the transpose of $`𝐔`$.

### II.3 The partition function

In order to calculate the partition function of a whole chain using Eq. 1, the sum over all possible polymer graphs is necessary. For the double-stranded conformations, an efficient dynamic programming algorithm was developed in Ref. chen98 . The idea is to start with a short chain segment and elongate the segment by adding one monomer for each step, and then to calculate recursively the partition function of the longer segment using the result of the shorter one. The algorithms will be discussed as needed in following sections. More useful alternative expression of Eq. 1 is

$`𝒬_N(T)={\displaystyle \underset{E}{}}g_N(E)e^{E/k_BT},`$ (3)
where $`g_N(E)`$ is the density of states, or the number of conformations having energy $`E`$, which is defined as

$`g_N(E)={\displaystyle \underset{\mathrm{\Gamma }}{}}\mathrm{\Omega }(\mathrm{\Gamma })|_{E(\mathrm{\Gamma })=E}.`$ (4)
## III Force stretching double-stranded chain molecule: the hairpin and RNA secondary structure conformations

Fig. 3 depicts the situation studied in this paper: one end of double-stranded chain molecule is fixed by a pin, and a constant force $`𝐅=f𝐱_o`$ is exerted on the other end, where $`𝐱_o`$ is a given unit vector; average extension $`x`$ along $`𝐱_o`$ is recorded as function of $`f`$. To include the contribution of force stretching energy, we extend the partition function Eq. 3 to $`𝒬_N(T;f)`$ at given force $`f`$ as

$`𝒬_N(T;f)={\displaystyle \underset{E}{}}{\displaystyle \underset{\mathrm{\Delta }}{}}g_N(E;\mathrm{\Delta })e^{\beta (Ef\mathrm{\Delta })},`$ (5)
where $`\mathrm{\Delta }`$ is EED of the chain along $`𝐱_o`$, and $`g_N(E;\mathrm{\Delta })`$ is the number of conformations whose energy and EED are $`E`$ and $`\mathrm{\Delta }`$, respectively. Define an auxiliary partition function $`𝒢_N(E;f)`$ as

$`𝒢_N(E;f)={\displaystyle \underset{\mathrm{\Delta }}{}}g_N(E;\mathrm{\Delta })e^{\beta f\mathrm{\Delta }}.`$ (6)
Then instead of Eq. 5, we have

$`𝒬_N(T;f)={\displaystyle \underset{E}{}}𝒢_N(E;f)e^{\beta E}.`$ (7)
Apparently, the following relations are satisfied: $`_\mathrm{\Delta }g_N(E;\mathrm{\Delta })=g_N(E)`$ or $`𝒢_N(E;f=0)=g_N(E)`$. The measured extension $`x`$ at force $`f`$ is calculated from partition function as

$`x=k_BT{\displaystyle \frac{}{f}}\mathrm{log}𝒬_N(T;f).`$ (8)
The calculation of $`g_N(E;\mathrm{\Delta })`$ is key in force stretching problem. Due to the new parameter, EED $`\mathrm{\Delta }`$, the dynamic programming algorithmchen98  must be modified and extended carefully. In following sections, we show how to compute $`g_N(E;\mathrm{\Delta })`$ of hairpin and RNA secondary structure conformations respectively.

### III.1 Force stretching hairpin conformations

As one of the simplest elements in secondary structure, hairpin exists in a large class of biomolecules, such as RNA hairpin, peptide $`\beta `$-hairpin and dsDNA hairpin. Recent works showed that the hairpin conformations have remarkable thermodynamic and kinetic behaviorsmunoz ; chen98 . One may hope that they also have remarkable mechanical characteristic. On the other hand, more precise formula and the essence of theory for complex RNA secondary structure conformations make us to explore them stretched by force independently.

The polymer graphs of hairpin conformations are that every curved link bears a nested relationship with respect to every other curved link. The polymer graph is composed of two parts: two non-self-contacting tail chains $`(s,s_0)`$ and $`(l_0,l)`$ and one double-stranded region $`(s_0,l_0)`$, which is defined as Closed Graph (CG)chen98 ; see Fig. 1. The number of conformations for a given graph without force equals a multiplication of the number of double-stranded conformations and the number of two tails conformations. In terms of four types of the outermost faces, the graphs are classified into four types. (LR type is excluded from hairpin conformations). To sum over all polymer graphs, two matrices, the Closed Graph Count Matrix (CGCM), $`𝐆_{}^{}{}_{t}{}^{}[E,s_0,l_0]`$ and diagonal matrix $`\omega [s,s_0;l_0,l]`$ have been defined: $`(𝐆_{}^{}{}_{t}{}^{}[E,s_0,l_0])_{ij}`$ is the sum over the number of conformations for all possible $`t`$ type graphs having energy $`E`$, given that the outmost link spans from vertex $`l_0`$ to vertex $`s_0`$, and innermost and outermost links are in $`i`$ and $`j`$ types conformations; diagonal matrix element $`(\omega [s,s_0;l_0,l])_{ii}`$ is the number of conformations of two tails $`(s,s_0)`$ and $`(l_0,l)`$ that are spatially compatible with type $`i`$ conformations. State density $`g_N(E)`$ then can be written as:

$`g_N(E)=𝐔{\displaystyle \underset{s_0,l_0,t}{}}\omega [s,s_0;l_0,l]𝐆_{}^{}{}_{t}{}^{}[E,s_0,l_0]𝐔^t,`$ (9)
here $`1t4`$chen98 .

When force is nonzero, $`g_N(E)`$ is extended to $`g_N(E;\mathrm{\Delta })`$. $`\omega [s,s_0;l_0,l]`$ is also nondegenerated into $`\omega [s,s_0;l_0,l|\mathrm{\Delta }]`$. From Eq. 9, we note that force stretching work just depends on EED of the chain, which does not involve any detailed double-stranded structure directly, i.e., only the outmost link $`(s_0,l_0)`$ contributes $`\pm 1`$ or $`0`$ to EED $`\mathrm{\Delta }`$ of the chain along $`𝐱_0`$ direction; see Fig. 4. Though, in real nucleic acid, a typical distance of a hydrogen bond is about $`3`$ times than the nucleotide distance, it does not make sense for our description on a coarse-grained level. We can exactly enumerate EED of the whole chain with three steps: first, to fix one conformation of a given CG on lattice and grow two tails which are compatible with the graph type; then, to enumerate EED of the two tails and combine them with EED of the graph according to force direction; and finally, to distribute the full conformation to lattice plane by eight square symmetry transformations ($`D_4`$ group). We define a new diagonal matrix $`\mathrm{\Pi }[s,s_0;l_0,l|f]`$,

$`\mathrm{\Pi }[s,s_0;l_0,l|f]={\displaystyle \underset{\mathrm{\Delta }}{}}\omega [s,s_0;l_0,l|\mathrm{\Delta }]e^{\beta f\mathrm{\Delta }}.`$ (10)
So we obtain

$`𝒢_N(E;f)=𝐔{\displaystyle \underset{s_0,l_0,t}{}}\mathrm{\Pi }[s,s_0;l_0,l|f]𝐆_{}^{}{}_{t}{}^{}[E,s_0,l_0]𝐔^t.`$ (11)
The basic idea is very simple, however, it is very cumbersome to complete this process. Rather than to give a detailed analysis here, we illustrate the process in Appendix A by a special case.

Though the enumeration method gives quite accurate EED of the hairpin conformations, unfortunately, as applied to longer molecule chain, it is hard to give exactly analytical formula about $`n_x^t(N;x)`$ and $`n_y^t(N;y)`$, which are defined in Appendix A. A more practical approximation is proposed in next section.

### III.2 Force stretching secondary structure conformations

Secondary structure can be seen as a tree-like structure into which four basic structure elements, helix, loops, bulges and junctions compose through self-similarity arrangement. Each graph of secondary structure conformations is divided into two levels: the first lever is a combination of unrelated double-stranded regions connected with single-stranded chains; the second level is that each CG may be viewed as an independent secondary structure except that two end monomers of the region contact. We first simply review how to calculate the state density $`g(E)`$ without force. Necessary definitions are introduced.

To be different from hairpin conformations, the calculation of full state density of secondary structure is more complex. First all possible CGCMs, $`𝐆_{}^{}{}_{t}{}^{}[E,a,b]`$ (Note that LR-type faces are included) are computed, where $`(a,b)`$ is the outmost link of the CG. Any CG is composed of smaller unrelated subclosed graphs, auxiliary matrices, $`𝐊_t[E,l,a,b]`$ counting a combination of conformations of all subgraphs are introduced, where $`\mathrm{𝑐𝑦𝑐𝑙𝑒𝑙𝑒𝑛𝑔𝑡ℎ}`$ $`l`$ is total number of monomers of the single-stranded chains in $`[a,b]`$chen98 . CGCMs can be obtained by multiplication of matrix $`𝐊_t[E,l,a,b]`$ and the number of conformations of the single-stranded part. Then in order to combine conformations of unrelated graphs into whole, matrices $`𝐆_t[E,s,a]`$ having energy $`E`$ for full polymer graph $`[s,a]`$ are defined, where $`t=0,1,2`$ represent three full graph types which are classified depending on whether $`0`$, $`1`$ and $`2`$ existing links are connected to the rightmost monomer, see Fig. 6. Their element $`(𝐆_t[E,s,a])_{ij}`$ is the number of conformations for graphs in which the outmost link (of rightmost subgraph) and the innermost link (of the leftmost subgraph) are in $`j`$ and $`i`$ conformations, respectively. All matrices mentioned above are calculated by dynamics programming algorithmchen98 . The full density of states of the secondary structure conformations is written as chen98

$`g_N(E)=𝐔{\displaystyle \underset{t=0}{\overset{2}{}}}𝐆_t[E,s,l]𝐔^t.`$ (12)
When force is exerted on the chain, exact knowledge about EED is required. In principle, we can generalize the enumeration method in hairpin case to secondary structure: to enumerate all EEDs of each single-stranded chains separated by CGs, and to combine them with contribution of the CGs, $`\pm 1`$ and $`0`$ to get the full EED. Apart from several special cases, e.g., only one or two CGs allowed in polymer graph, however, the approach is not practicable for general graphs from the experience in counting EEDs of hairpin conformations. We propose an alternative method. Because that the detailed structures of the CGs are not involved into the whole EED directly, CGs can be viewed as effective covalent bonds connecting left and right parts of conformations. Hence, the secondary structure conformations are identified into OSAWs with the reduced monomers; the effective bonds bear the excluded volume interaction of CGs with other units in the chain. The EEDs are computed by enumeration and extrapolation method. Numerical results of the latter can be found in Appendix B. Considering that conformations of monomers neighboring CGs are almost “frozen” by excluded volume interaction, it may be more reasonable to include more monomers into the effective bonds such as illustrated in Fig. 5. In present paper the simplest case is discussed. Though the effective chain approach (ECA) may overestimate the number of conformations for partially considering excluded volume interaction, it is still valuable before better methods are given.

To realize the ECA, we modify the matrices $`𝐆_t[E,s,a]`$ to $`𝐆_t[E,n,s,a]`$, where new parameter $`n`$ is the number of effective chain monomers. The recursive relations in Ref. chen98  (also see Fig. 6) then are extended into following relations:

$`𝐆_0[E,n,s,a]`$ $`=`$ $`𝐆_0[E,n1,s,a1]+𝐆_1[E,n1,s,a1],`$ (13)
$`𝐆_1[E,n,s,a]`$ $`=`$ $`\delta _{n,1}{\displaystyle \underset{t=L,M,I}{}}𝐆_{}^{}{}_{t}{}^{}[E,s,a]`$ (14)
$`+`$ $`{\displaystyle \underset{0<b<a}{}}{\displaystyle \underset{E_1}{}}𝐆_0[EE_1,n1,s,b]{\displaystyle \underset{t=L,M,I}{}}𝐆_{}^{}{}_{t}{}^{}[E_1,b,a]`$
$`+`$ $`{\displaystyle \underset{0<b<a}{}}{\displaystyle \underset{E_1}{}}𝐆_1[EE_1,n1,s,b]{\displaystyle \underset{t=M,I}{}}𝐆_{}^{}{}_{t}{}^{}[E_1,b,a],`$
$`𝐆_2[E,n,s,a]`$ $`=`$ $`\delta _{n,1}{\displaystyle \underset{t=R,LR}{}}𝐆_{}^{}{}_{t}{}^{}[E,s,a]`$ (15)
$`+`$ $`{\displaystyle \underset{0<b<a}{}}{\displaystyle \underset{E1}{}}𝐆_0[EE_1,n1,s,b]{\displaystyle \underset{t=R,LR}{}}𝐆_{}^{}{}_{t}{}^{}[E_1,b,a]`$
$`+`$ $`{\displaystyle \underset{0<b<a}{}}{\displaystyle \underset{E_1}{}}𝐆_1[EE_1,n1,s,b]𝐆_{}^{}{}_{R}{}^{}[E_1,b,a].`$
Correspondingly, $`g_N(E)`$ is extended to $`𝒢_N[E;f]`$ as

$`𝒢_N(E;f)={\displaystyle \underset{n=1}{\overset{N}{}}}{\displaystyle \underset{m}{}}C_x(n;m)e^{\beta fm}{\displaystyle \underset{t}{}}𝐔𝐆_t[E,n,s,l]𝐔^𝐭,`$ (16)
where $`C_x(n;m)`$ is the number of conformations of $`n`$-step OSAWs whose final $`x`$ coordinates are $`m`$. Eq. 16 clearly separates contribution of the unrelated CGs from the single-stranded parts.

Because hairpin is one of four elements in secondary structure, ECA should be suitable to it. We also note that exact enumeration still is available when the lengths of two tails are smaller than a given value ($`27`$ in this paper), a mixture of two methods is applied in hairpin case. We believe that it is reasonable, since as CG region $`(s_o,l_o)`$ is so large, or tails are so short that tails conformations almost are “frozen” by excluded-volume interactions, enumeration method may be accurate. While the lengths of tails are relatively longer than the CG size, or the excluded interactions between CG and tails is weaker, the ECA may be available.

### III.3 Test of ECA without force

We compute the density of states for hairpin and secondary structure conformations respectively. Comparing with Chen and Dill results, we find that our approximation gives considerable consistent results; see Figs. 7. The results can be tested against exact enumeration of shorter chains, which have been showed in Figs. 9 and 14 in Ref. chen98 .

## IV Tests of the model: EFCs, population distributions and re-entering phenomena

We have described above how the statistical mechanical model of double-stranded chain molecules is extended to force stretching problem. In this section, we make use of the model to explore EFCs, contact and population distribution, and re-entering phenomena. To account for specific monomer sequences, sequences are divided into two classes: homogeneous and RNA-like chain. Each contact of the homogeneous chain contributes a sticking energy $`\epsilon `$. While RNA-like chain has a specific sequence of four types of monomers: A, U, C and G, resembling the 4 types of bases of an RNA; only A-U pair or C-G pair contributes an sticking energy $`\epsilon `$. In following, we take the lattice spacing $`b`$.

### IV.1 Extension-force curves

We first investigate how temperature $`T`$ and monomer number $`(N+1)`$ affect EFCs of homogeneous chains restricted to hairpin and secondary structure conformations, respectively; see Fig. 8. We have following results: at given temperature, EFCs of any class conformations have similar shapes for any $`N`$-value; shapes of the EFCs become smooth when temperature increases. In particular, extensions in hairpins increase sharply to full length during a narrow force range when force reaches critical value such as $`f=0.48\epsilon /b`$ at $`T=0.1\epsilon /k_B`$. It is a typical all-or-none behavior. Whereas at the same temperature, extensions in secondary structures increase at force $`0.48\epsilon /b`$, but they reach the full length slowly and smoothly until force is about $`0.75\epsilon /b`$. In addition, the EFCs of hairpins are more stable with temperature than that of secondary structures. Generally, the dramatical jumps in hairpin EFCs show that contacts in this conformations have higher cooperativity. In contrast, for secondary structure conformations, there is a more larger number of suboptimal states, by which dramatical EED changes may be absorbed; lower cooperativity is reflected from smoothness of EFCs. The same reason may also explain the different thermal stability of EFCs with temperature observed in different structure restriction.

To illustrate effects of specific monomer sequence to EFCs, we compute EFCs of biomolecule sequences P5ab, P5abc$`\mathrm{\Delta }`$A, and P5abc; see Fig. 9(a). These sequences come from force unfolding RNA experiment studied by Liphardt et al. recentlyliphardt . Thought we do not attempt to fit experimental EFCs in this paper, EFCs of secondary structure conformations derived from our model do reveal similar information with observed in experiment: at $`T=0.1\epsilon /k_B`$, unfoldings of P5ab and P5abc start at different forces; EFC of P5abc has an apparent intermediate, i.e., after a large jump without intermediates at $`f=0.37\epsilon /b`$ there follows a inflection.

We also calculate EFCs of random, alternating and designed sequences; see Fig. 9(b). The designed is composed of four bases, A, U, C and G arranged by $`A\mathrm{}ACCCCU\mathrm{}UC\mathrm{}CAAAAG\mathrm{}G`$, where the dots represent 15-mer A, U, C and G consecutively. Comparing the EFCs of the artificial sequences with curves got before, we find the EFC of the random sequence is very similar with the EFCs of biomolecular sequences in Fig. 9(a); the EFC of the alternating sequence is almost the same with the results of homogeneous chain of secondary structure; while the EFC of the designed sequence like the EFCs of homogeneous hairpins. From designed sequence arrangement, we think that the chain prefers to form two hairpins simultaneously. So if the force can unzip one of two hairpins, it is also large enough to unzip another at the same time, which can be reflected from the large jump of the extension from small to full length abruptly.

### IV.2 Monomer-monomer contact distribution and population distribution

Our model not only can predict EFCs of different sequences, but also it relates the EFCs with structures. We calculate contact distributions $`p(i,j;T,f)`$ for every possible contact pair $`(i,j)`$chen98 . The equilibrium $`p(i,j;T,f)`$ is determined by the ration of the conditional partition function $`𝒬_N(i,j;T,f)`$ for all the conformations that contain contact $`(i,j)`$ and the full partition function $`𝒬_N(T;f)`$ defined in Eq. 5, $`p(i,j;T,f)=𝒬_N(i,j;T,f)/𝒬_N(T;f)`$. From the contact distribution, the structure of molecule is derived at given force and temperature. As an illustration, density diagrams of the distributions of the designed sequence restricted to secondary structure are shown in Fig. 10. It is apparent to see that two hairpins almost present or disappear simultaneously according to the force values, smaller or larger than $`0.42\epsilon /b`$.

Although contact distribution provides the information of chain structure, it is insufficient to see globally how force unzips complex chain molecules, or what role cooperativity plays. In thermal melting RNA, energy landscape $`F(n,nn)`$ was used to provided insights into cooperativity for the melting transition, where $`n`$ and $`nn`$ are the number of native contacts and non-native contacts respectivelychen00 . For our purpose, the EED $`\mathrm{\Delta }`$ and total contacts $`q`$ of conformations are taken as “order parameters”. Instead using free energy surface, we define population distribution

$`p(\mathrm{\Delta },q;T,f)={\displaystyle \frac{1}{𝒬_N(T;f)}}𝒬_N(\mathrm{\Delta },q;T,f),`$ (17)
where $`𝒬_N(\mathrm{\Delta },q;T,f)`$ is conditional partition function for all conformations whose EEDs and the number of contacts are $`\mathrm{\Delta }`$ and $`q`$. The probability $`p(\mathrm{\Delta },q;T,f)`$ carry the same information as the free energy. As an application, the population landscapes of the designed sequence with six typical forces at the same temperature $`0.37\epsilon /k_B`$ are shown successively in Fig. 11. At a force smaller than $`0.33\epsilon /b`$, virtually all conformations populate the “native” structure, in which two hairpins coexisting. Increasing force above $`0.33\epsilon /b`$, the conformations separate into two classes: one class is that conformations are still two hairpins coexisting; the other class is that only one hairpin exists in all conformations. Increasing force further (Fig. 11(d)), three classes of conformations present: except two mentioned above, the unfolded or full stretched conformations form a new class. Finally Fig. 11(f) shows that all molecules are full unfolded as force is larger than $`0.43\epsilon /b`$.

According to Eq. 8, the derivative of extension $`x`$ with respect to force $`f`$ can be expressed as $`d<x>/df=(<x^2><x>^2)/k_BT`$. The formula is almost the same with the heat capacity $`C(T)`$, except that extension $`x`$ is replaced by energy, and force is replaced by temperature. The analogy of heat capacity and EFCs can give lots of insights. E.g., in thermal melting proceeding, because the peaks in $`C(T)`$ curves mark molecular structure transitions, the sudden increases in EFCs might be also understood as structure transitions driven by external forces, while the counterparts of the melting temperatures at peaks in $`C(T)`$ are the “melting” forces at large jumps in EFCs. Based on the same reason, conformations between two extension jumps could be called in “intermediate” states like in thermal melting case, such as predicted in Fig. 9 and observed in experimentliphardt  about P5abc. So, just like the important roles of energy played in heat capacity concepts, we believe that EEDs may play the similar roles in force stretching problem, at least in nucleic acids. The most direct proof is that the population distribution $`p(\mathrm{\Delta },q;T,f)`$ of EED and energy in conformations are closely-related, as is showed in Fig. 11.

The corresponding of the EED and energy is not needless. To show its usefulness, we explore the EED distribution of the chains of hairpin and secondary structure conformations. The EED distributions $`p(\mathrm{\Delta };T,f)`$ can be given by a sum over contacts $`q`$ of $`p(\mathrm{\Delta },q;T,f)`$,

$`p(\mathrm{\Delta };T,f)={\displaystyle \underset{q}{}}p(\mathrm{\Delta },q;T,f).`$ (18)
Fig. 12 is the distributions of the 70-mer homogeneous double-stranded chains: for the secondary structure conformations, only one maximum is observed at any force; while in the hairpin conformations, at some forces, here near the melting force $`0.46\epsilon /b`$ at temperature $`0.25\epsilon /k_B`$, two distinct populations of shorter and longer extensions present, but no conformations with other lengths in between. We name two structure transitions as “one-state” and “two-state”, respectively. The transition definition is borrowed from thermodynamics of polymers, in which the distribution is about molecular microscopic energytostesen . The one-state and two-state transitions are resemble a critical point and “first-order transitions” in thermodynamic limit. The difference of EED distributions of two classes conformations warns that the simple EFCs may cover critical physical informations. Like the importance of the width of melting peaks in determining thermal transitions type, the third-order derivative of the average extension with respect to force is essential to determine the transitions types in the force “melting” double-stranded chain molecules.

To be different from ordinary polymers, mechanical properties of biomolecules are more complex. The main reason is the effect of monomer sequencelubensky ; muller225 ; muller348 . Hence, the types of structure transitions driven by force is inevitable to be affected by sequence. The apparent example is the EED distribution of the designed sequence, which can be see explicitly from Fig. 11: at temperature $`0.37\epsilon /k_B`$, three isolated peaks in the distribution at force $`0.38\epsilon /b`$, two peaks at force $`0.36`$ and $`0.40\epsilon /b`$. In fact, EFC of the sequence at the temperature seems very trivial; the curve is only more smoother than the EFC at $`0.1\epsilon /b`$.

### IV.3 Re-entering phenomena

Finally, we explore the extension dependence on temperatures as external force fixed. The studies are related with the phenomena about re-entering transitions happened at low temperature in stretching dsDNA case marenduzzo ; orlandini . It means that if one fixes external force at a value in a finite range and decreases temperature, beginning with stretched state, chain molecule will first collapse to a globular state; while temperature is lowered further, chain will re-enter stretched state again. One can judge the existence of the phenomenon from extension-temperature curves directly, i.e., there are dips at lower temperature on the curves for some force values. Recently M$`\ddot{u}`$ller pointed out that the re-entering phenomenon also present in force stretching RNA moleculesmuller225 . It is interesting to see whether our model can predict the transitions. We check it by numerical computing for homogeneous chains of secondary structure and hairpin conformations respectively. These results are shown in Fig. 13, The numerical data demonstrate that the re-entering transitions only present in hairpin conformations, but are not been found in secondary structure.

## V Summary

In this paper, we modify and extend the statistical model of double-stranded chain molecules developed by Chen and Dill to force stretching case. Unlike several theoretical models proposed previously, specific monomers sequence and excluded volume interaction are exactly included. We use our model to investigate how EFCs depend on chain length, sequence, and structure. In addition, our model can relate EFCs with molecular structure directly. And we also investigate re-entering phenomenon in hairpin and RNA conformations. Because our aim in this paper is just to illustrate, our model is largely simplified, e.g., the chain is restricted on 2D lattice, the chain stiffness (bending) and elasticity are neglected, and no stacking interaction is involved which may be important in force stretching theoryzhang . However, our results still reveal many important and interesting physical results. In future work, more realistic requirements, such as eliminating of lattice restrictionchen98  and adding of stacking energy could be considered.

###### Acknowledgements.

We are indebted Dr. Shi-jie Chen, Tao Xu, and Prof. H.-W. Peng for many helpful discussions in the work process.

## Appendix A Calculating end-to-end distance distribution of hairpin conformations

For hairpin conformations, the conformations of tails are determined by two factors: the outermost link type and the length of CG. To account for excluded interactions, the tails have been classified into two types, stiff (s) and flex (f)chen98 . Correspondingly, the number of conformations is $`n^{t_i}(N_i)`$, where $`t_i`$ (s or f) are the type of $`i`$th (1st or 2nd) tail with length $`N_i`$, or $`N_i`$-step OSAW. In order to calculate EEDs of whole chain , we introduce additional definitions: $`n_x^{t_i}(N_i;m)`$, the number of conformations for $`i`$th OSAW of type $`t_i`$ whose final $`x`$ coordinates are $`m`$; and $`n_y^{t_i}(N_i;m)`$ is similar except $`y`$ coordinate. It is ease to find $`n^{t_i}(N_i)=_mn_x^{t_i}(N_i;m)=_mn_y^{t_i}(N_i;m)`$. Here we illustrate the calculation of diagonal matrix element $`𝚷(l,l_0;s_0,s)_{11}`$ when the outmost link is of type $`1`$.

### A.1 $`s=s_0`$, $`ll_0`$

In this case, the length of the 2nd tail is $`N_2=ll_o`$. There are six classes of possible conformations for CG-tail complexeschen98 . In Fig. 14, we show one of them. Since the length of CG also affects exclude volume interaction, we distinguish different CG sizes.

(1) $`l_os_o=3`$. In this case, all six classes of conformations are viable. Therefore,

$`\mathrm{\Pi }[s,s_0;l_0,l|f]_{11}`$ $`=`$ $`2\times 2[{\displaystyle \underset{y}{}}n_y^{f_2}(N_2;y)(\mathrm{cosh}\beta fy+\mathrm{cosh}\beta f(y+1))`$ (19)
$`+{\displaystyle \underset{x}{}}n_x^{f_2}(N_2;x)\left(\mathrm{cosh}\beta fx+\mathrm{cosh}\beta f(x+1)\right)`$
$`+{\displaystyle \underset{y}{}}n_y^{s_2}(N_2;y)\left(2\mathrm{cosh}\beta fy+\mathrm{cosh}\beta f(y+1)+\mathrm{cosh}\beta f(y1)\right)`$
$`+{\displaystyle \underset{x}{}}n_x^{s_2}(N_2;x)\left(2\mathrm{cosh}\beta fx+2\mathrm{cosh}\beta f(x+1)\right)`$
$`(1+\mathrm{cosh}\beta f+\mathrm{cosh}\beta fN_2+\mathrm{cosh}\beta f(N_2+1))],`$
where one of coefficients $`2`$ is degeneracy degree along force direction, and another is from combining exponentials to hyperbolic cosine function, the negative part is to eliminate the straight tail conformations, which are counted repeatedly.

(2) $`l_os_o>3`$. There are five classes of conformations are involved. We have

$`\mathrm{\Pi }[s,s_0;l_0,l|f]_{11}`$ $`=`$ $`2\times 2[{\displaystyle \underset{y}{}}n_y^{f_2}(N_2;y)(\mathrm{cosh}\beta fy+\mathrm{cosh}\beta f(y+1))`$ (20)
$`+{\displaystyle \underset{x}{}}n_x^{f_2}(N_2;x)\left(\mathrm{cosh}\beta fx+\mathrm{cosh}\beta f(x+1)\right)`$
$`+{\displaystyle \underset{y}{}}n_y^{s_2}(N_2;y)\left(\mathrm{cosh}\beta fy+\mathrm{cosh}\beta f(y+1)+\mathrm{cosh}\beta f(y1)\right)`$
$`+{\displaystyle \underset{x}{}}n_x^{s_2}(N_2;x)\left(2\mathrm{cosh}\beta fx+\mathrm{cosh}\beta f(x+1)\right)`$
$`(\mathrm{cosh}\beta f+\mathrm{cosh}\beta fN_2)].`$
### A.2 $`ss_0`$, $`ll_0`$

A large number of CG-tail complexes have been listed in Ref. chen98 . Since there are two tails, we define $`N_1=s_os`$ and $`N_2=ll_o`$. According the length of CG part, we then have following results.

(1) $`l_os_o=3`$. This case is more cumbersome than above cases. We divide $`\mathrm{\Pi }[s,s_0;l_0,l|f]_{11}`$ into the sum over $`\mathrm{\Pi }^{(j)}`$, $`j=1,\mathrm{},5`$ such as

$`\mathrm{\Pi }[s,s_0;l_0,l|f]_{11}=2\times 2{\displaystyle \underset{j=1}{\overset{5}{}}}\mathrm{\Pi }^{(j)},`$ (21)
where

$`\mathrm{\Pi }^{(1)}`$ $`=`$ $`{\displaystyle \underset{y_1}{}}{\displaystyle \underset{y_2}{}}n_y^{s_1}(N_1;y_1)n_y^{s_2}(N_2;y_2)\left(2\mathrm{cosh}\beta f(y_1+y_2)+\mathrm{cosh}\beta f(y_1y_2)\right)`$ (22)
$`+{\displaystyle \underset{y_1}{}}{\displaystyle \underset{x_2}{}}n_y^{s_1}(N_1;y_1)n_x^{s_2}(N_2;x_2)(2\mathrm{cosh}\beta f(y_1+x_2+1)+\mathrm{cosh}\beta f(y_1x_21)`$
$`+\mathrm{cosh}\beta f(y_1+x_2)+\mathrm{cosh}\beta f(y_1x_2))`$
$`+{\displaystyle \underset{x_1}{}}{\displaystyle \underset{y_2}{}}n_x^{s_1}(N_1;x_1)n_y^{s_2}(N_2;y_2)(2\mathrm{cosh}\beta f(x_1+y_2+1)+\mathrm{cosh}\beta f(x_1y_2+1)`$
$`+2\mathrm{cosh}\beta f(x_1+y_2)+\mathrm{cosh}\beta f(x_1y_2))`$
$`+{\displaystyle \underset{x_1}{}}{\displaystyle \underset{x_2}{}}n_x^{s_1}(N_1;x_1)n_x^{s_2}(N_2;x_2)\left(3\mathrm{cosh}\beta f(x_1+x_2+1)\right),`$
$`\mathrm{\Pi }^{(2)}`$ $`=`$ $`{\displaystyle \underset{y_1}{}}{\displaystyle \underset{y_2}{}}n_y^{s_1}(N_1;y_1)n_y^{f_2}(N_2;y_2)\left(\mathrm{cosh}\beta f(y_1+y_2)\right)`$ (23)
$`+{\displaystyle \underset{y_1}{}}{\displaystyle \underset{x_2}{}}n_y^{s_1}(N_1;y_1)n_x^{f_2}(N_2;x_2)(\mathrm{cosh}\beta f(y_1+x_2+1)+\mathrm{cosh}\beta f(y_1+x_2)`$
$`+\mathrm{cosh}\beta f(y_1x_2))`$
$`+{\displaystyle \underset{x_1}{}}{\displaystyle \underset{y_2}{}}n_x^{s_1}(N_1;x_1)n_y^{f_2}(N_2;y_2)\left(2\mathrm{cosh}\beta f(x_1+y_2+1)+\mathrm{cosh}\beta f(x_1y_2)\right)`$
$`+{\displaystyle \underset{x_1}{}}{\displaystyle \underset{x_2}{}}n_x^{s_1}(N_1;x_1)n_x^{f_2}(N_2;x_2)\left(\mathrm{cosh}\beta f(y_1+y_2)\right),`$
$`\mathrm{\Pi }^{(3)}`$ $`=`$ $`{\displaystyle \underset{y_1}{}}{\displaystyle \underset{y_2}{}}n_y^{f_1}(N_1;y_1)n_y^{s_2}(N_2;y_2)\left(\mathrm{cosh}\beta f(y_1+y_2)\right)`$ (24)
$`+{\displaystyle \underset{y_1}{}}{\displaystyle \underset{x_2}{}}n_y^{f_1}(N_1;y_1)n_x^{s_2}(N_2;x_2)\left(2\mathrm{cosh}\beta f(y_1+x_2+1)+\mathrm{cosh}\beta f(y_1x_2)\right)`$
$`+{\displaystyle \underset{x_1}{}}{\displaystyle \underset{y_2}{}}n_x^{f_1}(N_1;x_1)n_y^{s_2}(N_2;y_2)(\mathrm{cosh}\beta f(x_1+y_2+1)+\mathrm{cosh}\beta f(x_1+y_2)`$
$`+\mathrm{cosh}\beta f(x_1y_2))`$
$`+{\displaystyle \underset{x_1}{}}{\displaystyle \underset{x_2}{}}n_x^{f_1}(N_1;x_1)n_x^{s_2}(N_2;x_2)\left(\mathrm{cosh}\beta f(x_1+x_2+1)\right),`$
$`\mathrm{\Pi }^{(4)}`$ $`=`$ $`{\displaystyle \underset{y_1}{}}{\displaystyle \underset{x_2}{}}n_y^{f_1}(N_1;y_1)n_x^{f_2}(N_2;x_2)\left(\mathrm{cosh}\beta f(y_1+x_2+1)+\mathrm{cosh}\beta f(y_1x_2)\right)`$ (25)
$`+{\displaystyle \underset{x_1}{}}{\displaystyle \underset{y_2}{}}n_x^{f_1}(N_1;x_1)n_y^{f_2}(N_2;y_2)\left(\mathrm{cosh}\beta f(x_1+y_2+1)+\mathrm{cosh}\beta f(x_1y_2)\right),`$
$`\mathrm{\Pi }^{(5)}`$ $`=`$ $`[{\displaystyle \underset{y_1}{}}n_y^{s_1}(N_1;y_1)(\mathrm{cosh}\beta f(y_1+N_2+1)+\mathrm{cosh}\beta f(y_1+N_2)+\mathrm{cosh}\beta fy_1)`$ (26)
$`+{\displaystyle \underset{x_1}{}}n_x^{s_1}(N_1;x_1)\left(\mathrm{cosh}\beta f(x_1+N_2+1)+\mathrm{cosh}\beta f(x_1+1)+\mathrm{cosh}\beta fx_1\right)`$
$`+{\displaystyle \underset{y_2}{}}n_y^{s_2}(N_2;y_2)\left(\mathrm{cosh}\beta f(y_2+N_1+1)+\mathrm{cosh}\beta f(y_2+N_1)+\mathrm{cosh}\beta fy_2\right)`$
$`+{\displaystyle \underset{x_2}{}}n_x^{s_2}(N_2;x_2)\left(\mathrm{cosh}\beta f(x_2+N_1+1)+\mathrm{cosh}\beta f(x_2+1)+\mathrm{cosh}\beta fx_2\right)`$
$`+{\displaystyle \underset{y_1}{}}n_y^{f_1}(N_1;y_1)\left(\mathrm{cosh}\beta f(y_1+N_2+1)\right)+{\displaystyle \underset{x_1}{}}n_x^{f_1}(N_1;x_1)\left(\mathrm{cosh}\beta fx_1\right)`$
$`+{\displaystyle \underset{y_2}{}}n_y^{f_2}(N_2;y_2)\left(\mathrm{cosh}\beta f(y_2+N_1+1)\right)+{\displaystyle \underset{x_2}{}}n_x^{f_2}(N_1;x_2)\left(\mathrm{cosh}\beta fx_2\right)].`$
(2) $`l_os_o>3`$.

$`\mathrm{\Pi }[s,s_0;l_0,l|f]_{11}`$ $`=`$ $`2\times [{\displaystyle \underset{y_1}{}}{\displaystyle \underset{x_2}{}}n_y^{s_1}(N_1;y_1)n_x^{s_2}(N_2;x_2)(\mathrm{cosh}\beta f(y_1+x_2+1)+\mathrm{cosh}\beta f(y_1x_2+1))`$ (27)
$`+{\displaystyle \underset{x_1}{}}{\displaystyle \underset{y_2}{}}n_x^{s_1}(N_1;x_1)n_y^{s_2}(N_2;y_2)\left(\mathrm{cosh}\beta f(x_1y_2)+\mathrm{cosh}\beta f(x_1+y_2+1)\right)`$
$`+{\displaystyle \underset{y_1}{}}{\displaystyle \underset{x_2}{}}n_y^{s_1}(N_1;y_1)n_x^{f_2}(N_2;x_2)\left(\mathrm{cosh}\beta f(y_1+x_2+1)+\mathrm{cosh}\beta f(y_1x_2)\right)`$
$`+{\displaystyle \underset{x_1}{}}{\displaystyle \underset{y_2}{}}n_x^{s_1}(N_1;x_1)n_y^{f_2}(N_2;y_2)\left(\mathrm{cosh}\beta f(x_1y_2)+\mathrm{cosh}\beta f(x_1+y_2+1)\right)`$
$`+{\displaystyle \underset{y_1}{}}{\displaystyle \underset{x_2}{}}n_y^{f_1}(N_1;y_1)n_x^{s_2}(N_2;x_2)\left(\mathrm{cosh}\beta f(y_1x_2)+\mathrm{cosh}\beta f(y_1+x_2+1)\right)`$
$`+{\displaystyle \underset{x_1}{}}{\displaystyle \underset{y_2}{}}n_x^{f_1}(N_1;x_1)n_y^{s_2}(N_2;y_2)\left(\mathrm{cosh}\beta f(x_1+y_2+1)+\mathrm{cosh}\beta f(x_1y_2)\right)`$
$`+{\displaystyle \underset{x_1}{}}{\displaystyle \underset{y_2}{}}n_x^{f_1}(N_1;x_1)n_y^{f_2}(N_2;y_2)\left(\mathrm{cosh}\beta f(x_1y_2)+\mathrm{cosh}\beta f(x_1+y_2+1)\right)`$
$`+{\displaystyle \underset{y_1}{}}{\displaystyle \underset{x_2}{}}n_y^{f_1}(N_1;y_1)n_x^{f_2}(N_2;x_2)(\mathrm{cosh}\beta f(y_1+x_2+1)+\mathrm{cosh}\beta f(y_1x_2))].`$
These equations are so boring that we are not ready to list all $`\mathrm{\Pi }[s,s_o;l_o,l|f]`$ (added up to 15 cases). If one sets force vanishing, all equations will return to Chen and Dill results for hairpin casechen98 .

## Appendix B extrapolation technique for one dimension end-to-end distance distribution of N-step open self-avoiding walk

The knowledge about the EED of N-step OSAW is important for our studies, especially, when N is larger enough, the exact enumeration can not bear. Fortunately, similar issues has been discussed early, where an asymptotic formula about the distribution of SAW was fitteddomb ; fisher66 , In this paper, we follows the similar approach to find the distribution density of the OSAW.

Let $`C_N`$ be the total number of N-steps OSAWs, and let $`C_x(N;m)`$ be the number of walks which reach $`m`$ at x coordinate. Then EED density is defined as $`\rho _N^x(m)=C_x(N;m)/C_N`$. According to the density form suggested by Domb et al., $`\rho _N^x(m)`$ is written as

$`\rho _N^x(m)={\displaystyle \frac{0.32070}{<x_N^2>^{1/2}}}\mathrm{exp}({\displaystyle \frac{0.11423x^4}{<x_N^2>^2}})`$ (28)
where $`<x_N^2>`$ is one dimensional mean-square EED. The mean square $`<x_N^2>`$ can be solved by extrapolation techniquefisher61 . Here we only give our numerical result, $`<x_N^2>=0.75N^{1.45}`$. To computer $`C_x(N;m)`$, we need $`C_N`$ from the definition of $`\rho _N^x`$, whose asymptotic formula has been provided by Chen and Dillchen98 . Several examples about EED density of OSAW are shown in Fig. 15

We note that there is a dip in the value of enumeration at the origin arising from the restriction of no returns to the origin. It might be possible to provide more refined approximationchay , but in view of the general errors involved in the extrapolation process we have not considered this worth pursuingdomb . In our studies, because chain molecule would like to extend under the force, the approximation is also reasonable.