text:
# Cyclic Spin Exchange in Cuprate Spin Ladders

\[

## Abstract

We investigate the influence of a cyclic spin exchange $`J_{cyc}`$ on the one- and two-triplet excitations of an undoped two-leg $`S`$=1/2 ladder, using the density matrix renormalization group (DMRG). The dispersion of the $`S`$=0 two-triplet bound state is dramatically reduced by $`J_{cyc}`$ due to a repulsion between triplets on neighboring rungs. In (La,Ca)14Cu24O41 a consistent description of both the spin gap and the $`S`$=0 bound state requires $`J_{cyc}/J_{}0.20`$–0.27 and $`J/J_{}1.25`$–1.35. With these coupling ratios the recently developed dynamical DMRG yields an excellent description of the entire $`S`$=0 excitation spectrum observed in the optical conductivity, including the continuum contribution.

\]

The antiferromagnetic parent compounds of the high-$`T_c`$ cuprates are thought to be the best representatives of the two-dimensional $`S`$=1/2 square-lattice Heisenberg model. Understanding their magnetic properties is of utmost importance due to the intimate relationship of magnetic correlations and high-$`T_c`$ superconductivity. Recently, the question how to set up a minimal model which accounts for these magnetic properties has been readdressed. High-resolution inelastic neutron scattering experiments performed on the two-dimensional $`S`$=1/2 antiferromagnet $`\mathrm{La}_2\mathrm{CuO}_4`$  exhibit a magnon dispersion at the zone boundary which cannot be obtained within a nearest-neighbor Heisenberg model. It has been argued that the inclusion of a cyclic spin exchange term of about 20% would reproduce this dispersion . This cyclic spin exchange emerges as a correction to the nearest-neighbor Heisenberg Hamiltonian in order $`t^4/U^3`$ from a $`t/U`$-expansion of the one-band Hubbard model . It is expected to be the dominant correction within a more realistic three-band description of the $`\mathrm{CuO}_2`$-planes because there the cyclic permutation of 4 spins on a plaquette can take place without double occupancy . Similar cyclic spin exchange processes have proven to be significant in other systems, e.g. in cuprate spin chains a ferromagnetic 2-spin cyclic exchange is responsible for the unusually strong exchange anisotropy .

Cuprate spin ladders offer an alternative approach to decide about the existence and potential implications of a cyclic spin exchange term. They are composed of the same corner-sharing Cu-O plaquettes as the 2D cuprates, thus similar exchange couplings are expected for all spin products. In fact the inclusion of a cyclic spin exchange has been suggested in order to explain the smallness of the ladder spin gap observed in (La,Ca)14Cu24O41 . However, it is impossible to extract a unique set of coupling parameters or even to settle the existence of a cyclic spin exchange term from an analysis of the spin gap only. Here, the optical conductivity $`\sigma (\omega )`$  can provide the missing information. Magnetic excitations can be observed in $`\sigma (\omega )`$ via the simultaneous excitation of a phonon . In (La,Ca)14Cu24O41, two peaks in $`\sigma (\omega )`$ were identified as $`S`$=0 bound states of two triplets .

Here we show that a consistent evaluation of $`\sigma (\omega )`$ is achieved by inclusion of $`J_{cyc}`$. The dispersion of the $`S`$=0 bound state is very sensitive to the magnitude of $`J_{cyc}`$. The combined analysis of both the $`S`$=0 bound state and the spin gap establishes a unique set of exchange couplings, namely $`J_{cyc}/J_{}0.20`$-0.27, $`J/J_{}1.25`$-1.35 and $`J_{}=950`$-1100 cm-1 for (La,Ca)14Cu24O41. For these values, we calculate the entire spectrum of $`\sigma (\omega )`$ using the dynamical DMRG  and find excellent agreement with new experimental data which yield an improved estimate of the continuum contribution. This clearly confirms the above value for $`J_{cyc}`$. One might speculate that such a sizeable cyclic spin exchange term will have significant consequences for superconductivity in the doped spin ladders  since magnetic and pairing correlations are considered as closely related phenomena .

We study an isolated two-leg $`S`$=1/2 ladder characterized by an antiferromagnetic Heisenberg Hamiltonian with an additional cyclic spin exchange term:

$`H`$ $`=J_{}{\displaystyle \underset{i}{}}𝐒_{i,l}𝐒_{i,r}+J{\displaystyle \underset{i}{}}\left(𝐒_{i,l}𝐒_{i+1,l}+𝐒_{i,r}𝐒_{i+1,r}\right)`$ (1)
$`+`$ $`J_{cyc}{\displaystyle \underset{i}{}}{\displaystyle \frac{1}{4}}(P_{(i,l),(i,r),(i+1,r),(i+1,l)}+P_{(i,l)(i,r)(i+1,r)(i+1,l)}^1)`$ (2)
where $`J_{}`$ and $`J`$ denote the rung and leg couplings, the index $`i`$ refers to the rungs, and $`l`$, $`r`$ label the two legs. The cyclic permutation operator $`P_{1234}`$ for 4 spins on a plaquette is given by:

$`P`$ $`{}_{1234}{}^{}+P_{1234}^1=`$ (5)
$`4(𝐒_1𝐒_2)(𝐒_3𝐒_4)+4(𝐒_1𝐒_4)(𝐒_2𝐒_3)4(𝐒_1𝐒_3)(𝐒_2𝐒_4)`$
$`+𝐒_1𝐒_2+𝐒_3𝐒_4+𝐒_1𝐒_4+𝐒_2𝐒_3+𝐒_1𝐒_3+𝐒_2𝐒_4.`$
The influence of $`J_{cyc}`$ can easily be understood in terms of the representation of the spin operators by rung singlets $`s_i^{}`$ and triplets $`t_{i,\alpha }^{}`$ ($`\alpha =x,y,z`$) :

$`H={\displaystyle \underset{i}{}}`$ $`\{{\displaystyle \frac{9}{16}}J_{cyc}{\displaystyle \frac{3}{4}}(J_{}+{\displaystyle \frac{1}{2}}J_{cyc})s_i^{}s_i`$ (12)
$`+{\displaystyle \frac{1}{4}}(J_{}{\displaystyle \frac{11}{2}}J_{cyc}){\displaystyle _\alpha }t_{i,\alpha }^{}t_{i,\alpha }`$
$`+{\displaystyle \frac{1}{2}}(J+{\displaystyle \frac{1}{2}}J_{cyc}){\displaystyle _\alpha }[s_i^{}s_{i+1}t_{i+1,\alpha }^{}t_{i,\alpha }+\mathrm{H}.\mathrm{c}.]`$
$`+{\displaystyle \frac{1}{2}}(J{\displaystyle \frac{1}{2}}J_{cyc}){\displaystyle _\alpha }[s_i^{}s_{i+1}^{}t_{i+1,\alpha }t_{i,\alpha }+\mathrm{H}.\mathrm{c}.]`$
$`+{\displaystyle \frac{1}{4}}(J+{\displaystyle \frac{1}{2}}J_{cyc}){\displaystyle _{\alpha ,\beta }}(1\delta _{\alpha \beta })\times `$
$`[t_{i,\alpha }^{}t_{i+1,\beta }^{}t_{i+1,\alpha }t_{i,\beta }t_{i,\alpha }^{}t_{i+1,\alpha }^{}t_{i+1,\beta }t_{i,\beta }+\mathrm{H}.\mathrm{c}.]`$
$`+J_{cyc}{\displaystyle _{\alpha ,\beta }}t_{i,\alpha }^{}t_{i,\alpha }t_{i+1,\beta }^{}t_{i+1,\beta }\}.`$
The consideration of $`J_{cyc}`$ primarily renormalizes the coupling strengths of the individual terms in the original Heisenberg Hamiltonian. The only new contribution is the last term, a repulsive interaction between triplets on neighboring rungs. With respect to the renormalization of the coupling constants the strongest effect is a reduction of the triplet on-site energy. Consequently the spin gap as well as the energy of the entire one-triplet dispersion are reduced. In Fig. 1(a) we show the one-triplet dispersion, obtained by the Lanczos-vector method  for a N=80-site ladder with $`J=J_{}`$ and different values of $`J_{cyc}`$. The shape of the one-triplet dispersion remains qualitatively unchanged. The strongest suppression affects the minimum ($`p_x=\pi `$) and the maximum (near $`p_x\pi /3`$). Nevertheless, there still is a local minimum at $`p_x`$=0 even for $`J_{cyc}/J=0.3`$, in contrast to the exact diagonalization results for a N=24-site ladder in Ref. . Extrapolating to an infinite ladder , we find that the spin gap $`\mathrm{\Delta }_s=E_0(N=\mathrm{},S=1)E_0(N=\mathrm{},S=0)`$ decreases exponentially (inset of Fig. 1(a)), but we cannot conclusively decide from our present data, whether the spin gap indeed vanishes at $`J_{cyc}/J1.2`$, as has been proposed in Ref..

Similar considerations also apply to the two-triplet excitations. With increasing $`J_{cyc}`$ the lower edge of the two-triplet continuum (open symbols in Fig. 1(b)) shifts corresponding to the one-triplet dispersion. In addition, however, a repulsive interaction between neighboring triplets is introduced by $`J_{cyc}`$, as mentioned above (see Eq. (12)). This reduces the binding energy of the $`S`$=0 two-triplet bound state. As a consequence the width of the bound-state dispersion decreases dramatically (full symbols in Fig. 1(b)). Evidently the dispersion of the $`S`$=0 two-triplet bound state is very sensitive to the magnitude of $`J_{cyc}`$ and thus provides a suitable quantity to probe the existence and magnitude of a cyclic spin exchange.

In (La,Ca)14Cu24O41 the $`S`$=0 bound state has been observed in the optical conductivity $`\sigma (\omega )`$ via phonon-assisted two-triplet absorption . Although the total momentum of the excitation has to be zero, the simultaneous excitation of a phonon with arbitrary momentum permits to measure a weighted average of the $`S`$=0 excitations from the entire Brillouin zone. In Fig. 3 we show the magnetic contribution to $`\sigma (\omega )`$ in La5.2Ca8.8Cu24O41 . For polarization of the electrical field parallel to the legs, the two peaks in $`\sigma _{\mathrm{leg}}(\omega )`$ at $`\omega _1`$=2140 cm-1 and $`\omega _2`$=2780 cm-1 correspond to van Hove singularities resulting from the maximum $`\omega _{\mathrm{max}}`$ and minimum $`\omega _{\mathrm{min}}`$ of the $`S`$=0 bound-state dispersion at $`p_x\pi /2`$ and $`p_x=\pi `$, respectively (see Fig. 1(b)). The broad features at higher energies represent the continuum contribution.

In order to determine the exchange couplings for La5.2Ca8.8Cu24O41 we have calculated $`\omega _{\mathrm{max}}`$ and $`\omega _{\mathrm{min}}`$ for different coupling ratios $`J/J_{}`$ and $`J_{cyc}/J_{}`$ using the Lanczos-vector method . By extrapolation to an infinite ladder  we obtain the values displayed in Fig. 2(a). The corresponding results for the spin gap $`\mathrm{\Delta }_s`$ are shown in Fig. 2(b). In the range depicted in Fig. 2 these values may be approximated within 5% by $`\omega _{\mathrm{min}}1.64J_{}0.54J_{cyc}`$, $`\omega _{\mathrm{max}}0.61J_{}1.87J_{cyc}+1.53J`$ and $`\mathrm{\Delta }_s0.48J_{}0.84J_{cyc}`$. In principle the three experimental quantities $`\omega _1`$, $`\omega _2`$ and $`\mathrm{\Delta }_s280`$ cm-1  should suffice to determine $`J_{cyc}`$, $`J`$ and $`J_{}`$. However, $`\omega _1`$ and $`\omega _2`$ depend also on the frequency $`\omega _{ph}`$ of the phonon involved in the optical absorption process, i.e., $`\omega _1=\omega _{\mathrm{min}}+\omega _{ph}^{p_x=\pi }`$ and $`\omega _2=\omega _{\mathrm{max}}+\omega _{ph}^{p_x=\pi /2}`$. Making use of the relative insensitivity of $`\mathrm{\Delta }_s`$ and $`\omega _{\mathrm{min}}`$ on $`J/J_{}`$, we can determine the magnitude of the cyclic spin exchange without consideration of $`\omega _2`$ and hence omit any effect of a possible phonon dispersion. For a conservative estimate of the Cu-O bond-stretching phonon frequency $`\omega _{ph}^{p_x=\pi }=600\pm 100`$cm-1 we obtain $`J_{cyc}/J_{}=0.200.27`$ with $`J_{}9501100`$cm-1. The dispersion of the phonon enters only in the determination of the coupling ratio $`J/J_{}`$. Allowing for a dispersion of $`\omega _{ph}^{p_x=\pi /2}\omega _{ph}^{p_x=\pi }\pm 50`$ cm-1 we obtain $`J/J_{}1.251.35`$.

Up to now we have derived a set of exchange couplings which consistently describes the spin gap and the $`S`$=0 two-triplet bound state. Now we compare the spectral density calculated for this parameter set with the line shape of the entire spectrum observed in $`\sigma (\omega )`$. This provides a thorough test whether we indeed have found the minimal model which captures all relevant properties. We consider an isolated $`\mathrm{Cu}_2\mathrm{O}_3`$-ladder, for which the major contribution to $`\sigma (\omega )`$ results from the simultaneous excitation of two neighboring spin-flips plus a Cu-O bond-stretching phonon:

$$\sigma (\omega )\omega \underset{p_x}{}\underset{p_y=0,\pi }{}f_p\mathrm{Im}\delta B_p;\delta B_p_{(\omega \omega _{ph})}$$
(13)
where $`\delta B_p^{\mathrm{leg}}=\frac{1}{N}_i_{k=l,r}e^{ipr_{i,k}}\left(𝐒_{i,k}𝐒_{i+1,k}𝐒_{i,k}𝐒_{i+1,k}\right)`$ and $`\delta B_p^{\mathrm{rung}}=\frac{1}{N}_ie^{ipr_i}\left(𝐒_{i,l}𝐒_{i,r}𝐒_{i,l}𝐒_{i,r}\right)`$ are the spin-flip operators for polarization of the electrical field along the legs and the rungs, respectively. To lowest order (4th order in the Cu-O hopping $`t_{pd}`$), the dominant contribution to the phonon form factor $`f_p^{\mathrm{leg}}`$ comes from the in-phase and the out-of-phase stretching modes of the O-ions on the legs, whereas for $`f_p^{\mathrm{rung}}`$ the out-of-phase stretching mode and the vibration of the O-ion on the rung are taken into account:

$$f_p^{\mathrm{leg}}=8\mathrm{sin}^4(\frac{p_x}{2}),f_p^{\mathrm{rung}}=8\mathrm{sin}^2(\frac{p_x}{2})+4$$
(14)
So far we have used the Lanczos-vector method (LVM)  to calculate the dispersion of the bound state. This was justified by the fact that the LVM uses the first Lanczos vectors as target states besides the ground state and consequently represents the bound state quite accurately. However, in order to calculate the high-frequency range of $`\sigma (\omega )`$ the LVM is not reliable. We use the correction-vector method , because here the density matrix is optimized to represent the correction vector for each frequency point separately.

In Fig. 3 we compare $`\sigma (\omega )`$ for polarization along the legs and along the rungs, obtained by the correction vector method for a ladder with N=80-sites, with the experimental result of $`\mathrm{La}_{5.2}\mathrm{Ca}_{8.8}\mathrm{Cu}_{24}\mathrm{O}_{41}`$ . Here we have used $`J/J_{}`$=1.3 and $`J_{cyc}/J_{}`$=0.2 with $`J_{}=1000\mathrm{cm}^1`$ . For polarization along the legs, $`\sigma (\omega )`$ contains two contributions, where the two legs are excited in-phase ($`p_y=0`$) or out-of-phase ($`p_y=\pi `$). The in-phase mode contains two- and four-triplet excitations, as can be inferred from the rung-triplet representation of the spin-flip operator $`𝐒_{i,l}𝐒_{i+1,l}+𝐒_{i,r}𝐒_{i+1,r}`$ in Eq. (12). It includes the $`S`$=0 two-triplet bound state as discussed above. The out-of-phase mode reflects excitations of three triplets due to $`𝐒_{i,l}𝐒_{i+1,l}𝐒_{i,r}𝐒_{i+1,r}=\frac{i}{2}_{\alpha ,\beta ,\gamma }ϵ_{\alpha \beta \gamma }\{t_{i,\alpha }^{}t_{i+1,\beta }^{}t_{i+1,\gamma }s_i+t_{i,\beta }^{}t_{i+1,\alpha }^{}s_{i+1}t_{i,\gamma }\mathrm{H}.\mathrm{c}.\}`$ and hence contributes only to the high-energy continuum excitations . For polarization along the rungs, $`\sigma (\omega )`$ contains only two-triplet excitations (see Eq. (12) for the rung-triplet representation of $`𝐒_{i,l}𝐒_{i,r}`$) and the lower bound state is suppressed due to a selection rule .

Comparing the DMRG results for polarization along the legs with the experimental spectra we find excellent agreement within the range of the bound state (which has been used in order to fix the coupling constants), and also with respect to the overall line shape of the high-energy continuum excitations. For this polarization the experimental determination of the precise spectral weight of the continuum is difficult due to an electronic background . The small differences between experiment and theory in the continuum range are certainly within the experimental error bars. For polarization along the rungs, where experimentally the magnetic contribution to $`\sigma (\omega )`$ can be identified unambiguously and with high precision , the continuum is reproduced nearly perfectly. This provides strong evidence that a nearest-neighbor Heisenberg Hamiltonian with $`J/J_{}`$=1.3 and an additional cyclic spin exchange of $`J_{cyc}/J_{}`$=0.2 represents the magnetic excitations of the spin-ladder compound (La,Ca)14Cu24O41.

In spite of the overall impressive consistency of experiment and theory, some discrepancies are visible: (i) In $`\sigma _{\mathrm{leg}}(\omega )`$ the upper bound state (for $`p_x\pi /2`$) is suppressed much stronger in the DMRG calculation than in the experimental spectrum. (ii) In $`\sigma _{\mathrm{rung}}(\omega )`$ the experimental data displays a pronounced shoulder at about 2200 cm-1, i.e., slightly above the suppressed lower bound state, which is not present in our DMRG calculation. We propose that these discrepancies are related to different phonon form factors which arise due to the actual arrangement of the ladders within the trellis lattice, and in higher order (with respect to Cu-O hopping) in the external photon-electron-phonon vertex. Considering, e.g., an additional term $`\mathrm{sin}^2(p_x/2)`$ in $`f_p^{\mathrm{leg}}`$ will enhance the weight of the upper bound state in $`\sigma _{\mathrm{leg}}(\omega )`$ with respect to the lower one. Such considerations, however, depend on details of the non-local excitation process in these ladder systems and are therefore not of fundamental relevance.

In summary, we have demonstrated that the inclusion of a cyclic spin exchange $`J_{cyc}`$ is necessary to obtain a consistent interpretation of optical experiments and the magnitude of the spin gap in (La,Ca)14Cu24O41. The dispersion of the $`S`$=0 bound state is very sensitive to the magnitude of $`J_{cyc}`$. Quantitative investigation yields $`J_{cyc}/J_{}0.20`$-0.27 and $`J/J_{}1.25`$-1.35 with $`J_{}950`$-1100 cm-1. With these parameters we obtain excellent agreement with the experimental spectra of $`\mathrm{La}_{5.2}\mathrm{Ca}_{8.8}\mathrm{Cu}_{24}\mathrm{O}_{41}`$, both for the bound state and the continuum. This gives strong evidence that we indeed have identified the minimal model which captures all relevant magnetic properties.

We would like to thank M. Greiter, A. P. Kampf and G. S. Uhrig for stimulating discussions. This project was supported by the DFG (SFB 484 and SFB 608) and by the BMBF (13N6918/1).


format: text

text:
# Constraining the variation of the coupling constants with big bang nucleosynthesis

## I Introduction

Big bang nucleosynthesis (BBN) is one of the most important tools for probing the early universe. The standard big bang nucleosynthesis (SBBN) predicts the primordial abundances of light elements (D, 3He, 4He and 7Li) with the only free parameter, the baryon-to-photon-ratio $`\eta `$. For $`\eta O(10^{10})`$, they are remarkably consistent with the observations of the abundance which extend over 9 digits. Requiring not to vitiate this consistency, the stringent constraints on various theories which affect cosmology have been obtained.

Among such theories constrained by BBN, there are some in which values of coupling constants may vary. The SBBN prediction assumes, along with three light neutrino species and no lepton asymmetry, that physical parameters involved in the calculation (the fine structure constant $`\alpha _{EM}`$, the Fermi constant $`G_F`$, the electron mass $`m_e`$, etc.) are the same for now and the BBN time. However, this is not an obvious choice because there are numbers of ways in which they could have varied. In addition to the theoretical possibility, the recent analysis of quasar absorption lines found possible evidence for variation in $`\alpha _{EM}`$ .

In this paper, we consider a model taken from string theory where these coupling constants are all related to the expectation values of dilaton field $`\mathrm{\Phi }`$ and could in principle vary with time. For example, very general arguments indicate that the dilaton cannot be stabilized at a value we would characterise as corresponding to weak coupling . Also, in certain popular models for stabilising the dilaton using gaugino condensates, the cosmological evolution would almost inevitably tend to overshoot the desired minimum of the dilaton potential and run off to an anti-de Sitter vacuum .

The notion that the coupling constants are determined by a single scalar field not only motivates the investigation of their time variation as a probe for physics beyond the Standard Model, but also makes the analysis of simultaneous changes in the couplings simple and concrete. In addition to the free parameter of the SBBN $`\eta `$, we have only another parameter $`{\displaystyle \frac{\mathrm{\Delta }\mathrm{\Phi }}{\mathrm{\Phi }}}{\displaystyle \frac{\mathrm{\Phi }_{BBN}\mathrm{\Phi }_{now}}{\mathrm{\Phi }_{now}}}`$, the fractional variation of the dilaton field. After we calculate the primordial abundances with various values of these two quantities, we search a parameter region not excluded by the observation. Thus obtained limit on $`\mathrm{\Delta }\mathrm{\Phi }/\mathrm{\Phi }`$ is readily translated in the limits on the coupling constants.

Our analysis is positioned as the extension to Ref. . There, the dilaton dependence of the coupling constants is determined by the action of the heterotic string in the Einstein frame. The constraint is obtained from 4He whose abundance and coupling dependence can be estimated without recourse to a numerical calculation. Their method is easy to calculate and appropriate for an order estimation. But we can extract more information by performing the numerical calculation, whose advantages are: 1. the abundances of the light elements other than 4He are calculated so we can use more observational data, especially D’s, to put a constraint, 2. we can take into account of $`\eta `$ dependence of the abundances, 3. we use realistic values for the reaction rates so we avoid the rough estimation of the weak reaction rate (both absolute value and temperature dependence) based on dimensional analysis, 4. by Monte-Carlo simulation, we can estimate theoretical uncertainty, 5. by calculating a statistical measure (such as $`\chi ^2`$), we can objectively quantify the constraint. Of course, as the SBBN becomes more and more precise owing to the progress in the studies of nuclear reaction rates and of primordial abundances, numerical computation is necessary in general to reflect the recent development and to produce a result close to the truth.

In the next section, we review briefly the SBBN and see what kind of physical quantities are needed to make the prediction. In Sec. III, we estimate their coupling dependences. In Sec. IV, we introduce a model containing the dilaton and we investigate the dilaton dependence of the coupling constants. In Sec. V, we put a bound on $`\mathrm{\Delta }\mathrm{\Phi }/\mathrm{\Phi }`$ by calculating $`\chi ^2`$ with observational and theoretical uncertainty from the recent data (the former very much dominates). In terms of the limit on $`\mathrm{\Delta }\alpha _{EM}/\alpha _{EM}`$, we find

$$6.0\times 10^4<\frac{\mathrm{\Delta }\alpha _{EM}}{\alpha _{EM}}<1.5\times 10^4,$$
(1)
two orders more restrictive than the one found in Ref.  where only the variation of $`\alpha _{EM}`$ is taken into account. We discuss this originates from the unified gauge couplings, which relates electromagnetic sector to strong sector. At the energy scale of BBN, manifest quantities are $`\alpha _{EM}`$ and $`\mathrm{\Lambda }_{QCD}`$ which appears instead of $`\alpha _{strong}`$ through dimensional transmutation. The former linearly depend on the dilaton but the latter does exponentially. 4He abundance is determined by the magnitude of neutron-proton mass difference $`\mathrm{\Delta }m`$ which is in turn determined by $`\mathrm{\Lambda }_{QCD}`$ and this dominates the constraints by BBN.

## II Physical quantities in the SBBN

In order to see how the effects of changing coupling constants arise, we summarise the main points of the SBBN calculation. Following Ref. , we divide it into three stages. The actual calculation is performed by computer which runs the code to solve the set of ordinary differential equations so, of course, does not distinguish such stages but we see the important quantities which determine the primordial abundances of the light elements appear in this brief review.

1st Stage (statistical equilibrium)($`T1`$ MeV; $`t1`$ sec)

The energy density and the number density are dominated by the relativistic electron $`e^{}`$, positron $`e^+`$, neutrino $`\nu `$, anti-neutrino $`\overline{\nu }`$ and photon $`\gamma `$. There are three types of neutrino and anti-neutrino. There is only tiny fraction of protons and neutrons, just $`10^9`$ number fraction. In this period, all these particles make elastic scattering so frequently that they are in thermal equilibrium and has equal temperatures. Especially, weak interaction cross sections are large and neutrinos are in thermal equilibrium with the other particles. In addition, there are weak interaction processes interchanging protons and neutrons,

$`n+\nu _ep+e^{}`$ (2)
$`n+e^+p+\overline{\nu }_e`$ (3)
$`np+e^{}+\overline{\nu }_e,`$ (4)
so they are also in the chemical equilibrium. Then, the number ratio of proton to neutron at temperature $`T`$ is $`\frac{n_n}{n_p}=e^{(m_nm_p)/T}e^{(\mu _n\mu _p)/T}`$ where $`m_n(m_p)`$ is neutron (proton) mass and $`\mu _n(\mu _p)`$ is its chemical potential, but the last factor $`e^{(\mu _n\mu _p)/T}`$ is almost unity if the universe has no lepton asymmetry, which is the assumption of SBBN. So as long as equilibrium holds with temperature $`T`$, the neutron-to-proton ratio is determined by the neutron-proton mass difference $`\mathrm{\Delta }m=m_nm_p`$,

$$\frac{n_n}{n_p}=e^{\mathrm{\Delta }m/T}.$$
(5)
At this stage, the nuclear reaction rates are also fast enough for light nuclei to be in chemical equilibrium. Then the abundances of the nuclear species with mass number $`A`$ and atomic number $`Z`$ is

$`Y_A`$ $`=`$ $`g_A[\zeta (3)^{A1}\pi ^{(1A)/2}2^{(3A5)/2}]A^{5/2}\left({\displaystyle \frac{T}{m_N}}\right)^{3(A1)/2}`$ (7)
$`\times \eta ^{A1}Y_p^ZY_n^{AZ}\mathrm{exp}(B_A/T),`$
where $`\zeta (3)1.202`$, $`m_N`$ is the nucleon mass, $`B_A`$ is the binding energy and $`g_A`$ is the number of degrees of freedom. Since $`\eta `$ is small ($`10^{10}`$) and $`B_AT`$, nuclei abundance is negligible at this temperature.

2nd Stage (neutron-proton freeze-out) ($`T0.7`$ MeV; $`t2`$ sec)

As the universe expands and cools, the rates of interactions involving neutrinos decreases. Especially, there are some point when the reactions (2)-(4) practically stop and the chemical equilibrium break down. After this point, the numbers of neutrons and protons do not change (actually, small amount of neutrons turns into protons by the beta decay) and their ratio is fixed (froze out) at the value determined from Eq.(5) with the temperature of that epoch. We call this ”freeze-out temperature” and denote $`T_f`$.

The freeze-out occurs when the rate of the reaction $`\mathrm{\Gamma }(T)`$ and the expansion rate of the universe $`H(T)`$ becomes equal:

$$H(T_f)\mathrm{\Gamma }_{weak}(T_f).$$
(8)
The expansion rate is determined by the energy density $`\rho `$ through the Friedman equation,

$$H^2=\frac{8\pi G\rho }{3}.$$
(9)
At this stage, $`\rho `$ is dominated by photons, neutrinos and relativistic electrons ($`T>m_e=0.511`$ MeV) and hence the cosmic density is given by $`\rho =10.75(\pi ^2/30)T^4`$. The rate of weak interaction is

$`\mathrm{\Gamma }_{weak}`$ $`=`$ $`\mathrm{\Gamma }(np)`$
$`=`$ $`\mathrm{\Gamma }(n+\nu _ep+e^{})`$
$`+\mathrm{\Gamma }(n+e^+p+\overline{\nu }_e)+\mathrm{\Gamma }(np+e^{}+\overline{\nu }_e).`$
Using the Fermi theory of the weak interaction, this can be written as

$`\mathrm{\Gamma }(np)`$
$`=`$ $`A\left[{\displaystyle _0^{\mathrm{}}}+{\displaystyle _{\mathrm{}}^{\mathrm{\Delta }mm_e}}+{\displaystyle _{\mathrm{\Delta }m+m_e}^0}\right]dx`$
$`\times \left(1{\displaystyle \frac{m_e^2}{(x+\mathrm{\Delta }m)^2}}\right)^{\frac{1}{2}}(x+\mathrm{\Delta }m)^2x^2`$
$`\times {\displaystyle \frac{1}{1+e^{x/T_\nu }}}{\displaystyle \frac{1}{1+e^{(x+\mathrm{\Delta }m)/T_e}}}.`$
We calculate the normalization factor $`A`$ by using $`\mathrm{\Gamma }(np+e^{}+\overline{\nu }_e)|_{T=0}=\tau _n^1`$ where $`\tau _n`$ is the neutron lifetime. Integration gives $`A=0.05606\tau _n^1`$ MeV-4. From these expression and Eq.(8), we find $`T_f0.7`$ MeV and, from Eq. (5), freeze-out ratio $`\left(n_n/n_p\right)_f=e^{\mathrm{\Delta }m/T_f}0.158`$.

For the nuclear reaction, the situation did not change from the previous stage. Their abundances are still very small.

3rd Stage (light-element synthesis) (0.7 MeV $`>T>`$ 0.05 MeV; 3 sec $`<t<`$ 6 min)

After the second stage, the electron-positron pair annihilation completed so there remains the photons and the neutrinos as relativistic particles.

Since the early universe has lower density than inside of stars, there practically occurs only two body reactions. Thus unless the deuterons are synthesized by the reaction $`p+n\mathrm{D}+\gamma `$, larger nuclei are not synthesized. This reaction does not proceed until the photon density is low enough not to photodissociate D or, in other words, the factor $`\eta \mathrm{exp}(B_\mathrm{D}/T)`$ in Eq.(7) becomes $`O(1)`$. It occurs at about $`T=0.06`$ MeV. After that, mainly charged-particle reactions such as $`\mathrm{D}+\mathrm{D}t+p`$ and $`t+\mathrm{D}{}_{}{}^{4}\mathrm{He}+n`$ proceed synthesizing almost all of the neutrons into 4He because it has the largest binding energy among the light elements. Therefore, 4He abundance (conventionally expressed by mass ratio) can be estimated by the frozen n-p ratio found at the 2nd stage:

$$Y_{{}_{}{}^{4}\mathrm{He}}\frac{2}{1+(n_p/n_n)_f}=\frac{2}{1+e^{\mathrm{\Delta }m/T_f}}$$
(10)
Besides 4He, small amount of D and 3He and very small amount of 7Li are synthesized. Their abundances are determined by $`\eta `$ and the reaction rates.

In summary, to perform BBN calculation, we need to input following values at BBN time:

* neutron-proton mass difference $`\mathrm{\Delta }m`$
* neutron lifetime $`\tau _n`$
* nuclear reaction rates

Mainly, the abundance of 4He are affected by the first two as seen from Eqs. (8) and (10), and D, 3He and 7Li by the last one.

## III The coupling constants dependence of the BBN inputs

We know the present value of the BBN input quantities pointed out in the previous section with some experimental uncertainty. To estimate their values during the BBN, we try to express these quantities in terms of the coupling constants.

### A neutron-proton mass difference $`\mathrm{\Delta }m`$

The origin of the neutron-proton mass difference is traced to the electromagnetic self-energy difference and the d-u quark mass difference (QCD chiral symmetry breaking by the mass terms). The former makes proton heavier than neutron but the latter does the inverse and the total is measured accurately to be $`\mathrm{\Delta }mm_nm_p=1.2933318\pm 0.0000005`$ MeV .

As argued in Ref.  the largest contribution to the electromagnetic part comes from the Born term of the Cottingham formula  and hence it can be calculated with relatively less uncertainty. This formula expresses the nucleon self-energy in the first order of $`\alpha _{EM}`$ by its electric and magnetic form factor. Using the formula and experimental data for the form factors, the difference in the neutron and the proton is $`0.76`$MeV. This value is proportional to (of course) $`\alpha _{EM}`$ and to $`\mathrm{\Lambda }_{QCD}`$ to have a proper dimension.

On the other hand, the absolute values of u and d quark masses are not well known, $`m_u=15`$ MeV and $`m_d=39`$ MeV . But knowing the electromagnetic contribution, the quark contribution should be 2.05 MeV. This is proportional to the Yukawa couplings $`y_u,y_d`$ and the Higgs expectation value $`H`$.

Now we can write

$`\mathrm{\Delta }m`$ $`=`$ $`a\alpha _{EM}\mathrm{\Lambda }_{QCD}+b(y_dy_u)H,`$ (11)
where $`a`$ and $`b`$ are constants we assume not to depend on any coupling constants. Thus, we estimate the mass difference at BBN epoch $`\mathrm{\Delta }m_{bbn}`$ as

$`\mathrm{\Delta }m_{bbn}`$ (12)
$`=`$ $`a\alpha _{EMbbn}\mathrm{\Lambda }_{QCDbbn}+by_{bbn}H_{bbn}`$ (13)
$`=`$ $`0.76{\displaystyle \frac{\alpha _{EMbbn}}{\alpha _{EMnow}}}{\displaystyle \frac{\mathrm{\Lambda }_{QCDbbn}}{\mathrm{\Lambda }_{QCDnow}}}+2.05{\displaystyle \frac{y_{bbn}}{y_{now}}}{\displaystyle \frac{H_{bbn}}{H_{now}}}[\mathrm{MeV}].`$ (14)
### B neutron lifetime

Neutron $`\beta `$-decay $`np+e+\overline{\nu }`$ is very well approximated by the one-point interaction of four particles: neutron, proton, electron and neutrino. Its coupling constant is denoted $`G_F`$ called the Fermi coupling constant. Using this theory,

$`\tau _n^1`$ $``$ $`G_F^2{\displaystyle d^3p_ed^3p_\nu \delta (E_e+E_\nu +m_pm_n)}`$ (16)
$`=`$ $`G_F^216\pi ^2{\displaystyle _{m_e}^{\mathrm{\Delta }m}}𝑑E_eE_e\sqrt{E_e^2m_e^2}(\mathrm{\Delta }mE_e)^2`$ (17)
$`=`$ $`G_F^216\pi ^2m_e^5{\displaystyle \frac{1}{60}}\{\sqrt{q^21}(2q^49q^28)`$ (19)
$`+15q\mathrm{log}(q+\sqrt{q^21})\}`$
where we defined $`q\frac{\mathrm{\Delta }m}{m_e}`$. Denoting the factor $`\frac{1}{60}\{\mathrm{}\}f(q)`$,

$`\tau _n`$ $``$ $`G_F^2m_e^5f(q)^1`$ (20)
$`=`$ $`H^4(y_eH)^5f(q)^1=H^1y_e^5f(q)^1.`$ (21)
The first equality follows from $`G_F=\frac{g_2^2}{M_W^2}=\frac{g_2^2}{(g_2H)^2}=\frac{1}{H^2}`$ where $`g_2`$ is the SU(2) coupling constant and $`M_W`$ is the weak boson mass. Therefore, we obtain

$`\tau _{n,bbn}=\left[{\displaystyle \frac{y_{e,bbn}}{y_{e,now}}}\right]^5\left[{\displaystyle \frac{H_{bbn}}{H_{now}}}\right]^1\left[{\displaystyle \frac{f(q_{bbn})}{f(q_{now})}}\right]^1,`$ (22)
where $`y_e`$ is the electron Yukawa coupling and $`f(q_{now})=1.63615`$.

### C charged-particle induced reaction rates

Most of the reaction rates involved in the BBN calculation are charged-particle induced reaction rates. Their $`\alpha _{EM}`$ dependence is considered in Ref.  for the rates discussed in Ref. . We implemented the $`\alpha _{EM}`$ dependence in the same manner updating the reaction rates recently compiled by Angulo et al. . Since only a few reaction rates have resonance terms which is considered to be dependent of the strong coupling constant, or in other words since the charged-particle reaction rates influential for BBN are practically determined by Coulomb barrier penetrability, we expect that neglecting their strong coupling dependence will not affect the results a lot.

### D neutron-induced reaction rates

The cross section for $`n+p\mathrm{D}+\gamma `$ is calculated at energies relevant to BBN using the effective field theory that describes the two-nucleon sector . The rate obtained by thermal averaging this theoretical cross section reproduces the abundance calculation using the rate shown in Ref. . We exploit this theoretical formula to estimate the coupling dependence by assuming the parameters which have the dimension \[length\]n in the formula proportional to $`m_\pi ^n`$. The coupling dependence of pion mass $`m_\pi `$ is known from Gell-Mann-Oakes-Renner relation to be $`m_\pi ^2m_{quark}\mathrm{\Lambda }_{QCD}`$. The other parameters in the formula are the nucleon mass $`m_N`$ and the deuteron binding energy $`B_D`$. Their dependence is estimated by $`m_N\mathrm{\Lambda }_{QCD}`$ and $`B_D\frac{m_\pi ^2}{\mathrm{\Lambda }_{QCD}}m_{quark}`$. This $`B_D`$’s dependence is estimated considering a simple square well potential for deuteron. It is also estimated by uncertainty principle: $`p\frac{1}{\mathrm{\Delta }x}m_\pi `$, and $`B_D+\frac{p^2}{2m_N}0`$.

For the other neutron-induced reactions, 3He$`(n,p)t`$ and 7Be$`(n,p)^7`$Li, we use the rate fitted from the experimental data as found in Ref.  for the former and Ref.  for the latter, and they are assumed to be independent of the couplings because there are no theoretical derivation for these. Because 3He$`(n,p)t`$ affects 3He abundance which is not used for our analysis, our result does not change by neglecting its coupling dependence. On the other hand, 7Be$`(n,p)^7`$Li affects 7Li abundance a lot when $`\eta `$ is large so its effect may be large. However, since there is large observational uncertainty in 7Li, it is expected that our result does not change but we have to know its coupling dependence to extract a reliable constraint from 7Li data when the uncertainty of 7Li abundances decreases in the future.

## IV The dilaton dependence of the coupling constants

Next, we introduce the action which governs the variation of the coupling constants as is used in Ref. . This is the tree level low energy action of the heterotic string in the Einstein frame.

$`S={\displaystyle }d^4x\sqrt{g}({\displaystyle \frac{1}{2\kappa ^2}}R{\displaystyle \frac{1}{2}}_\mu \mathrm{\Phi }_\mu \mathrm{\Phi }{\displaystyle \frac{1}{2}}_\mu \varphi _\mu \varphi `$ (23)
$`\mathrm{\Omega }^2V(\varphi )\overline{\psi }\gamma ^\mu D_\mu \psi \mathrm{\Omega }^1m_\psi \overline{\psi }\psi `$ (24)
$`{\displaystyle \frac{\alpha ^{}}{16\kappa ^2}}\mathrm{\Omega }^2F_{\mu \nu }F^{\mu \nu }).`$ (25)
where $`\mathrm{\Phi }`$ is the dilaton field, $`\varphi `$ is an arbitrary scalar field and $`\psi `$ is an arbitrary fermion. $`D_\mu `$ is the gauge covariant derivative corresponding to gauge fields with field strength $`F_{\mu \nu }`$. $`\kappa ^2=8\pi G`$ and $`\mathrm{\Omega }=e^{\kappa \mathrm{\Phi }/\sqrt{2}}`$ which is the conformal factor used to move from the string frame. Powers of $`\mathrm{\Omega }(\mathrm{\Phi })`$ multiplying terms in the action indicate the dilaton dependence of the coupling constants and masses.

More concretely, $`\varphi `$ is Higgs field and $`V(\varphi )`$ is its potential which we assume to be given by hand (as is done in the Standard Model). The all over $`\mathrm{\Omega }`$ factor before the scalar potential means that the Higgs vacuum expectation value $`H`$ is independent of the dilaton so has same value for now and BBN time. $`H`$ is taken to be constant in our calculation.

$`F_{\mu \nu }`$ is the gauge field with gauge group including $`SU(3)\times SU(2)\times U(1)`$. We define its Lagrangian density as $`\frac{1}{4g^2}F_{\mu \nu }F^{\mu \nu }`$ where $`g`$ is the unified coupling constant. Comparing with Eq. (25),

$$\frac{1}{g^2(M_p^2)}=\frac{\alpha ^{}e^{\sqrt{2}\kappa \mathrm{\Phi }}}{4\kappa ^2}\frac{\alpha ^{}S}{4\kappa },$$
(26)
where we define

$$S\frac{e^{\sqrt{2}\kappa \mathrm{\Phi }}}{\kappa }.$$
(27)
For each value of $`S`$, we can calculate the gauge coupling constants at the low energy using renormalization group equations. $`\alpha _{EM}`$ almost does not run so,

$$\alpha _{EM}(M_{BBN})\alpha _{EM}(M_p)=\frac{g(M_p)^2}{4\pi }=\frac{\kappa }{\pi \alpha ^{}}S^1.$$
(28)
Therefore,

$$\frac{\alpha _{EM,bbn}}{\alpha _{EM,now}}=\left[\frac{S_{bbn}}{S_{now}}\right]^1=\frac{1}{1+v_S}.$$
(29)
where we define fractional $`S`$ variation,

$$v_S=\frac{S_{bbn}S_{now}}{S_{now}}.$$
(30)
From the solution of 1-loop RGE for the SU(3) coupling constant $`g_3`$, of which integration constant is determined by $`g_3(\mathrm{\Lambda }_{QCD})=\mathrm{}`$,

$$\frac{\kappa }{\pi \alpha ^{}}S^1=\frac{g_3(M_p)^2}{4\pi }\frac{12\pi }{27\mathrm{log}(M_p^2/\mathrm{\Lambda }_{QCD}^2)}$$
(31)
or

$$\mathrm{\Lambda }_{QCD}=M_p\mathrm{exp}\left(\frac{2\pi ^2\alpha ^{}S}{9\kappa }\right).$$
(32)
Therefore, we obtain

$`{\displaystyle \frac{\mathrm{\Lambda }_{QCD,bbn}}{\mathrm{\Lambda }_{QCD,now}}}`$ $`=`$ $`\mathrm{exp}\left({\displaystyle \frac{2\pi ^2\alpha ^{}[S_{bbn}S_{now}]}{9\kappa }}\right)`$ (33)
$`=`$ $`\mathrm{exp}\left({\displaystyle \frac{8\pi ^2}{9g(M_p)^2}}v_S\right),`$ (34)
where we use $`g(M_p)^2=0.1`$.

Finally, $`\psi `$’s are the ordinary SM leptons and quarks. As we take $`H=const.`$, Yukawa couplings $`y`$ depend on dilaton as $`e^{\kappa \mathrm{\Phi }/\sqrt{2}}`$. In terms of $`S`$,

$$\frac{y_{bbn}}{y_{now}}=\frac{1}{\sqrt{1+v_S}}.$$
(35)
## V Constraints on the variation of the coupling constants

Using the model described in the previous section, we can express the coupling dependence of the BBN input parameters considered in Sec. III with the fractional variation of the dilaton $`v_S`$ defined in Eqs. (27) and (30). In order to quantify how much variation is consistent with the observation, we calculate the abundances for different values of $`v_S`$ in addition to $`\eta `$ with the standard BBN code . Then we calculate $`\chi ^2(\eta ,v_S)`$ as,

$$\chi ^2=\underset{i}{}\frac{(a_i^{th}a_i^{obs})^2}{(\sigma _i^{th})^2+(\sigma _i^{obs})^2}$$
(36)
where $`i`$ is the type of the element with which we try to put a constraint. To estimate theoretical errors, we have performed 1000 Monte-Carlo simulations using the values of Ref.  for the nuclear reaction rates uncertainty and Particle Data Group  for the neutron lifetime $`885.7\pm 0.8`$ s. For the observational errors for D, 4He and 7Li, we adopt:

$`(\mathrm{D}/\mathrm{H})^{obs}`$ $`=`$ $`(3.0\pm 0.4)\times 10^5,`$ (37)
$`Y^{obs}`$ $`=`$ $`0.238\pm (0.002)_{stat}\pm (0.005)_{syst},`$ (38)
$`\mathrm{log}_{10}[(^7\mathrm{Li}/\mathrm{H})^{obs}]`$ $`=`$ $`9.76\pm (0.012)_{stat}\pm (0.05)_{syst}`$ (40)
$`\pm (0.3)_{add}.`$
Eq. (37) is taken from Ref.  and Eq. (38) from Ref.  where the first error is the statistical uncertainty and the second error is the systematic one. Eq. (40) is from Ref.  with the error we have added for the uncertainty in chemical evolution .

The results are shown in Fig. 1. The shape of the contour lines are easily understood. The contours drawn from the $`\chi ^2`$ of three elements (Fig. 1 (a)) are just the product set of D (Fig. 1 (b)) and 4He (Fig. 1 (c)) because 7Li (Fig. 1 (d)) does not give much constraint.

D and 7Li do not vertically constrain much and 4He alone constrain $`v_S`$. The reason is D and 7Li abundances are determined mainly by charged-particle reactions while 4He by freeze-out ratio. The former is determined by Coulomb barrier penetrability or electro-magnetic coupling $`\alpha _{EM}`$ and the latter by neutron-proton mass difference $`\mathrm{\Delta }m`$. $`\mathrm{\Delta }m`$ is depend linearly on $`\mathrm{\Lambda }_{QCD}`$ as described in Sec. III A. $`\alpha _{EM}`$ and $`\mathrm{\Lambda }_{QCD}`$ are related to each other through Eqs. (29) and (31) which came from general notion of the gauge coupling unification. This relation tells that $`{\displaystyle \frac{\mathrm{\Delta }\mathrm{\Lambda }_{QCD}}{\mathrm{\Lambda }_{QCD,now}}}{\displaystyle \frac{80\pi ^2}{9}}{\displaystyle \frac{\mathrm{\Delta }\alpha _{EM}}{\alpha _{EM,now}}}`$ where $`\mathrm{\Delta }\alpha _{EM}=\alpha _{EM,bbn}\alpha _{EM,now}`$ etc. so small increase in $`\alpha _{EM}`$ is accompanied by great increase in $`\mathrm{\Lambda }_{QCD}`$ and $`\mathrm{\Delta }m`$. This is the reason why 4He constrain $`v_S`$ much more than D and 7Li do in our model.

To make sure, we explain the trend of the 4He contour. As $`v_S`$ increases, $`\mathrm{\Lambda }_{QCD}`$ decreases and hence $`\mathrm{\Delta }m`$ increases. This makes $`\frac{n_n}{n_p}|_{freezeout}`$ decrease, so 4He abundance decreases. Since the 4He abundance is an increase function of $`\eta `$, the increase in $`v_S`$ relaxes the constraint for higher $`\eta `$. This shows up in the trend that the contour goes up in the direction of increasing $`\eta `$.

The constraint on the dilaton field variation is obtained from Fig. 1 (a). For 95% confidence level,

$$1.5\times 10^4<v_S<6.0\times 10^4,$$
(41)
or using $`S=e^{\sqrt{2}\kappa \mathrm{\Phi }}/\kappa `$ with $`\kappa =(8\pi G)^{1/2}=(2.43\times 10^{18}`$GeV$`)^1`$, we obtain

$$1.0\times 10^{14}\mathrm{GeV}<\mathrm{\Delta }\mathrm{\Phi }<2.6\times 10^{14}\mathrm{GeV}.$$
(42)
## VI Conclusion

In summary, we have considered the BBN calculation with the coupling constants taking different values from the present ones and whose variation is determined by the dilaton field $`\mathrm{\Phi }`$ as in the action (25). The set up is the same as Ref.  but we have performed full nuclear reaction network calculation with the BBN code so the abundances of D and 7Li in addition to 4He are available to compare with the observation and make a constraint. This is done by calculating $`\chi ^2`$ as a function of the baryon-to-photon ratio $`\eta `$ and the fractional variation of the dilaton $`v_S`$ defined in Eqs. (27) and (30), and we have obtained the constraint (41), which is the same order as the value found in Ref. . It might be interesting that the standard BBN ($`v_S`$=0) is outside the allowed region with 70% confidence level as seen in Fig. 1 (a) although it is not statistically significant.

The constraint on $`v_S`$ is readily converted to that on $`\alpha _{EM}`$ using Eq.(29), which becomes $`\mathrm{\Delta }\alpha _{EM}/\alpha _{EM}=v_S`$ when $`v_S`$ is small, and gives

$$6.0\times 10^4<\frac{\mathrm{\Delta }\alpha _{EM}}{\alpha _{EM}}<1.5\times 10^4,$$
(43)
which is two orders stricter than what found in Ref.  where only $`\alpha _{EM}`$ is varied. This feature should be common in the models in which variation of the gauge couplings are related to each other through gauge unification so the variation in $`\mathrm{\Lambda }_{QCD}`$ is much larger than that of $`\alpha _{EM}`$.

For comparison, we draw a contour diagram when only $`\alpha _{EM}`$ is varied. The result is Fig. 2. The shape is similar but noticing the vertical scale, it is two orders larger than the Fig. 1 (a). The constraint becomes $`5.0\times 10^2<\mathrm{\Delta }\alpha _{EM,only}/\alpha _{EM,only}<1.0\times 10^2`$ and this is similar to the result obtained in Ref.  (factor difference is attributed to the difference in the adopted observational data). The contours in the Fig. 2 completely covers the contours in the Fig. 1 (a). This is another demonstration that the constraint is not obtained by the variation of $`\alpha _{EM}`$ directly but by $`\mathrm{\Lambda }_{QCD}`$.

Finally, we notice that the limit (43) is consistent and has favourable trend with the result found in Ref.  from the quasar absorption systems, that is, $`\mathrm{\Delta }\alpha _{EM,only}/\alpha _{EM,only}=(0.72\pm 0.18)\times 10^5`$. This shows $`\alpha _{EM}`$ was smaller at $`z=O(1)`$. So if the variation is monotone, it should be further smaller at the BBN time. This is just the tendency of Eq. (43) (the allowed region is wider for the negative side). More quantitative comparison needs a model for time evolution of the coupling constants (for the model considered here, it is specifying dilaton potential) and the analysis of the quasar data under the situation where the couplings other than $`\alpha _{EM}`$ vary.

We have not reached a stage where such a constraint as (42) can be used to say something about particle physics. In this paper, we bear string theory in mind and treat the time variation of the coupling constants as caused by the dilaton, scalar field found in the theory. However, in cosmology, there flourishes scalar fields to attack important problems (e.g. inflation, quintessence) so the similar analysis (with different actions) can be used to constrain such scenarios. In order to find useful application, the analysis itself has to be sharpened because there still remains unsatisfactory estimates concerning the nuclear force. This sector will be improved with the field theoretic (QCD based) understanding of the nuclear force.


text:
# Statistical Geometry of Packing Defects of Lattice Chain Polymer from Enumeration and Sequential Monte Carlo Method

## I Introduction

Soluble proteins are well-packed, their packing densities may be as high as that of crystalline solids . Yet there are numerous packing defects or voids in protein structures, whose size distributions are broad . The volume ($`v`$) and area ($`a`$) of protein does not scale as $`va^{3/2}`$, which would be expected for models of tight packing. Rather, $`v`$ and $`a`$ scale linearly with each other . In addition, the scaling of protein volume and cluster-radius  is characteristic of random sphere packing. Such scaling behavior indicates that the interior of proteins is more like Swiss cheese with many holes than tightly packed jigsaw puzzles .

What effects do void have? Proteins are often very tolerant to mutations , which may suggest potentially stabilizing roles of voids in proteins. Voids in proteins are also often associated with protein function. The binding sites of proteins for substrate catalysis and ligand interactions are frequently prominent voids and pockets on protein structures . However, the energetic and kinetic effects of maintaining specific voids in proteins are not well-understood, and the shape space of voids of folded and unfolded proteins are largely unknown.

In this paper, we examine the details of the statistical nature of voids in simple lattice polymers. Lattice models have been widely used for studying protein folding, where the conformational space of simplified polymers can be examined in detail . Despite its simplistic nature, lattice model has provided important insights about proteins, including collapse and folding transitions , influence of packing on secondary structure formation , and designability of lattice structures . However, one drawback is that lattice model is not well-suited for studying void-related structural features, such as protein functional sites, since it is not easy to model the geometry of voids.

In this article, we first define voids as topological defects and describe a simple algorithm for void detection in two-dimensional lattice. We then enumerate exhaustively the conformations for all $`n`$-polymers up to $`n=22`$, and analyze the relationship of probability of void formation, expected packing density and compactness, as well as expected wall interval of void with chain length. To study statistical geometry of long chain polymers, we describe a Monte Carlo sampling strategy under the framework of Sequential Importance Sampling, and introduce the technique of resampling. The results of simulation of long chain polymers up to $`N=200`$ for several geometric parameters are then presented. We further explore the conformation reduction factor $`R`$ of void formation, and describe the significant end-effect of void formation, as well as the scaling law of $`R`$ and wall interval of voids. In the final section, we summarize our results and discuss effective sampling strategy for studying the conformational space of voids.

## II Lattice Model and Voids

Lattice polymers are self-avoiding walks (SAWs), which can be obtained from a chain-growth model . Specifically, an $`n`$-polymer $`P`$ on a two-dimensional square lattice $`Z^2`$ is formed by monomers $`n_i,i\{1,..,N\}`$. The location $`x_i`$ of a monomer $`n_i`$ is defined by its coordinates $`x_i=(a_i,b_i)`$, where $`a_i`$ and $`b_i`$ are integers. The monomers are connected as a chain, and the distance between bonded monomers $`x_i`$ and $`x_{i+1}`$ is $`1`$. The chain is self-avoiding: $`x_ix_j`$ for all $`ij`$. We consider the beginning and the end of a polymer to be distinct. Only conformations that are not related by translation, rotation, and reflection are considered to be distinct. This is achieved by following the rule that a chain is always grown from the origin, the first step is always to the right, and the chain always goes up at the first time it deviates from the $`x`$-axis. For a chain polymer, two non-bonded monomers $`n_i`$ and $`n_j`$ are in topological contact if they intersect at an edge that they share. If two monomers share a vertex of a square but not an edge, these two monomers are defined as not in contact.

(Figure 1 here.)

When the number of monomer $`n`$ is 8 or more, a polymer may contain one or more void (Figure 1a). We define voids as topological features of the polymer. The complement space $`Z^2P`$ that is not occupied by the polymer $`P`$ can be partitioned into disjoint components:

$`Z^2P=V_0\dot{}V_1\mathrm{}\dot{}V_k.`$Here $`V_0`$ is the unique component of the complement space that extends to infinity. We call this the outside. The rest of the components that are disjoint or disconnected to each other are voids of the polymer. Because non-bonded monomers intersecting at a vertex are defined as not in contact, they do not break up the complement space. As an example, the unfilled space contained within the polymer in Figure 1b is regarded as one connected void of size $`4`$ rather than two disjoint voids of size $`2`$. A simple algorithm for void detection can be found in Appendix. Figure 2 shows the only six conformations among all 301,100,754 conformations of 22-mer found to have 4 voids.

(Figure 2 here.)

## III Voids Distribution by Exact Enumeration

Probability of Forming Voids and Expected Number of Voids.  The number of conformations $`\omega (n)`$ for $`n`$-polymer up to $`n=25`$ by exhaustive enumeration is shown in Table I. The numbers of conformations for polymers up to $`n=15`$ are in exact agreement with those reported in Chan and Dill . Table I also lists the number of conformations $`\omega _k(n)`$ containing $`k=1,2,3`$, or $`4`$ voids.

The probability for a polymer to form one or more voids $`\pi _v`$ is calculated as:

$`\pi _v={\displaystyle \frac{_{i=1}^k\omega _k(n)}{\omega (n)}}.`$The expected number of voids $`\overline{n}_v`$ for a polymer is:

$`\overline{n}_v={\displaystyle \frac{_{i=1}^k\omega _k(n)k}{\omega (n)}}.`$As the chain length grows, it is clear that both $`\pi _v`$ and $`\overline{n}_v`$ increases (Figure 3a and Figure 3b).

(Figure 3 here.)

Void Size.  The total size $`v`$ of voids in a polymer is the sum of the sizes of all voids, namely, the total number of all unoccupied squares that are fully contained within the polymer. Let $`\omega _v(n)`$ be the number of conformations of $`n`$-polymer with total void size $`v`$. The expected total void size $`\overline{v}`$ for $`n`$-polymer is:

$`\overline{v}={\displaystyle \frac{_v\omega _v(n)v}{\omega (n)}}.`$Figure 3c shows that the expected void size $`\overline{v}`$ increases with chain length $`n`$.

Wall Size of Void.  For a void $`V`$ of size $`v`$, what is the required minimum length $`l(v)`$ for a polymer that can form such a void? Equivalently, what is the size of the wall of the polymer containing void $`V`$? Here we first restrict our discussion to voids formed only by strongly connected unoccupied sites, namely, any neighboring two sites of a void must be sharing at least one edge of the squares. We exclude voids containing weakly connected sites, where two neighboring sites are connected by only one shared vertex (Figure 1b). For $`v=1,2`$ and $`3`$, it is easy from the geometry of the voids to see that $`l(v)=8,10`$ and $`12`$, respectively. However, in general $`l(v)`$ also depends on the shape of the void. A void of size $`4`$ can have five different shapes. If the void is of the shape of a $`2\times 2`$ square, $`l(4)=12`$. For the other four shapes, $`l(4)=14`$.

For any strongly connected void, we find that the following general recurrence relationship for $`l(v)`$ holds:

$`l(v)=l(v1)+\{\begin{array}{cc}2,& \text{if }\mathrm{\Delta }V=3\\ 0,& \text{if }\mathrm{\Delta }V=2\\ 1,& \text{if }\mathrm{\Delta }V=1,\end{array}`$where $`V`$ represents the boundary edges of void $`V`$, and $`\mathrm{\Delta }V`$ represents the net gain in the number of boundary edges introduced by the newly added unoccupied site. Although the number and explicit shapes of strongly connected voids of size up to 5 can be found in , there is no general analytical formula known for the number of shapes of a void of size $`v`$. This is related to the problem of determining the number of polyominos or animals (as in percolation theory) of a given size.

When weakly connected voids are also considered, there are more possible wall sizes for void. For $`22`$-mer, the number of different wall sizes observed for a void, strongly or weakly connected, at various size are shown in Figure 4a. Voids of size 5 has the largest diversity in wall size. This is of course due to the fixed chain length. A short chain such as the 22mer has only a small number of ways for form large voids. Figure 4b shows the average wall size for various void size in 22-mer. The expected or average wall size $`\overline{w}(n)`$ for a void in a $`n`$-polymer can be calculated as:

$`\overline{w}(n)={\displaystyle \underset{v}{}}w{\displaystyle \frac{\omega _{v,w}(n)}{\omega _v(n)}},`$where $`v`$ is the void size, $`w`$ the wall size of the void, $`\omega _{v,w}(n)`$ is the number of $`n`$-polymers containing a void of size $`v`$ with wall size $`w`$, and $`\omega _v(n)`$ is the total number of $`n`$-polymers with a void of size $`v`$. Figure 3d shows that $`\overline{w}(n)`$ increases with chain length. Wall size and void size are analogous to the area and volume of voids in three dimensional space.

Packing Density.  An important parameter that describes how effectively atoms fill space is the packing density $`p`$. In proteins, it is defined by Richards and colleagues as the amount of the space that is occupied within the van der Waals envelope of the molecule, divided by the total volume of space that contains the molecule . It has been widely used by protein chemists as a parameter for characterizing protein folding . Following this original definition, the packing density $`p`$ for lattice polymer is:

$`p=n/(n+v),`$when a $`n`$-polymer has a total void size of $`v`$.

The expected packing density $`\overline{p}(n)`$ for a $`n`$-polymer can be calculated as:

$`\overline{p}(n)={\displaystyle \underset{p}{}}p{\displaystyle \frac{\omega _p(n)}{\omega (n)}},`$where $`\omega (n)`$ is the number of all conformations of $`n`$-mer, $`\omega _p(n)`$ the number of $`n`$-mers with packing density of $`p`$. The scaling of $`\overline{p}(n)`$ with the chain length $`n`$ decreases roughly linearly between $`n=7`$ and $`n=22`$ (Figure 5a). Because it takes at least two additional monomer to increase the size of a void by one, $`\overline{p}(n)`$ decreases only when $`n`$ is an odd number for short chains.

Although voids are packing defects, most conformations with voids have high packing density, namely, the total size of voids are small. Among all conformations of 22-mer containing one void, the number of conformation increases monotonically with packing density. The lowest packing density $`0.52`$ has only $`11`$ conformations, whereas the highest packing density $`0.92`$ has the largest number ($`6,756,751`$) of conformations (Figure 5c). Similar relationship is found among conformations with 2 and 3 voids (Figure5c).

Compactness.  Another important parameter that measures packing of lattice polymer is the number of nonbonded contacts $`t`$. It is related to the compactness parameter $`\rho `$, defined in  as $`\rho =t/t_{\mathrm{max}}`$, where $`t_{\mathrm{max}}`$ is the maximum number of nonbonded contact possible for a $`n`$-polymer. Compactness $`\rho `$ has been studied extensively in seminal works by Chan and Dill . Although $`p`$ is sometimes correlated with the compactness $`\rho `$, these two parameters are distinct. The relationship between compactness and expected packing density for chain polymer of length $`1422`$ is shown in Figure 5. For all chain lengths, both maximally compact polymer ($`\rho =1`$) and extended polymer ($`\rho =0`$) have maximal packing density ($`p=1`$). Polymers with $`\rho `$ between 0.4 and 0.6 have lowest packing density and therefore tend to have larger void size. The explanation is simple. An extended lattice chain polymer has no voids, it therefore achieves maximal packing density of $`p=1`$, but its compactness $`\rho `$ is 0. A maximally compact polymer with $`\rho =1`$ also contains no voids, its $`p`$ is $`1`$. On the other hand, non-maximally compact polymers can have a range of packing densities.

(Figure 5 here.)

## IV Obtaining Void Statistics for Long Chains via Importance Sampling

Sequential Importance Sampling. Geometrically complex and interesting features emerge only in polymers of sufficient length, which are not accessible for analysis by exhaustive enumeration, due to the fact that the number of possible SAWs increases exponentially with the chain length. Monte Carlo methods are often used to generate samples from all possible conformations and obtain estimates of feature statistics using those samples. However, when chain length becomes large, the direct generation of SAWs using rejection method (i.e., generate random walks on the lattice and only accept those that are self-avoiding) from the uniform distribution of all possible SAW’s becomes difficult. The success rate $`s_N`$ of generating SAWs decreases exponentially: $`s_NZ_N/(4\times 3^{N1})`$. For $`N=48`$, $`s_N`$ is only $`0.79\%`$ . To overcome this attrition problem, a widely used approach is the Rosenbluth Monte Carlo method of biased sampling . The task is to grow one more monomer for a $`t`$-polymer chain that has been successfully grown from 1 monomer after $`t1`$ successive steps without self-crossing, until $`t=n`$, the targeted chain length. In this method, the placement of the $`(t+1)`$-th monomer is determined by the current conformation of the polymer. If there are $`n_t`$ unoccupied neighbors for the $`t`$-th monomer, we then randomly (with equal probability) set the $`(t+1)`$-th monomer to any one of the $`n_t`$ sites. However, the resulting sample is biased toward more compact conformations and does not follow the uniform distribution. Hence each sample is assigned a “weight” to adjust for the bias. Any statistics can then be obtained from weighted average of the samples. In the case of Rosenbluth chain growth method, the weight is computed recursively as $`w_t=n_tw_{t1}`$.

Liu and Chen  provided a general framework of Sequential Monte Carlo (SMC) methods which extend the Rosenbluth method to more general setting. Sophisticated but more flexible and effective algorithms can be developed under this framework. In the context of growing polymer, SMC can be formulated as follows. Let $`(x_1,\mathrm{},x_t)`$ be the position of the $`t`$ monomers in a chain of length $`t`$. Let $`\pi _1(x_1),\pi _2(x_1,x_2),\mathrm{},\pi _t(x_1,\mathrm{},x_t)`$ be a sequence of target distributions, with $`\pi (x_1,\mathrm{},x_n)=\pi _n(x_1,\mathrm{},x_n)`$ being the final objective distribution from which we wish to draw inference from. Let $`g_{t+1}(x_{t+1}x_1,\mathrm{},x_t)`$ be a sequence of trial distributions which dictates the growing of the polymer. Then we have:

| Procedure SMC ($`n`$) |
| --- |
| Draw $`x_1^{(j)}`$, $`j=1,\mathrm{},m`$ from $`g_1(x_1)`$ |
| Set the incremental weight $`w_1^{(j)}=\pi _1(x_1^{(j)})/g_1(x_1^{(j)})`$ |
| for  $`t=1`$ to $`n1`$ |
|  | for  $`j=1`$ to $`m`$ |
|  |  | // Sampling for the $`(t+1)`$-th monomer for the $`j`$-th sample |
|  |  | Draw position $`x_{t+1}^{(j)}`$ from |
|  |  |  | $`g_{t+1}(x_{t+1}|x_1^{(j)}\mathrm{}x_t^{(j)})`$ |
|  |  | // Compute the incremental weight. |
|  |  | $`u_{t+1}^{(j)}\frac{\pi _{t+1}\left(x_1^{\left(j\right)}\mathrm{}x_{t+1}^{\left(j\right)}\right)}{\pi _t\left(x_1^{\left(j\right)}\mathrm{}x_t^{\left(j\right)}\right)g_{t+1}\left(x_{t+1}^{\left(j\right)}|x_1^{\left(j\right)}\mathrm{}x_t^{\left(j\right)}\right)}`$ |
|  |  | $`w_{t+1}^{(j)}u_{t+1}^{(j)}w_t^{(j)}`$ |
|  | endfor |
|  | Resampling |
| endfor |

At the end, the configurations of successfully generated polymers $`\{(x_1^{(j)},\mathrm{},x_n^{(j)})\}_{j=1}^m`$ and their associated weights $`\{w_n^{(j)}\}_{j=1}^m`$ can be used to estimate any properties of the polymers, such as expected void size, compactness, and packing density. That is, the objective inference $`\mu _h=E_\pi [h(x_1,\mathrm{},x_n)]`$ is estimated with

$$\widehat{\mu }_h=\frac{_{j=1}^mh(x_1^{(j)},\mathrm{},x_n^{(j)})w_n^{(j)}}{_{j=1}^mw_n^{(j)}},$$
(1)
for any integrable function $`h`$ of interests.

The critical choices that affect the effectiveness of the SMC method are: (1) the approximating target distribution $`\pi _t(x_1\mathrm{}x_t)`$, (2) the sampling distribution $`g_{t+1}(x_{t+1}|x_1\mathrm{}x_t)`$, and (3) the resampling scheme. In this study, we are interested in sampling from the uniform distribution $`\pi _n(x_1\mathrm{}x_n)`$ of all geometrically feasible conformations of length $`n`$, which we call the final objective distribution. It can also be chosen to be the Boltzmann distribution when energy function such as the HP model  is introduced.

The Rosenbluth method  is a special case of SMC. Its target distributions $`\pi _t(x_1\mathrm{}x_t)`$ is the uniform distribution of all SAWs of length $`t`$. Its sampling distribution $`g_{t+1}(x_{t+1}|x_1\mathrm{}x_t)`$ is the uniform distribution among all $`n_1(x_1,\mathrm{},x_t)`$ unoccupied neighboring sites of the last monomer $`x_t`$, and the weight function is

$`w(x_1,\mathrm{},x_t,x_{t+1})=w(x_1,\mathrm{},x_t)n_1(x_1,\mathrm{},x_t).`$When there is no unoccupied neighboring sites ($`n_1(x_1,\mathrm{},x_t)=0`$), there is no place to place the $`(t+1)`$-th monomer. In this case, the chain runs into a dead end and we declare the conformation dead, with weight assigned to be $`0`$. In the case of Rosenbluth method, no resampling is used.

Similarly, the $`k`$-step look ahead algorithm  chooses $`\pi _{t+1}(x_1,\mathrm{},x_{t+1})`$ being the marginal distribution of $`\pi _{t+k}^{}(x_1,\mathrm{},x_{t+k})`$, the uniform distribution of all SAW’s of length $`t+k`$. Hence $`\pi _{t+1}`$ is closer to the final objective distribution – the uniform distribution of all SAW’s of length $`n`$. Specifically,

$`\pi _{t+1}(x_1,\mathrm{},x_{t+1})`$ $`=`$ $`{\displaystyle \underset{x_{t+2},\mathrm{},x_{t+k}}{}}\pi _{t+k}^{}(x_1,\mathrm{},x_{t+1},x_{t+2},\mathrm{},x_{t+k})`$
$``$ $`n_k(x_1,\mathrm{},x_{t+1})`$
where $`n_k(x_1,\mathrm{},x_{t+1})`$ is the total number of SAWs of length $`t+k`$ “grown” from $`(x_1,\mathrm{},x_{t+1})`$ \[i.e. with the first $`(t+1)`$ positions at $`(x_1,\mathrm{},x_{t+1})`$.\] In the $`k`$-step look-ahead algorithm, the sampling distribution is

$`g_{t+1}(x_{t+1}=xx_1,\mathrm{},x_t)={\displaystyle \frac{n_k(x_1,\mathrm{},x_t,x)}{n_{k+1}(x_1,\mathrm{},x_t)}}.`$It chooses the next position according to what will happen $`k`$ steps later. Namely, the probability of placing the $`t+1`$-th monomer at $`x`$ is determined by the ratio of the total number of SAWs of length $`t+k`$ grown from $`(x_1,\mathrm{},x_t,x)`$ and the total number of SAWs of the same length $`t+k`$ grown from one step earlier $`(x_1,\mathrm{},x_t)`$. The corresponding weight function is

$`w(x_1,\mathrm{},x_t,x_{t+1})`$ $`=`$ $`{\displaystyle \frac{n_k(x_1,\mathrm{},x_{t+1})}{n_k(x_1,\mathrm{},x_t)\frac{n_k(x_1,\mathrm{},x_{t+1})}{n_{k+1}(x_1,\mathrm{},x_t)}}}`$
$`=`$ $`{\displaystyle \frac{n_{k+1}(x_1,\mathrm{},x_t)}{n_k(x_1,\mathrm{},x_t)}}.`$
Although it has higher computational cost, it usually produces better inference on the final objective distribution, with less “dead” conformations. The standard Rosenbluth algorithm is $`1`$ step look ahead algorithm.

To compare geometric properties estimated from sequential Monte Carlo method and those obtained by exhaust enumeration, we examine the expected number of voids and expected void size for polymer from chain length $`14`$ to $`22`$. Figure 6 shows that sequential Monte Carlo can provide very accurate estimation of these geometric properties of voids. Here 2-step look ahead is used, with Monte Carlo sample size of 100,000 and no resampling.

(Figure 6 here.)

The resampling step is one of the key ingredient of the SMC . There are many cases where resampling is beneficial. First, note that it is unavoidable to have some dead conformations during the growth. These chains need to be replaced to maintain sufficient Monte Carlo sample size. Second, the weight of some chains may become relatively so small that their contribution in the weighted average (1) is negligible. When the variance of the weights is large, the effective Monte Carlo sample size becomes small . Third, for a specific function $`h`$, its value may become too small (even zero) for some sampled conformations. In all these cases, efficiency can be gained by replacing those conformations with “better” ones. This procedure is called “resampling”. There are many different ways to do resampling. One approach is rejection control , which regenerates the replacement conformations from scratch. An easier approach is to duplicate the existing and good conformations . Specifically,

| Procedure Resampling |
| --- |
| //$`m`$: number of original samples. |
| //$`\{(x_1^{(j)},\mathrm{},x_t^{(j)}),w^{(j)}\}_{j=1}^m`$: original properly weighted samples |
| for $`j=1`$ to $`m`$ |
|  | Set resampling probability of $`j`$th conformation $`\alpha ^{(j)}`$ |
| endfor |
| for $`j=1`$ to $`m`$ |
|  | Draw $`j`$th sample from original samples $`\{(x_1^{(j)},\mathrm{},x_t^{(j))}\}_{j=1}^m`$ |
|  |  | with probabilities $`\{\alpha ^{(j)}\}_{j=1}^m`$ |
|  | //Each sample in the newly formed sample is assigned a new weight. |
|  | //$`j`$-th chain in new sample is a copy of $`k`$-th chain in original sample. |
|  | $`w^{(j)}w^{(k)}/\alpha ^{(k)}`$ |
| endfor |

In the resampling step, the $`m`$ new samples $`\{(x_1^{(j)},\mathrm{},x_t^{(j)}\}_{j=1}^m`$ can be obtained either by residual sampling or by simple random sampling. In residual sampling, we first obtain the normalized probability $`\stackrel{~}{\alpha }^{(j)}=\alpha ^{(j)}/\alpha ^{(j)}`$. Then $`[m\stackrel{~}{\alpha }^{(j)}]`$ copies of $`j`$-th sample are made deterministically for $`j=1,\mathrm{},m`$. For the remaining $`m[m\stackrel{~}{\alpha }^{(j)}]`$ samples to be made, we randomly sample from the original set with probability proportional to $`m\stackrel{~}{\alpha }^{(j)}[m\stackrel{~}{\alpha }^{(j)}]`$.

The choice of resampling probability proportional to $`\alpha ^{(j)}`$ is problem specific. For general function $`h`$, such as the end-to-end extension $`x_nx_1`$, it is common to use $`\alpha ^{(j)}=w_t^{(j)}`$. In this case, all the samples in the new set have equal weight. When the function is irregular, a carefully chosen set of $`\alpha ^{(j)}`$ will increase the efficiency significantly.

The method of pruning and enriching of Grassberger  is a special case of the residual sampling, with $`\alpha ^{(j)}=0`$ for the $`k`$ chains with zero weight (dead conformations), $`\alpha ^{(j)}=2`$ for the top $`k`$ chains with largest weights, and $`\alpha ^{(j)}=1`$ for the rest of the chains. Residual sampling on this set of $`\alpha `$ is completely deterministic. The resulting sample consists of two copies of the top $`k`$ conformations (each of them having half of their original weight) and one copy of the middle $`n2k`$ chains with their original weight. The $`k`$ dead conformations are removed.

In our study of the relationship between compactness and packing density, we use a more flexible resampling method. Our focus is on the packing density among all conformations with certain range of compactness. In this case, our object target distribution is the uniform distribution among all possible SAW’s with compactness measure falling within a certain interval, i.e., a truncated distribution. Although compactness changes slowly as the chain grows, to grow into a long chain it is possible that the compactness of a chain evolve and cover a wide range during growth. Hence we choose the uniform distribution of all possible SAWs’ of length $`t`$ as our target distribution at $`t`$, and only select those with the desire compactness at the end for our estimation of the packing density. In order to have higher number of usable samples (i.e.,, to achieve better acceptance rate) at the end, we encourage growth of chains with desirable compactness through resampling. Specifically,

| Procedure Resampling ($`m,d,c_t`$) |
| --- |
| // $`m`$: Monte Carlo sample size, $`d`$: steps of looking-back. |
| // $`c_t`$: targeting compactness. |
| $`k`$ number of dead conformations. |
| Divide $`mk`$ samples randomly into $`k`$ groups. |
| for  group $`i=1`$ to $`k`$ |
|  | Find conformations not picked in previous $`d`$ steps. |
|  |  | //Pick the best conformation $`P_j`$, for example |
|  |  | $`P_j`$ polymer with $`\mathrm{min}|cc_t|`$ |
|  | Replace one of $`k`$ dead conformations with $`P_j`$ |
|  | Assign both copies of $`P_j`$ half its original weight. |
| endfor |

Here $`d`$ is used to maintain higher diversity for resampled conformations.

Most polymers sampled by sequential Monte Carlo without resampling are well-extended with few voids, as shown in Figure 7(a) and (b). They have small compactness (less than 0.5) and large packing density. As a result, a small number of samples are accepted at the end whose compactness falls within the desired interval of higher than 0.5. By using the resampling step described above, we were able to generate more samples near the desired compactness value of 0.6 (Figure 7c). Figure 7c is a pure histogram of compactness in the observed samples, without regarding the weight of the samples. Figure 7d shows that the resampling technique is also very effective in shifting the samples to small packign density values, hence improve the inferences.

(Figure 7 here.)

## V Voids Distribution of Long Chains

We apply the techniques of sequential Monte Carlo with resampling to study the statistical geometry of voids in long chain polymers. Figure 8a shows that the probability of void formation increases with the chain length. At chain length 105–110, about half of the conformations contain voids. The expected number of voids (Figure 8b) increases linearly with chain length. Similar linear scaling behavior is also observed in proteins . The expected wall size of void and void size also increase with chain length (Figure 8c and Figure 8d).

The expected packing density is found to decrease with chain length, which is consistent with the scaling relationship of void size and chain length shown in Figure 8c. The compactness $`\rho `$ of chain polymer has been the subject of several studies . The asymptotic value of $`\rho `$ we found is 0.18, slightly different from that reported in  ($`\rho =0.16`$), and is within the range of 0.16 – 0.24 reported in .

(Figure 8 here.)

To explore the relationship of packing density $`p`$ and compactness $`\rho `$, we use sequential Monte Carlo with 2-step look-ahead to sample 200,000 conformations, each with appropriate weight assigned. This is repeated 20 times, and the weighted average values of packing density at various compactness for chains with 60–100 monomers are plotted (Figure 9). The compactness value corresponding to the minimum packing density seems to have shifted from $`0.462`$ for 22-mer by enumeration to above 0.5 for 100-mer by sampling. However, the overall pattern of $`p`$ and $`\rho `$ found by Monte Carlo is very similar to the pattern found by enumeration for polymers with $`N<22`$.

The accuracy of geometric properties of long chain polymers estimated by Monte Carlo can be assessed by the variance obtained from multiple Monte Carlo runs.

(Figure 9 here.)

## VI End Effects of Void Formation

What is the effect of void formation on the size of conformational space? We consider the conformational reduction factor of voids. Following , we define the conformational reduction factor due to the constraint of a void as:

$`R(n;i,j)={\displaystyle \frac{\omega (n;i,j)}{\omega (n)}},`$where $`\omega (n;i,j)`$ is the number of conformations that contains a void beginning at monomer $`(i)`$ and ending at monomer $`(j)`$, and $`\omega (n)`$ is the total number of conformations of $`n`$-polymers. $`R(n;i,j)`$ reflects the restriction of conformational space due to the formation of a void with wall interval of $`k=|ij|`$. Figure 10a shows a 24-mer with one void that starts at $`i=4`$ and $`k=19`$. Unlike self-contacts or self-loops, which was subject of detailed studies by Chan and Dill , all conformations analyzed here must contain a void. The polymer shown in Figure 10 with a large loop has no void, and such polymers do not contribute to the numerator of $`R`$.

(Figure 10 here.)

Figure 11a shows the reduction factor $`R`$ calculated by enumeration for voids at different starting positions with wall intervals $`k=7,9`$ and $`11`$. There are clearly strong end-effects: The reduction factor of voids of the same wall interval depends on where the void is located. $`R`$ decreases rapidly as the void moves from the end of chain towards the middle. Void formation is much more preferred at the end of chain. Similar end effects of void formation are also observed for 55-mer sampled by sequential Monte Carlo (Figure 11b).

(Figure 11 here.)

The end-effect of voids has the same origin as the end-effect of self-contact, which has been extensively studied by Chan and Dill . Because of the effect of excluded volume, sterically it is less hindering to form a void at the end of a polymer. When a void is formed, the conformational space of the $`k+1`$ monomers between monomer $`i`$ and $`j`$, as well as the two tails become restricted. When void is formed at chain end, only one tail is subject to conformational restriction.

Void formation is different from self-contact. When monomer $`i`$ and $`j`$ form self-contact, it may involve the formation of a void, but it is also possible that there will be no unfilled space between $`i`$ and $`j`$. When a void is formed beginning at monomer $`i`$ and ending at monomer $`j`$, some monomers between $`i`$ and $`j`$ will have unsatisfied contact interactions. Compare to non-bonded self-contact, the effect of conformation reduction is more pronounced for void formation. For two-dimensional lattice, the ratios between reduction factors of self-contact at chain end and mid chain of a sufficiently long polymer are $`1.3,1.4,1.5`$ and $`1.6`$ for $`k=3,5,7`$ and $`9`$, respectively , whereas the ratios for voids at chain end and the symmetric midpoint of $`N=22`$ polymer are $`3.4,4.0`$, and $`4.4`$. for $`k=7,9`$ and $`11`$. The conformational reduction factor $`R(i,j)`$ for voids at various beginning positions $`i`$ and various ending position $`j`$ can be summarize in a two-dimensional contour plot as shown in Figure 11b.

We now consider the power-law dependence of $`R(N;i,j)`$ on the wall interval $`k=|ij|`$. In the studies of self-contacting loops by Chan and Dill , the scaling exponent $`\nu `$ of the reduction factor $`R`$ and loop length $`k=|ij|`$ for $`R(N;i,j)k^\nu `$ is found to be dependent both on $`k`$ and the location of the cycle in the chain. The values of $`\nu `$ for self-contact range from 1.6 when $`k=N`$ to 2.4 when the loop is in the middle of a long chain with two long tails. Because void formation involves at least 8 monomers, its scaling behavior is less amenable to exhaust enumeration, and application of Monte Carlo sampling is essential. Based on estimations from Monte Carlo simulation of void formation in 50-mer, $`\nu `$ ranges from $`1.4\pm 0.2`$ for $`l_0=1`$ to $`3.0\pm 0.2`$ for void initiation position $`l_0=8`$ (Figure 11c). Our results show that the scaling exponent of $`R`$ with $`k=|ij|`$ for void formation is similar to that of self-contacting loop. This scaling exponent also depends on the location of the void.

## VII Conclusion

In this work, we have studied the statistical geometry of voids as topological features in two-dimensional lattice chain polymers. We define voids as unfilled space fully contained within the polymer, and have developed a simple algorithm for its detection. We have explored the relationship of various statistical geometric properties with the chain length of the polymer, including the probability of void formation $`\pi _v`$, the expected number of voids $`\overline{n}_v`$, the expected void size $`\overline{v}`$, the expected wall size of voids $`\overline{w}`$, packing density $`p`$, and the expected compactness $`\rho `$. Our results show that for chains of $`>`$105-110 monomers, at least half of the conformations contain a void. At about 150 monomers, there will be at least one void expected in a polymer. The expected wall size scale linearly with the chain length, and about 10% of the monomers participate in the formation of voids. We formalize the concept of packing density for lattice polymers. We found that both the packing density and compactness decrease with chain length. The asymptotic value of compactness $`\rho `$ is estimated to be 0.18.

We have also characterized the relationship of packing density and compactness, two parameters that have been used frequently for studying protein packing. Our results indicate that packing density reaches minimum values between compactness 0.4 – 0.6. The entropic effects of voids are studied by analyzing the conformational reduction factor $`R`$ of void formation. We found that there is significant end-effect for void formation: the ratio of $`R`$ at chain end and at mid chain may be twice as large as that of the $`R`$ factor for contact loops, where the formation of voids is not required.

In this study, we have applied sequential Monte Carlo sampling and resampling techniques to study the statistical geometry of voids. Sequential Monte Carlo sampling and resampling is essential for exploring the geometry of long chain polymers. This is a very general approach that allows the generation of increased number of conformations with interesting characteristics. For example, we can replace dead conformations with existing conformations of highest weight, or conformations with highest compactness, or with smallest radius of gyration. Figure 12a shows the histograms of conformation of 100-mer at different packing density generated without resampling. Figure 12c shows the histograms of conformations when resampling by weight and resampling by compactness $`\rho `$ are used. Other resampling schemes are possible, e.g. resampling by radius-of-gyration, by packing density. During resampling, the number $`k`$ of dead conformations at each step of growth is identified and these are replaced with conformations of interest from $`k`$ randomly divided groups. These conformations must have not been resampled in previous 4 steps of the growth process to maintain sample diversity. Both histograms where resampling is used deviate from that of Figure 12a. Resampling by weight shifts the peak of the conformations to below 0.2, and resampling by compactness turns the histogram into bi-modal. The latter produces a lot more conformations with compactness $`\rho >0.4`$.

SMC sampling and resampling use biased samples since conformations are generated with probability different from that of the target distribution. The bias is ictated by different method of resampling and different choices of the number of steps of look-ahead in sequential Monte Carlo. An essential component of a successful biased Monte Carlo sampling is the appropriate weight assignment to each sample conformation. This is necessary because we need to estimate the expected values of parameter such as packing density and void size under the target distribution of all geometrically feasible conformations. In Figure 12a where each of the 200,000 starting conformations is generated by two-step look-ahead without resampling, not every conformation is generated with the same probability and therefore is assigned different weight accordingly. Figure 12b shows the weight-adjusted histogram, which is indicative of the probability density function at different compactness for the population of all geometrically feasible 100-mers.

Figure 12d shows that when weights are incorporated and the area of the histogram normalized to the final number of surviving conformations, the weighted distributions of conformations using different resampling techniques have excellent agreement with the weighted distribution when no resampling is used (Figure 12b). This example shows that by incorporating weights, the target distributions can be faithfully recovered even when the sampling is very biased.

(Figure 12 here.)

Although sequential Monte Carlo sampling is very effective, the estimation of parameters associated with rare events remain difficult. In Figure 11 where conformational reduction factor $`R`$ is plotted at various void initiation position and wall interval length, voids starting at position of 1 but with odd wall intervals ($`k\{11,13,\mathrm{},25\}`$) are much rarer, and it is unlikely sequential Monte Carlo sampling with limited sample size can provide large enough effective sample size for the accurate estimation of scaling parameters $`\nu `$, where $`R(N;i,j)k^\nu `$.

In this study, we are interested in the statistics of void geometry, and our target distribution is the uniform distribution of all conformations of length $`n`$. With the introduction of appropriate potential function and alphabet of monomers such as the HP model , we can study the thermodynamics, kinetics, and sequence degeneracy of chain polymers when voids are formed in polymers. In these cases, our target distributions will be chain polymers under the Boltzmann distribution derived from the corresponding potential functions.

## VIII Acknowledgments

This work is supported by funding from National Science Foundation DMS 9982846, CMS 9980599, DMS 0073601, DBI0078270, and MCB998008, and American Chemical Society/Petroleum Research Fund.

## IX Appendix

To detect voids in a polymer, we use a simple search method. For an $`l\times l`$ lattice, we start from the lower-left corner. Once we found an unoccupied site $`u`$, we use the breadth-first-search (BFS) method to identify all other unoccupied sites that are connected to site $`u`$. These sites are grouped together and marked as “visited”. Collectively they represent one void in the lattice. We continue this process until all unoccupied sites are marked as visited:

| Algorithm VoidDetection (lattice, $`l`$) |
| --- |
| $`v=0`$ // Number of voids |
| for  $`i=1`$ to $`l`$ |
|  | for  $`j=1`$ to $`l`$ |
|  |  | if site($`i,j`$) is unoccupied and not visited |
|  |  |  | $`vv+1`$ |
|  |  |  | Mark ($`i,j`$) as visited. |
|  |  |  | BreadthFirstSearch(lattice, ($`i,j`$)) |
|  |  |  | Update the size of $`void(i,j)`$ |
|  |  | endif |
|  | endfor |
| endfor |

Details of BFS can be found in algorithm textbooks such as .