text:
# Optical measurements of the superconducting gap in MgB2

11institutetext: 1. Physikalisches Institut, Universität Stuttgart, Pfaffenwaldring 57, D-70550 Stuttgart, Germany 22institutetext: Laboratory of Solid State Physics, Materials Science Centre, Nijenborgh 4, 9747 AG Groningen, The Netherlands 33institutetext: Department of Physics, Aoyama-Gakuin University, 6-16-1 Chitsedai, Setagaya-ku, Tokyo 157, Japan (Received: date 7 March 2001 / Revised version: )

## Abstract

Far-infrared reflectivity studies on the polycrystalline intermetallic compound MgB2 with a superconducting transition temperature $`T_c=39`$ K were performed at temperatures 20 K to 300 K. We observe a significant raise of the superconducting-to-normal state reflectivity ratio below $`70\mathrm{cm}^1`$, with a maximum at about $`2530\mathrm{cm}^1`$, which gives a lower estimate of the superconducting gap of $`2\mathrm{\Delta }(0)`$$`34`$ meV.

Recently superconductivity in the binary intermetallic compound MgB2 with a $`T_c`$ close to 40 K was reported by Akimitsu et al. Akimitsu01 . Most studies performed on this compound yet indicate that MgB2 consistently behaves as a phonon mediated superconductor within the framework of the BCS theory, probably in a strong coupling limit Budko01 ; Finnemore01 ; Jung01 ; Kotegawa01 . Up to now no optical investigations of this material were reported, and this is in part due to the poor quality of the presently available samples, which are sintered polycrystals. It is well known, however, that optical spectroscopy is an extremely powerful method for studying the materials in the superconducting (SC) state and is able of providing information on such important parameters as the SC energy gap, penetration depth, coherence effects, scattering mechanism, etc. Tinkham96 ; Timusk89 .

In this short note we report on far-infrared reflectivity measurements performed on a sintered pellet of MgB2. We are fully aware that these experiments bear inherent problems, mainly connected with the surface roughness, which we cannot overcome at this point. Nevertheless, we have attempted to extract important information on the intrinsic properties of MgB2 from our reflectivity data.

MgB2 has a hexagonal structure with $`P6/mmm`$ symmetry. It crystallizes in the so-called AlB2 structure where the boron atoms are located at a primitive honeycomb lattice, consisting of graphite-type sheets. The structure of MgB2 is shown in Fig. 1; the dimensions of the unit cell are $`a`$=3.086 Å and $`c`$=3.524 Å, according to our X-ray analysis. The borons span hexagonal prisms; the large, almost spherical pores are filled by Mg which acts as a spacer. Similar to graphite, the distance between the boron planes is larger by a factor of two compared to the intraplanar B-B bonds, and hence the B-B bonding is strongly anisotropic (two-dimensional).

Powder of high purity MgB2 was pressed in a pellet and treated in Ar atmosphere at $``$900C; details on the sample preparation are reported in Ref. Akimitsu01 . The quality of the sintered samples was checked by X-ray analysis, resistivity, and susceptibility measurements. In Fig. 1 the temperature dependence of the dc resistivity and of the magnetic susceptibility are plotted. The SC phase transition is characterized by a width of less than 1 K from resistivity and of around 5 K from the susceptibility measurements. At $`T`$=300 K we find $`\rho `$=$`115\mu \mathrm{\Omega }`$cm, and just above the SC transition temperature $`T_c`$=39 K the resistivity has decreased by a factor of 2. The resistivity in the normal state can best be fitted by a power-law temperature dependence $`T^\alpha `$ with $`\alpha `$=2.5 up to 250 K; behaviors with $`\alpha `$=2 and 3 were reported earlier Finnemore01 ; Jung01 . The susceptibility (inset of Fig.1) is decreasing rapidly in the temperature range 39 - 34 K by 0.034 emu/g, whereas for lower temperatures it decreases only slowly to the value -0.04 emu/g. Magnetization measurements allow us to estimate a SC fraction of more than 50%.

For the optical investigations the sintered pellet of MgB2 was cut to a piece of 5$`\times `$5$`\times `$2 mm3 size and polished. We note that it had a smooth but not shiny surface, with a remaining roughness due to pores. The scanning electron microscopy analysis revealed some traces of oxygen on the sample surface which may correspond to a layer of the MgO on it.

Here we report the far-infrared reflectivity spectra measured in the grazing (80 degrees) incidence geometry using a FT-IR spectrometer. The grazing incidence geometry has been chosen because in this case the reflectivity is more sensitive to changes of optical conductivity due to opening of the superconducting gap. This experimental technique has been previously used to measure the $`s`$-wave gap of NbN and the $`d`$-wave gap of La1.85Sr0.15CuO4 Schutzmann97 . In addition, normal incidence spectra were collected in the same frequency range, which were qualitatively consistent with the grazing reflectance data reported in this paper.

In our measurements, we have observed that the copper block on which the sample is mounted in the cryostat and the sample surface had quite different temperatures. In order to get a precise surface temperature, an additional Pt thermoresistor was later attached to the sample face. It turned out that the temperature difference is strongly enhanced below 40 K, suggesting that the large thermal gradient is due to strong lowering of the sample thermal conductivity in the SC state. In particular, the lowest surface temperature we could achieve was only about 20 K, while the cryostat cold finger was at liquid helium temperature. This additional temperature calibration has been used in the analysis presented in this paper. However, we should note that the systematic error of this measurement on this rather thick sample is about 2-3 K.

Fig. 2 presents the intensity reflected from the sample at various temperatures normalized to the intensity reflected at 45 K, i. e., slightly above $`T_c`$. We clearly see a rise of the reflectivity ratio (RR) above 1, starting below $`T_c`$. This is seen also in Fig. 3, where the temperature dependences of the RR is plotted for fixed frequencies. The frequency below which the RR starts to increase above 1 shifts to higher values with decreasing $`T`$; for the lowest temperature $`T`$=20 K the increase starts at around $`\nu `$=70 cm-1. As expected, the RR spectra in the SC state reveal maxima, meaning that RR=1 at $`\nu `$=0.

The observed variation of the reflectivity is reminiscent of the behavior found in conventional superconductors and, to a certain extent, in the high-temperature cuprates. There, the condensation of carriers into pairs and the opening of the energy gap $`2\mathrm{\Delta }`$ in the density of states reveal themselves as a decrease of the conductivity starting around the corresponding frequency $`\nu `$=$`2\mathrm{\Delta }/hc`$, and as the response in the dielectric constant of $`ϵ1/\nu ^2`$, both leading to pronounced changes of reflectivity and a maximum in the RR at around $`\nu `$=$`2\mathrm{\Delta }/hc`$ Tinkham96 ; Timusk89 . Accordingly, we associate the changes we observe for MgB2 by entering the SC state with the response of condensed superconducting pairs, which is only in part obscured by the surface scattering. (We note that the rise in reflection for an essentially two dimensional system is not always an indication of the SC gap; for instance, the gap-like feature in the reflectivity spectrum of polycrystalline La1.85Sr0.15CuO4 Sherwin88  was later explained by the manifestation of the $`c`$-axis Josephson plasmon Tamasaki92 ).

A straight-forward assignment based on an isotropic $`s`$-wave gap would associate the gap energy $`2\mathrm{\Delta }`$ with the frequency where the RR reaches its maximum (around $`2530\mathrm{cm}^1`$, according to Fig.2), providing an unrealistically small estimate of $`2\mathrm{\Delta }`$ = 3 - 4 meV. On the other hand it has been anticipated that superconducting gap is not uniform and has different values for different four sheets of the Fermi surface Liu01 . With this interpretation the RR maximum should correspond to the minimum value of the anisotropic gap. This is also consistent with a large variation of $`2\mathrm{\Delta }`$ values obtained from tunneling, Raman, photoemission, and NMR techniques Kotegawa01 ; Rubio01 ; Sharoni01 ; Kaparetrov01 ; Schmidt01 ; Takahashi01 ; Chen01  ranging from 4 to 17 meV

In conclusion, we have investigated the optical properties of sintered MgB2 in the normal and superconducting state by measuring the ratio of the reflected intensity in the superconducting state to that in the normal state. An increase of the ratio in the superconducting state is reliably detected and interpreted as the signature of the superconducting gap with the lowest value of $`2\mathrm{\Delta }(0)`$$`34`$ meV.


format: text

text:
# Subsystem purity as an enforcer of entanglement

## Abstract

We show that entanglement can always arise in the interaction of an arbitrarily large system in any mixed state with a single qubit in a pure state. This small initial purity is enough to enforce entanglement even when the total entropy is close to maximum. We demonstrate this feature using the Jaynes-Cummings interaction of a two level atom in a pure state with a field in a thermal state at an arbitrarily high temperature. We find the time and temperature variation of a lower bound on the amount of entanglement produced and study the classical correlations quantified by the mutual information.

Entanglement is a key resource in quantum information processing . In practice, due to decoherence, it is very difficult to generate and maintain entanglement in the form of a pure state. With this motivation, there has been a lot of recent interest in understanding and quantifying entanglement of mixed states . As entanglement is an important resource, we need to investigate whether and how it can be generated in severe conditions of mixedness . Interestingly, if we had a qubit in a completely pure state $`|0_1`$ and another qubit in a completely mixed state $`(1/2)(|00|+|11|)_2`$ , we could let them interact through

$`|0_1|0_2`$ $``$ $`|0_1|0_2`$ (1)
$`|0_1|1_2`$ $``$ $`|\psi ^+_{12},`$ (2)
where $`|\psi ^+_{12}=\frac{1}{2}(|01_{12}+|10_{12})`$, to generate an entangled state. Surprisingly, this also holds for a pure qubit interacting with a fully mixed state of an arbitrary dimensional (even macroscopic) system. In this letter, we demonstrate this in the context of a pure state atom interacting with mixed state quantum field (an infinite dimensional system) through the Jaynes-Cummings model.

The Jaynes Cummings model (JCM)  is the simplest model that describes the interaction between light and matter. It consists of a two level atom interacting with a near-resonant quantized mode of the field. The atom and field become dynamically entangled by their interaction. This model provides direct evidence for the quantum mechanical nature of the electromagnetic field by predicting collapses and revivals of Rabi oscillations of the atom . These have been tested experimentally . The JCM has analytical solutions for arbitrary coupling constants and can easily be extended to include a variety of initial conditions , dissipation and damping , multiple atoms , multilevel atoms , multimode fields  and more elaborate interactions . For certain pure initial states of the atom and the cavity mode, the JCM produces entanglement which oscillates with time . In this paper we show that even if the cavity field is initially in a thermal state, entanglement still arises in course of the JCM interaction. The result holds irrespective of the temperature of the field. This goes counter to the folklore that entanglement, being a very quantum attribute, should automatically disappear at high temperatures. We give a lower bound on the entanglement produced as a function of the temperature of the cavity field and time. We also compare the entanglement, which quantifies the quantum part of the correlations, to the total correlations as quantified by the mutual information.

Entanglement is well understood for pure states of bipartite systems such as the JCM. The general state of a system involving two subsystems (say atom and field) can be written as a superposition of the products of individual states. If the state is pure it can always be written in the Schmidt form

$$|\mathrm{\Psi }=\underset{n}{}g_n|u_n^{}_a|v_n^{}_f$$

where $`|u_n^{}_a`$ and $`|u_n^{}_f`$ are orthonormal bases for the atomic and field subsystems respectively. The correlations of the two systems are then fully displayed. In this case (pure bipartite states), entanglement is quantified by the entropy of the reduced density matrix of either of the subsystems defined as $`S(\rho )=Tr(\rho _a\mathrm{log}(\rho _a))`$. This quantification of entanglement has been used in all earlier studies of entanglement in JCM with pure initial states . For the JCM with both pure atomic and cavity field states, the cavity field can be considered as an effective two level system .

Here we are interested in studying entanglement for an initial mixed state of the cavity field, as mixed states are the true representation of the state of the field at a finite temperature. Entanglement for mixed states is difficult to define. This is because we cannot easily define an analogue of the Schmidt decomposition for a general mixed state of a composite system. Such a mixed state can be expanded in terms of pure states in infinitely many different ways and it is not clear which, if any decomposition should be favored. Mixing two entangled pure states could result in a mixed state with entanglement much less than the average entanglement of the states mixed. Mixed state entanglement is thus a very different entity to either correlations or pure state entanglement. At least three different measures have been used to quantify entanglement for a mixed state. One of these measures, the relative entropy of entanglement , is defined for mixed state $`\rho _{af}`$ of a composite system (such as the atom-field system in the JCM) as

$$E_{re}(\rho _{af})=\underset{\sigma \mathrm{D}}{\mathrm{min}}(Tr(\rho _{af}(\mathrm{log}\rho _{af}\mathrm{log}\sigma _{af}))),$$

where $`\mathrm{D}`$ is the set of disentangled (separable) states of the system. A disentangled state can be written in the form $`_ip_i\rho _a^i\rho _f^i`$. This measure tells us how difficult it is to distinguish the given entangled state from its closest approximation to the set of separable states. The other measures of entanglement are associated with formation and distillation of entangled states. Consider the number, n, of copies of a non maximally entangled state $`\rho _{af}`$ that can be created by using only correlated local actions (i.e. through correlated actions on the field state and atomic state alone) on a number, m, of maximally entangled states. Entanglement of formation is the asymptotic conversion ratio, $`m/n`$ in the limit of infinitely many copies ,

$$E_f(\rho _{af})=\mathrm{min}\underset{i}{}p_iS(\rho _a^i)$$

where the minimum is taken over all the possible realizations of the state $`\rho _{af}=_ip_i|\mathrm{\Psi }_{af}^i\mathrm{\Psi }_{af}^i|`$. A measure stemming from the opposite process (distillation) is the entanglement of distillation. It is the asymptotic rate $`m/n`$ of converting $`m`$ copies of a non-maximally entangled state $`\rho _{af}`$ into $`n`$ copies of a maximally entangled state by means of correlated local actions. The entanglement of distillation is in general smaller than formation. All the different measures of entanglement are related to each other through the amount of available classical information about the decomposition of the state .

It is not easy to compute the value of entanglement from the measures. The entanglement of formation is the only measure for which an analytical method exists for calculating the entanglement, but this is specific to the case a state of $`2\times 2`$ systems. However, in our case there is a two level atom interacting with a cavity field, which is an infinite dimensional system. For such states, we can still give a lower bound on the entanglement from the known result in the case of $`2\times 2`$ systems. We first project the entire atom-field state onto a subspace equivalent to a $`2\times 2`$ system. We can then compute the entanglement of formation for each of the outcomes. This particular projection onto a $`2\times 2`$ system, as we will show, can be done by local actions alone. Being local, such an action cannot increase the entanglement on average . If we compute the average of the entanglement over all possible outcomes, the result will thus be a lower bound on the entanglement in the initial $`2\times \mathrm{}`$ state of the atom and the cavity field. We will also look at the total correlations between the atom and the mixed field as quantified by the mutual information ,

$$I=S(\rho _a)+S(\rho _f)S(\rho _{af})$$

where $`\rho _a`$ and $`\rho _f`$ are the reduced density matrices of the atom and field. This measure quantifies how much the correlated systems know about the state of each other. As the mutual information indicates the total correlations, it should be larger than the lower bound on entanglement we compute. We will compare the two quantities (i.e. the lower bound on the entanglement and the mutual information) to understand how much the purely quantum correlations contribute to the total correlations.

We consider the field in our JCM example initially to be in a thermal state at some temperature $`T`$ with the probability distribution $`P_n`$ for number states $`|n`$ being given by

$$P_n=\frac{1}{1+n}\left(\frac{n}{1+n}\right)^n$$
(3)
where $`n=\{e^{\beta \mathrm{}\omega }1\}^1`$ is the mean photon number, $`\beta =1/k_BT`$ and $`k_B`$ is Boltzmann’s constant, $`T`$ is the temperature, $`\omega `$ the frequency of the optical mode and $`\mathrm{}`$ the Planck’s constant. The two-level atom is initially taken to be in the excited state $`|e_a`$ (the ground state being $`|g_a`$). The JCM interaction between the atom and the field is given by

$$H_{\text{JCM}}=g(|eg|_aa_f+a_f^{}|ge|_a),$$
(4)
where $`a_f`$ and $`a_f^{}`$ are the annihilation and creation operators of the field mode respectively. The joint density matrix for the atom-field system evolves with time as

$$\rho _{af}=\underset{n=0}{\overset{\mathrm{}}{}}P_n\rho _n,$$
(5)
where

$`\rho _n`$ $`=`$ $`\mathrm{cos}^2\left({\displaystyle \frac{\mathrm{\Omega }_nt}{2}}\right)|e,ne,n|_{af}`$ (6)
$``$ $`i\mathrm{cos}\left({\displaystyle \frac{\mathrm{\Omega }_nt}{2}}\right)\mathrm{sin}\left({\displaystyle \frac{\mathrm{\Omega }_nt}{2}}\right)|e,ng,n+1|_{af}`$ (7)
$`+`$ $`i\mathrm{cos}\left({\displaystyle \frac{\mathrm{\Omega }_nt}{2}}\right)\mathrm{sin}\left({\displaystyle \frac{\mathrm{\Omega }_nt}{2}}\right)|g,n+1e,n|_{af}`$ (8)
$`+`$ $`\mathrm{sin}^2\left({\displaystyle \frac{\mathrm{\Omega }_nt}{2}}\right)|g,n+1g,n+1|_{af}`$ (9)
where $`\mathrm{\Omega }_n=2g\sqrt{n+1}`$ is the Rabi frequency.

We cannot exactly compute the total atom-field entanglement in the mixed state given by the above Eq.(5). However, we can obtain an estimate a lower bound on the entanglement using the projection to a $`2\times 2`$ subspace discussed earlier. To see this, consider the field state being projected into the subspace spanned by $`|n_f`$ and $`|n+1_f`$. As this is a local action on the field, it cannot increase the entanglement. The resulting state is

$$\rho _{af}^n=\left(\begin{array}{cccc}P_{n1}S_{n1}^2& 0& 0& 0\\ 0& P_nS_n^2& P_niC_nS_n& 0\\ 0& P_niC_nS_n& P_nC_n^2& 0\\ 0& 0& 0& P_{n+1}C_{n+1}^2\end{array}\right)$$

where

$$C_n=\mathrm{cos}\left(\frac{\mathrm{\Omega }_nt}{2}\right),S_n=\mathrm{sin}\left(\frac{\mathrm{\Omega }_nt}{2}\right).$$
(10)
Before proceeding to the evaluation of entanglement based on the above formula, we would check the separability of the above state. To prove the inseparability of this matrix we compute the eigenvalues of the partially transposed matrix. If one of the eigenvalues is negative then $`\rho _{af}^n`$ is inseparable . The existence of the negative eigenvalue reduces to the condition

$`(P_nC_nS_n)^2>P_{n1}P_{n+1}(C_{n+1}S_{n1})^2.`$ (11)
Substituting $`P_n,P_{n+1}`$ and $`P_{n1}`$ from Eq.(3) in the above expression we obtain the condition

$`\mathrm{\Lambda }_n=(C_nS_n)^2(C_{n+1}S_{n1})^2>0,`$ (12)
which is independent of $`n`$ (i.e. of the temperature). We plot the expression $`\mathrm{\Lambda }_n`$ with time $`t`$ in Fig.1 for three values of $`n`$.

From Fig.1, we see that the expression $`\mathrm{\Lambda }_n`$ is positive for some intervals of time for each $`n`$ implying $`\rho _{af}^n`$ is entangled during those intervals of time. The entire atom-cavity state $`\rho _{af}`$ is thus also entangled in those intervals of time (otherwise, no local projection would have given an entangled outcome). Entanglement thus arises due to JCM interaction of an atom and a field in a thermal state irrespective of the temperature of the field. Note that this is a rigorous result because even if $`\mathrm{\Lambda }`$ is negative for just one specific value of $`(n,t)`$, then a local projection at that $`t`$ will result in the entangled outcome $`\rho _{af}^n`$ with a finite probability. This would mean that prior to the local measurement, entanglement was present in the atom-field state at that time.

Based on the plots in Fig.1, we heuristically justify the following conjecture: Entanglement is present at all times except at $`t=0`$. The conjecture relies on the observed behaviour that for higher values of $`n`$, oscillations of $`\mathrm{\Lambda }_n`$ (separability) with $`t`$ are faster. Consider an arbitrarily small interval of time $`\delta t<ϵ`$. By going to sufficiently high $`n`$, one can always find that $`\mathrm{\Lambda }_n`$ has a period smaller than $`\delta t`$ and thereby $`\rho _{af}^n`$ is entangled for an interval of time within $`\delta t`$. We can choose the time interval $`\delta t`$ smaller and smaller, so that eventually there is entanglement at all instants of time.

We now briefly comment on the alternative scenario where the atom is in the initial thermal state

$$\lambda |ee|+(1\lambda )|gg|,$$
(13)
where $`\lambda /(1\lambda )=\mathrm{exp}(\mathrm{\Delta }E/kT)`$ ($`\mathrm{\Delta }E`$ being the energy difference between $`|e`$ and $`|g`$ and $`T`$ the temperature). The field is assumed to start in a pure Fock state $`|n`$. As the JCM is completely symmetric between the atomic and field operators, we can exchange the field and the atom states and have the same entanglement. This new scenario is equivalent to just two levels of the field involved in interaction with a pure atom. As shown earlier this also leads to entanglement at all temperatures of the atom.

The case when both the atom and the field are in thermal states is much more complex. In this case we expect that there is a cut-off temperature above which there is no entanglement. At infinite temperature, of course, all entanglement disappears as the total state is just proportional to the identity. However, our projection method can only lead to lower bounds as the failure of our method to produce entanglement by local projections does not imply that entanglement does not exist in the original mixed state. So, until there an operational necessary and sufficient condition for separability of $`2\times N`$ density matrices, we cannot fully address this case.

As mentioned earlier, if we evaluate the entanglement of $`\rho _{af}^n`$ and take an average over all possible values of $`n`$ (weighted by the probability of obtaining $`\rho _{af}^n`$), we should get a lower bound on the entanglement of the state $`\rho _{af}`$. This holds because $`\rho _{af}^n`$ is the result of a local projection on $`\rho _{af}`$. We have computed this lower bound on entanglement and the mutual information for initial thermal states with various mean photon numbers (e.g. a plot for the mean photon number $`n=10`$ is shown in Fig. 2). We find, as we proved earlier, that JCM generates entanglement at any temperature of the cavity field, no matter how high. Moreover, though it starts from zero, it never completely vanishes. This supports our earlier conjecture that entanglement is present at all times except at $`t=0`$. Mutual information, on the other hand, quantifies total correlations and is therefore greater or equal to entanglement as seen in Fig. 2. Classical correlations are, of course, always seen to be present and this means that the total state is never of the product form at $`t>0`$ due to the JCM interaction. The mutual information is a good measure of classical correlations when entanglement is absent, since then total correlations are equal to classical correlations.

We have shown that entanglement can be generated from a very small amount of purity existing in an overall very mixed state. In principle, a $`2\times N`$ system can be entangled even if the total entropy is as high as $`\mathrm{log}N`$ (the maximum being $`1+\mathrm{log}N`$). This is an important result among attempts to relate mixedness and entanglement . This means in principle, one could entangle a microscopic system in a pure state with a macroscopic system in a thermal state, as suggested by the results in . Our results should have implications for quantum computation with mixed states . This kind of entanglement is probably involved in the functioning of Shor’s algorithm with only one pure qubit . Further studies could involve other natural interactions. For example, one could study how the creation of entanglement is affected by the off resonant interaction of the atom and the field. It is also an interesting and open question how efficiently one can entangle two partially mixed systems in general.

This research has been partly supported by the European Union, The UK Engineering and Physical Sciences Research Council and Hewlett-Packard. I. F.-G. would like to thank Consejo Nacional de Ciencia y Tecnologia (Mexico) Grant no. 115569/135963 for financial support.

2


text:
# Period - magnitude relationships in BVIJHK-Bands for fundamental mode and first overtone Cepheids

## 1 Introduction

Over the last five years, a wealth of data for Cepheids was collected by several observational projects (HST key project, HIPPARCOS, microlensing experiments EROS, MACHO, OGLE, etc…). This fantastic amount of data allows unprecedented opportunities to understand the properties of these variable stars. One of the relevance of Cepheids, and the main goal of the HST key project (cf. Freedman 2000, and references therein), is the determination of local extragalactic distance scales and consequently of the Hubble constant $`H_0`$. The period - luminosity relationship (PL) of Cepheids is the conerstone of such purpose. An important contribution to the collection of data is provided by microlensing search projects toward the Magellanic Clouds (EROS: Renault et al. 1996; OGLE: Udalski, Kubiak and Szymanski 1997; MACHO: Alcock et al. 1997; and references therein). The recent release to the public domain of their Cepheid Catalog (EROS: Afonso et al. 2000; OGLE: Udalski et al. 1999a,b) provides large, homogeneous samples of high quality data which in principle allow a deep analysis of the period-luminosity relationships. In particular, a key question related to the PL relationship is its dependence on metallicity. This is of fundamental importance for extragalactic distance scales, usually based on an universal PL relationship. A controversy still exists regarding metallicity effects, both on the observational and theoretical viewpoints (see Alibert et al. 1999 and references therein). Another notorious debate concerns the distance to the LMC, which is an important source of uncertainty in the determination of $`H_0`$ (see Freedman 2000). Although PL relationships used for distance estimate are based on fundamental (F) pulsators, the afore-mentioned catalogs provide also statistically significant samples of first overtone (1H) Cepheids. Because of the distinction of F and 1H pulsator sequences in period - magnitude diagrams, 1H provides additional constraints on theoretical models and on distance determinations based on PL relationships.

In a recent paper, Alibert et al. (1999) performed self-consistent stellar evolution and linear stability analysis calculations for Cepheids. Their analysis was essentially devoted to fundamental mode pulsators and their period-magnitude-color relationships as a function of metallicity. The analysis of their results for first overtone pulsators and the comparison with observations were hampered at that time by the lack of significant samples of data. The aim of the present paper is to complete the work of Alibert et al. (1999) for first overtone pulsators and test their results against the observational constraints provided by both F and 1H pulsators. In section 2, we derive PL relationships for 1H and in §3 we compare the results to observations. In §4, we briefly discuss the uncertainties of our period - magnitude relationships and consequences on distance determinations.

## 2 Period - magnitude relationships for first overtone pulsators

The description of the stellar evolution and pulsation calculations are given in details in Alibert et al. (1999). We recall the main ingredients: (i) evolutionary models for Cepheids are constructed with the Lyon evolutionary code from 3 to 12 $`M_{}`$. We consider various initial compositions $`(Z,Y)=(0.02,0.28)`$, $`(0.01,0.25)`$ and $`(0.004,0.25)`$ 111$`Z`$ is the metal mass fraction and $`Y`$ the helium mass fraction, representative of respectively the Galactic, the Large Magellanic (LMC) and the Small Magellanic (SMC) clouds environments. Evolutionary calculations do not include core overshooting. (ii) A linear non-adiabatic stability analysis is performed on the complete evolutionary models along the evolutionary tracks. This provides consistent mass-age-period-luminosity relations. (iii) Static atmosphere models and their corresponding synthetic spectra are calculated for the same compositions used in (i), providing magnitudes and colors for a given ($`Z`$, $`M`$, $`L`$, $`T_{\mathrm{eff}}`$).

### 2.1 Period - luminosity and Period - magnitude relationships in different filters.

In order to derive statistical PL relationships, a mean position is assigned to each mass, according to the time spent in different locations in the instability strip (IS) and a linear least-square fit to these points is derived. As shown in Alibert et al. (1999), the mean position of a given stellar mass, accounting for its evolutionary time, is roughly located in the middle of the IS.

The minimum masses $`m_{\mathrm{min}}`$ undergoing a blue loop in the 1H instability strip are 3.25 $`M_{}`$, 4 $`M_{}`$ and 5 $`M_{}`$ for respectively $`Z`$ = 0.004, 0.01 and 0.02. In order to avoid biases due to the change of slope predicted near these minimum masses (see Alibert et al. 1999, their §3.4), we exclude in the present analysis $`m_{\mathrm{min}}`$ to derive the mean PL relationships. Because of the reduction of the blue loop extension toward lower masses, the slope of the PL relationships for 1H can be affected if the fit is derived down to $`m_{\mathrm{min}}`$, yielding steeper relations. Bauer et al. (1999) observed such a change of slope for fundamental pulsators in the SMC, but did not observe it for 1H, although predicted by the models (see the discussion by Alibert et al. 1999, §4.2.3). The relationships given in Tables 1 and 2 are then derived from 3.5 $`M_{}`$ for $`Z`$ = 0.004, 4.25 $`M_{}`$ for $`Z`$ = 0.01 and 5.5 $`M_{}`$ for $`Z`$ = 0.02. Note that Alibert et al. (1999) derived preliminary PL relationships for 1H based on the whole range of masses unstable in the IS (including $`m_{\mathrm{min}}`$) and obtained slightly different results than in Table 1.

The period - magnitude relationships in different filters are given in Table 2 for F and 1H modes. $`BVRI`$ magnitudes are based on the Johnson-Cousins system (Bessell 1990) and $`JHK`$ magnitudes are defined in the CIT system (Leggett 1992). Table 2 displays also the reddening-free Wesenheit index $`W_\mathrm{I}=I\mathrm{\hspace{0.17em}1.55}(VI)`$ (Madore and Freedman 1991), with the coefficient 1.55 ($`=A_\mathrm{I}/E(VI)`$) resulting from standard interstellar extinction curves (e.g., Schlegel et al. 1998).

### 2.2 Effects of metallicity

Figure 1 displays the mean $`\mathrm{log}P_1`$ \- $`\mathrm{log}L/L_{}`$ relationships (Fig. 1a) and the blue and red edges in a $`\mathrm{log}P_1`$ \- $`T_{\mathrm{eff}}`$ diagram (Fig. 1b) as a function of $`Z`$. Excluding the minimum mass undergoing a blue loop in the IS, which decreases with $`Z`$, there is no noticeable effect of $`Z`$ on the location and width of the IS in Fig. 1b. The periods corresponding to $`m_{\mathrm{min}}`$ in Fig. 1 are $`\mathrm{log}P_1`$ 0.3, 0.1 and 0 for respectively $`Z`$ = 0.02, 0.01 and 0.004. We note however a small effect of $`Z`$ on the PL relationship, with $`\mathrm{log}L`$ increasing by 0.04 for $`Z`$ increasing from 0.004 to 0.02 at a given $`P_1`$.

In the following, our analysis is restricted to period $`\mathrm{log}P_1`$ 0.3, since $`\mathrm{log}P_1`$ = 0.3 is the minimum period for $`Z`$=0.02 1H pulsators undergoing a blue loop (see Fig. 1). An inspection of the relationships given in Table 1 shows that the slopes hardly depend on $`Z`$, but the zero points are affected by metallicity. Although small, the increase of $`L`$ with $`Z`$ translates in brighter magnitudes at a given $`P_1`$, with the largest effects in the $`VIJ`$ bands and the maximum effect in the $`I`$ band. We note that $`M_\mathrm{V}`$, $`M_\mathrm{I}`$ and $`M_\mathrm{J}`$ are brighter by 0.1 - 0.15 mag, $`M_\mathrm{H}`$ and $`M_\mathrm{K}`$ are brighter by 0.8 - 0.13 mag at a given $`P_1`$ when $`Z`$ increases from 0.004 to 0.02. Between $`Z`$=0.004 and $`Z`$ = 0.01 the effect is $``$ 0.1 mag in the $`VIJHK`$ bands. Interestingly enough, the effect is the smallest in the $`B`$ band, with a difference of $`M_\mathrm{B}`$ less than 0.1 mag between $`Z`$ = 0.004 and $`Z`$ = 0.02. Opposite trends are found for fundamental mode PL relationships which show the largest $`Z`$ effect in the $`B`$ band, but small effects (less than 0.1 mag) at longer wavelengths (see Alibert et al. 1999, §5).

This effect of $`Z`$ on the PL relationship for 1H is related to the behavior of the period ratio $`P_1`$/$`P_0`$ as a function of $`Z`$. This ratio decreases as $`Z`$ increases (cf. Baraffe et al. 1998): we find for $`Z`$ = 0.004, $`P_1`$/$`P_0`$ $``$ 0.73 and for $`Z`$ = 0.02, $`P_1`$/$`P_0`$ $``$ 0.7. Since for a given $`L`$, $`P_0`$ is essentially independent of $`Z`$, $`P_1`$ decreases because of decreasing $`P_1`$/$`P_0`$ ratio for increasing $`Z`$. Conversely, for a given $`P_1`$, $`L`$ increases with $`Z`$.

The behavior of $`P_1`$/$`P_0`$ with $`Z`$ can be understood by analysing the regions in the star which contribute to the period of each mode. Such regions can be determined by the Epstein weight functions which are obtained by expressing the eigen frequency of each mode by an integral over the whole star (Epstein 1950; Cox 1980). As already noted by Epstein (1950), the greatest contribution to the fundamental mode period is located between $`r/R`$ 0.7 and $`r/R`$ 0.9. For the first overtone, two regions are important: a first inner region located between $`r/R`$ 0.4 and $`r/R`$ 0.7 and an outer zone with $`r/R>0.85`$. These regions are characteristic of the mass range of interest (5-10 $`M_{}`$). The inner region for 1H includes the opacity peak due to metals at $`T\mathrm{2\hspace{0.17em}10}^5`$ K and the outer region covers the H-He ionization zone. In the first region covering the metal opacity peak, an increase of the opacity $`\kappa `$, related to an increase of the metallicity, yields an increase of $`T`$ or a decrease of the density $`\rho `$ at fixed $`r`$ (or $`P`$). Note that this region is dominantly radiative and the temperature gradient depends directly on $`Z`$ since the radiative gradient is proportional to $`\kappa `$. Consequently, in this region, the sound speed $`C_\mathrm{s}(P/\rho )^{1/2}T^{1/2}`$ is more sensitive to the metallicity than in other outer regions and increases with $`Z`$. The integral of the sound travel time $`\tau _\mathrm{s}=dr/C_\mathrm{s}`$ through this zone is thus more sensitive to the metallicity, decreasing as $`Z`$ increases, than in the other regions located at $`r/R>0.7`$. The period ratio $`P_1`$/$`P_0`$ is related to the ratio of $`\tau _\mathrm{s}`$ through the regions relevant for 1H to $`\tau _\mathrm{s}`$ in the zones contributing to F. Because of the larger sensitivity to $`Z`$ of $`\tau _\mathrm{s}`$ for 1H, this ratio decreases as $`Z`$ increases. We checked our arguments by calculating the integral of $`\tau _\mathrm{s}`$ through the different regions of interest for different cases of $`M`$, $`L`$ and $`T_{\mathrm{eff}}`$ covering the instability strip and for different values of $`Z`$.

The behavior of $`L`$ as a function of $`P_1`$ with $`Z`$ explains the trends previously mentioned of the fluxes in different filters. In particular, it explains the low sensitivity of the $`\mathrm{log}P_1`$ \- $`M_\mathrm{B}`$ relationship to $`Z`$ because of the compensating effects of (i) higher $`L`$ with increasing $`Z`$ at a given $`P_1`$ and (ii) lower $`B`$-flux for a given $`L`$ and $`T_{\mathrm{eff}}`$ as $`Z`$ increases, because of more metallic line absorption in the $`B`$ bandpass (see §5 of Alibert et al. 1999).

Finally, because of the maximum effect in the $`I`$-bandpass, we note a non negligible effect of $`Z`$ on the index $`W_\mathrm{I}`$, which differs by 0.14 mag at long $`P_1`$ up to 0.18 mag at shorter $`P_1`$ between Z=0.004 and 0.02. When $`Z`$ increases from 0.004 to 0.01, $`W_\mathrm{I}`$ gets brighter by $``$ 0.13 mag.

## 3 Period - magnitude diagrams: comparison with observations

Comparison with observations are made with the EROS2 and OGLE2 data for the MC. EROS2 (Afonso et al. 2000) reports 239 1H Cepheids in the SMC and 113 in the LMC. OGLE2 provides $``$ 800 1H pulsators in the SMC and $``$ 500 in the LMC.

### 3.1 SMC Cepheids: $`VI`$-Bands

Figure 2 displays the comparison between models with metallicity $`Z`$ = 0.004 and the OGLE2 data. EROS2 data are not shown for the sake of clarity, but are comparable to the OGLE2 data. We adopt the dereddened data provided by the OGLE2 catalog, with different reddening corrections for each fields (see Udalski et al. 1999b for details) and the same distance modulus $`\mu _0=18.94`$ as used in Alibert et al. (1999). This value, suggested by Laney and Stobie (1994), yields general agreement between models and various observations for fundamental pulsators (cf. Alibert et al. 1999). Unstable models during the first crossing phase (large open circles) and the core He burning blue loop (filled circles) are indicated (see Alibert et al. 1999 for details). As shown in Fig. 2, a general good agreement is found for 1H pulsators. The fundamental mode sequence is indicated by the mean PL relationship, and compared to the relationship of Udalski et al. (1999c), recently revised by Udalski (2000). Note that these authors truncate their SMC sample at $`\mathrm{log}P<\mathrm{\hspace{0.17em}0.4}`$, to avoid biases at short periods. The agreement between predicted and observed F relationships is thus excellent for $`\mathrm{log}P>\mathrm{\hspace{0.17em}0.4}`$. As preliminary noted by Alibert et al. (1999) and confirmed by the data, observations for 1H do not extend above $`\mathrm{log}P\mathrm{\hspace{0.17em}0.6}`$, which corresponds to $`m\mathrm{\hspace{0.17em}5}M_{}`$ with standard (no overshooting) evolutionary models. Note that models including overshooting would yield a lower mass for such period. Our linear stability analysis however finds unstable 1H modes (as well as F modes) for masses $`m>5M_{}`$. Only non-linear calculations can determine the dominant mode of pulsation. In the same vein, models predict that 1H pulsators with $`\mathrm{log}P<\mathrm{\hspace{0.17em}\hspace{0.17em}0}`$ are in the first crossing phase. The significant number of such low period pulsators provides another strong constraint on non-linear calculations which could test the following scenarios: (1) models on the first crossing, characterised by a lower luminosity than during the He core burning blue loop for a given mass, should favor 1H as the dominant mode, and (2) models with $`m>\mathrm{\hspace{0.17em}\hspace{0.17em}5}M_{}`$ on the blue loop should oscillate predominantly in the fundamental mode.

We note however that in case (1), one may expect a change in the number density of objects with $`\mathrm{log}P<\mathrm{\hspace{0.17em}\hspace{0.17em}0}`$, given the much faster evolutionary timescale on the first crossing compared to the blue loop phase. Such a change is not displayed by current observations. A detailed statistical analysis taking into account evolutionary timescales, mass functions and eventual observational biases is required to investigate this point.

### 3.2 LMC Cepheids: $`VI`$-Bands

Comparison between models with Z=0.01 and observations in the LMC is shown in Fig. 3. Data are from EROS2 (Afonso et al. 2000) and OGLE2 (Udalski et al. 1999a). We adopt the same reddening corrections as in Udalski et al. (1999a) for the OGLE2 sample, which varies from field to field. For the EROS2 data, we adopt the same reddening correction $`E(BV)`$=0.10 as used in Alibert et al. (1999). We note that adopting a constant $`E(BV)`$=0.10 for the OGLE2 data has no noticeable effect when comparing models and data in Fig. 3. We adopt the same distance modulus 18.50 as in Alibert et al. (1999). The theoretical PL relations for F and 1H modes are compared to the observed relations derived by Udalski et al. (1999c) on the LMC sample limited to $`\mathrm{log}P>\mathrm{\hspace{0.17em}0.4}`$, and recently revised by Udalsky (2000). The same agreement as for the SMC data is found and the models reproduce the width and location of the instability strip satisfactorily. We however note discrepancies between models and observations for the mean $`P`$ \- $`M_\mathrm{V}`$ relationship for both F and 1H modes. This discrepancy appears also, but to a lesser extent, for the SMC (cf. Fig. 2) and is better illustrated in $`P`$ \- $`W_\mathrm{I}`$ diagrams (see next section).

As illustrated in Fig. 3, no 1H pulsators are observed above $`\mathrm{log}P0.8`$, corresponding to $`m>\mathrm{\hspace{0.17em}7}M_{}`$. Note that the lack of 1H pulsators at long periods could also be interpreted in terms of very small amplitudes below the level of detection of OGLE or EROS. However, there is no obvious decrease of the amplitudes of 1H pulsators in the LMC, as well as in the SMC (see Afonso et al. 1999), as the period increases which could support this interpretation. Below $`\mathrm{log}P\mathrm{\hspace{0.17em}0.1}`$, observations can be explained by first crossing models. We therefore find the same type of constraints as derived in the previous section from the SMC data on non-linear calculations for the determination of the dominant pulsating mode. Such comparisons suggest that the afore-mentioned properties in scenarios (1) and (2) (see §3.1) are intrinsic to 1H pulsators for $`Z`$ varying from 0.004 to 0.01. More observations at higher metallicities are currently required to determine if these properties apply also to Galactic Cepheids.

### 3.3 $`P`$ \- $`W_\mathrm{I}`$ diagram

Figures 4 and 5 display observations and models for both F and 1H in a $`\mathrm{log}P`$ \- $`W_\mathrm{I}`$ diagram, for respectively the SMC and the LMC. The Wesenheit index $`W_\mathrm{I}`$ is an useful quantity since it is rather insensitive to the reddening and can in principle remove part of the scatter due to differential reddening. For the LMC, the $`Z`$ = 0.01 models are also compared to observations from Gieren et al. (1998) for F Cepheids. Note that Alibert et al. (1999) found good agreement in different optical and near-infrared $`P`$ \- magnitude diagrams with the latter sample of data. As expected, the scatter of data in such diagram is smaller than in the $`V`$ and $`I`$ bands. A first inspection of Figs. 4 and 5 show a general agreement between models and observations for both SMC and LMC: the models corresponding to different masses cover the observed location of both F and 1H Cepheids. We note that the models are in better agreement with the observed width of the IS in the LMC, whereas the data in the SMC show a larger scatter. This is consistent with the recent analysis by Groenewegen (2000) who interprets the larger dispersion in the PL relation in $`W_\mathrm{I}`$ for the SMC compared to the LMC in terms of a larger intrinsic depth of the former Cloud.

For both the SMC and the LMC, the predicted $`\mathrm{log}P`$ \- $`W_\mathrm{I}`$ mean relationships are shifted compared to the observed relationships. Since the shift between observed and predicted relationships is almost constant on the whole range of $`P`$ and is the same between F and 1H relations, we find that an increase of the distance moduli $`\mu _0`$ by the same amount $``$ 0.15 - 0.2 mag for both clouds removes easily this discrepancy. This is illustrated in the insets of Figs. 4 and 5. Note that such an increase of $`\mu _0`$ does not alter the agreement found between observations and models in the $`I`$-band (see Fig. 2 and Fig. 3). It however increases the discrepancy in the $`V`$-band mentioned in the previous section. It is difficult to analyse the reason of such inconsistency. The choice of the reddening and of the extinction curve coefficients, affecting mostly the $`V`$-band, could be the reason of this problem. This emphasizes the large uncertainties intrinsic to observations in this filter and suggests to use rather $`W_\mathrm{I}`$ to derive a distance modulus. Note however that, independently of reddening and distance modulus, the models do predict the correct location of the observed 1H sequence relative to the F sequence, for both SMC and LMC.

We also note that the models yield a difference of distance moduli between SMC and LMC of $`0.45`$, in good agreement with previous determinations based on Cepheids (Laney and Stobie 1994; Udalski et al. 1999c), RR Lyr or Red Clump stars (cf. Udalski et al. 1999c; Udalski 2000). Finally, Caputo et al. (2000) also derived theoretical $`\mathrm{log}P`$ \- $`W_\mathrm{I}`$ relationships for the fundamental mode which are shallower than the present relationships. Alibert et al. (1999) already mentioned and discussed this difference for fundamental mode period - magnitude relationships. As illustrated in Figs. 4 and 5, the Caputo et al. (2000) $`\mathrm{log}P`$ \- $`W_\mathrm{I}`$ relationships are also shallower compared to the observed data in both SMC and LMC. For the distance moduli adopted in the insets of Figs. 4 and 5, the F sequence of Caputo et al. (2000) reproduce the F data at long period, but merges into the observed 1H pulsators at shorter periods. A variation of the distance modulus cannot clearly solve such discrepancy. This highlights the constraint on models and distance moduli provided by the combination of F and 1H observed sequences.

## 4 Uncertainties on distance determination

Alibert et al. (1999) analysed the effects of uncertainties inherent to stellar evolution models: convection treatment, overshooting, mass loss, initial helium abundance, etc… The main conclusion is that such uncertainties barely affect the PL relationships and the conclusion on metallicity effects. We however recall that such uncertainties can affect the comparison between models and observations if the mass is used as another constraint.

The main uncertainty in our calculations is due to the neglect of convection - pulsation coupling in our linear stability analysis. Convection is frozen in, which means that the perturbation $`\delta F_{\mathrm{conv}}`$ of the convective flux is neglected in the linearized energy equation. An arbitrary criterion is used to define a red edge for the instability strip (see Alibert et al. 1999 for details), since such approximation cannot yield naturally a red edge.

A time dependent non local theory of convection is required to take into account the effect of convection on pulsation. Such a theory is lacking and the current recipes include several free parameters (cf. Yecko et al. 1998 and references therein). Despite the lack of a robust theory, it has now become clear that convective energy transport is essential to describe the observed properties of pulsating stars such as light curves or Fourier coefficients (Yecko et al. 1998; Feuchtinger 1999; Feuchtinger et al. 2000).

In a systematic analysis based on linear stability analysis of Cepheid models, Yecko et al. (1998) have shown the high sensitivity of the position of the IS to these free parameters, which need to be calibrated according to astronomical observations. But as emphasized by Feuchtinger et al. (2000) even a correct description of the instability strip requires hydrodynamical calculations, since the red edge is determined by nonlinear effects. Even for the blue edge, Feuchtinger et al. (2000) report a shift toward higher $`T_{\mathrm{eff}}`$ by about 350 K for the fundamental mode when convection is included in the pulsation calculation.

In the present paper, we derive a rough estimate of the sensitivity of the period - magnitude relationships derived under the afore-mentioned approximations, and the effects on distance determination. A more detailed analysis and the quantification of uncertainties resulting from convection is addressed in a forthcoming paper (Alibert & Baraffe 2001). Since our results yield a general good agreement with the width of observed instability strips (see Alibert et al. 1999), a first estimate of such uncertainties can easily be derived by shifting arbitrarily both the blue and red edges of the instability strip toward cooler or hotter $`T_{\mathrm{eff}}`$, keeping the width unmodified. In the following, we focus on fundamental PL relationships, since they are mostly used for distance determinations. Inspired by the afore-mentioned results by Feuchtinger et al. (2000), we adopt a shift in $`T_{\mathrm{eff}}`$ of 350 K. Based on the models appropriate for the LMC ($`Z`$ = 0.01), a shift by 350 K toward hotter (cooler) $`T_{\mathrm{eff}}`$ yields brighter (fainter) magnitudes at a given $`P`$. The effect decreases toward near-IR filters: the resulting variation of the magnitude in the $`B`$-band reaches up to 0.6 mag, it does not exceed 0.3 mag in the $`I`$-band and remains below 0.2 mag at longer wavelengths. This is expected, given the decreasing sensitivity to $`T_{\mathrm{eff}}`$ of the flux toward longer wavelengths. In the $`V`$-band, $`M_\mathrm{V}`$ varies from 0.3 mag at short $`P`$ ($`\mathrm{log}P=0.5`$) up to 0.45 mag at longer $`P`$ ($`\mathrm{log}P=1.7`$). In the $`K`$ band, the effect is less than 0.1 mag on the whole range of periods. This simple test highlights the high sensitivity of distance determinations based on theoretical $`P`$ \- $`M_\mathrm{V}`$ relationships, since a variation of 350 K on the location of the blue or red edge can easily result from uncertainties due to convection (see Yecko et al. 1999; Feuchtinger et al. 2000; Alibert and Baraffe 2001).

Finally, since our simple test clearly shows that near-IR period - magnitude relationships are more reliable, it is worth to derive a distance modulus for LMC based on near-IR data. Recent observations in the near-IR (Madore, 2000 priv. comm.; Persson et al. 2000) complete data already available from Laney and Stobie (1994) and Gieren et al. (1998). The $`P`$ \- $`M_\mathrm{K}`$ relationship (Table 2) yields a distance modulus for LMC of 18.60- 18.70, as illustrated in Fig. 6. Interestingly enough, this distance modulus is in good agreement with the value determined in §3.3 based on $`W_\mathrm{I}`$ and the combination of F and 1H observed sequences.

## 5 Conclusion

We have derived period - magnitude relationships for 1H pulsators in different filters. In contrast to fundamental mode, we find a small effect of metallicity in PL relationships for first overtone pulsators. This effect is due to the dependence on metallicity of the ratio $`P_1`$/$`P_0`$. Our models reproduce the location of both F and 1H observed sequences of OGLE2 in the SMC and the LMC. Using the reddening free index $`W_\mathrm{I}`$, models and observations are in good agreement for a LMC distance modulus $`\mu _0`$ = 18.65 - 18.70. We note however an inconsistency since this value yields a significant discrepancy between predictions and data in the $`V`$-band. Such inconsistency may illustrate problems in the choice of the extinction coefficients and the reddening, which are important for the comparison between models and data in the $`V`$-band but not crucial for $`W_\mathrm{I}`$.

We show that an arbitrary shift of 350 K in $`T_{\mathrm{eff}}`$ of the location of the instability strip yields up to 0.45 mag effect on $`M_\mathrm{V}`$ at a given $`P`$. Toward near-infrared wavelengths, the effect is smaller and less than 0.1 mag in the $`K`$-band. Using recent Cepheid data in the $`K`$-band for fundamental pulsators in the LMC, the models predict a distance modulus for LMC $`\mu _0`$ = 18.60 - 18.70, in good agreement with the predictions based on $`W_I`$ and combined 1H and F data.

Although Cepheid observations in the $`V`$-band are widely used for distance determination, we stress the high uncertainties inherent to observations in this filter (reddening correction, intrinsic dispersion in the instability strip, etc…) and the high sensitivity of $`P`$ \- $`M_\mathrm{V}`$ relationships to theoretical uncertainties. We also emphasize that even if the coupling between convection and pulsation provides a large source of uncertainty in the present work and in all current theoretical calculations, predictions in the $`K`$-band are expected to be much less affected by such uncertainty.

###### Acknowledgements.

We are indebted to Andrzej Udalski for providing data and results prior to publication. We are grateful to Barry Madore for valuable discussions and for making available to us near-IR data for the LMC. Many thanks to our referee, A. Gautschy, who contributed to the improvement of the manuscript and provided excellent ideas. The calculations were performed using facilities at Centre d’Etudes Nucléaires de Grenoble.


text:
# On the Dominance of Trivial Knots among SAPs on a Cubic Lattice

## 1 Introduction

The Self-Avoiding Polygon (SAP) with fixed topology gives a simplified model of real ring polymers in solution that have a topological constraint as well as excluded volume. Throughout the time evolution, a circular polymer keeps the same knot which is given to it when it is made; it does not change its topology under any thermal fluctuations since no crossing through itself is allowed. On the other hand, SAP corresponds to the special case of the Self-Avoiding Walk (SAW) that returns to the origin. If we construct a set of SAPs, then their topological states may contain several different knots. Therefore, it is not trivial how to realize the topological constraint on an SAP. One possible method for assigning the topological constraint on an SAP is that after generating a large number of SAPs we select only such SAPs that have the same given knot. By this method, we can assign the topological constraint on any model of SAP. In this context, the probability that a given SAP has a fixed knot plays a central role, and we call it the knotting probability of the model of SAP for the knot. Among many different models, the SAP on the cubic lattice with fixed topology is one of the most fundamental models of SAP. It has an advantage that the definition is very simple. We expect that the model should be suitable for general and mathematical study. In fact, several rigorous results on knotting probability have been derived for the SAP on the cubic lattice . Thus, the main motivation of the present research is to characterize the knotting probability of the SAP on the cubic lattice through numerical simulations.

Let us discuss some previous numerical results on knotting probabilities of SAP . For a model of SAP (or random polygon) with $`N`$ steps, we denote by $`P_K(N)`$ the knotting probability of the model for knot $`K`$. For the Gaussian model of random polygon and the rod-bead model of SAP, knotting probabilities were evaluated through numerical simulations for the trivial knot $`K=\mathrm{}`$  and also for some non-trivial knots . It was found that for the models of SAP and random polygon, the knotting probability as a function of the step number $`N`$ is given by the following:

$$P_K(N)=C(K)N^{m(K)}\mathrm{exp}[N/N(K)],$$
(1.1)
where $`C(K)`$, $`m(K)`$ and $`N(K)`$ are fitting parameters. For large $`N`$, the formula (1.1) is consistent with the asymptotic expansion of the partition function of the SAP with the fixed knot type. For finite $`N`$, although it is not rigorous, it seems that the formula (1.1) fits numerical data well. From the numerical result it was conjectured that the parameter $`N(K)`$, which we call the characteristic length of knot $`K`$, should be given by the same value for any knot $`K`$ . Furthermore, it was also conjectured that the parameter $`m(K)`$, which we call the exponent of knot $`K`$, should be universal for different models of SAP or random polygon . We note that the fitting formula (1.1) together with the two conjectures are consistent with the standard asymptotic behavior expected for SAP or random polygon. Here we also note that the rod-bead model is an off-lattice model of SAP.

Recently, the SAP on the cubic lattice with fixed knot was studied through a numerical simulation using the BFACF algorithm, which generates SAPs with the same fixed knot but different step numbers $`N`$ . In the simulation, the exponent $`m(K)`$ and the growth constant for the number of allowed configurations of the SAP with knot $`K`$ has been estimated for some knots. Furthermore, it was shown that the knotting probability for the trivial knot decays “exponentially” for a face-centered cubic lattice  and for the cubic lattice . However, any precise estimate of the knotting probability or the characteristic length $`N(K)`$ has not been given for the cubic lattice. Thus, it is the primary purpose of this paper to evaluate the characteristic length $`N(K)`$ for the SAP on the cubic lattice.

In 1962, Delbrück  noticed that the topological constraint may be very important for polymers in biology and chemistry. Since then, the topological problem has been studied in several approaches in physics and mathematics. As one of the earliest studies, des Cloizeaux and Mehta  estimated through numerical simulations the critical exponent $`\nu `$ for the internal distance of the Gaussian random polygon and discussed some possible properties of random polygons under the topological constraint.

After the rediscovery of the pivot algorithm, many properties of SAW and SAP have been investigated not only in the field theory but also by computer simulations . It seems, however, that there are only a few works such as  where the gyration radius $`R_G^2_{SAP}`$ of the SAP on the cubic lattice is studied by numerical simulations. Let us denote by $`R_G^2_{SAW}`$ and $`R_G^2_{SAP}`$ the mean square of the gyration radius of SAW and SAP, respectively. It is interesting to evaluate the universal amplitude ratio between $`R_G^2_{SAP}`$ and $`R_G^2_{SAW}`$. The ratio has been evaluated only up to $`O(\epsilon )`$ in the $`\epsilon `$-expansion method . Furthermore, according to the scaling theory of polymers, the exponent $`\nu _{\mathrm{𝖲𝖠𝖯}}`$ for $`R_G^2_{SAP}`$ should be given by the exponent $`\nu _{\mathrm{𝖲𝖠𝖶}}`$ for $`R_G^2_{SAW}`$. The agreement is confirmed up to $`O(\epsilon )`$ by the renormalization group theory . However, it is nontrivial to confirm the agreement for the SAP on the cubic lattice through numerical simulations. Thus, the numerical study of the gyration radius for the SAP on the cubic lattice is another purpose of this research.

Hereafter in this introduction, we explain some of the main results of our numerical simulations. Employing the pivot algorithm, we construct a large number of SAPs for the SAP on the cubic lattice with given step number $`N`$. For the gyration radius, we have obtained the exponent: $`\nu _{\mathrm{𝖲𝖠𝖯}}=0.5867\pm 0.0017`$. This is indeed in good agreement with the estimate of the critical exponent of the SAW in the $`\epsilon `$-expansion: $`\nu _{\mathrm{𝖲𝖠𝖶}}=0.5882\pm 0.0011`$. Thus, our simulation in the paper also confirms the agreement of the two exponents.

Let us explicitly consider the method for the topological constraint on the SAP in the cubic lattice with given step number $`N`$. The pivot algorithm of the SAP can generate all the allowed configurations of the SAP with equal probability. Therefore, the set of SAPs generated by the algorithm may contain various knots. Suppose that we have constructed $`M`$ SAPs of step number $`N`$. Calculating some knot invariants for each of the SAPs, we effectively detect the knot types of the SAPs. We enumerate the number of such SAPs that have the same set of values of the knot invariants for knot $`K`$, and denote it by $`M(K)`$. The expectation value of a physical quantity under the topological constraint with fixed knot $`K`$ can be effectively calculated by taking the statistical average of the quantity for the $`M(K)`$ SAPs.

We now turn to the knotting probability. Let the symbol $`P_K(N)`$ denote the knotting probability of the SAP with $`N`$ steps on the cubic lattice for knot $`K`$. If there are $`M(K)`$ SAPs in the total $`M`$ SAPs, then we evaluate it by $`P_K(N)`$$`=M(K)/M`$. In our simulation, $`10^5`$ SAPs ($`M=10^5`$) are constructed for six different values of the step number from $`N=500`$ to $`N=3000`$. We have found that almost all SAPs are topologically equivalent to the trivial knot, and also that the resulting values of the knotting probability for the trivial knot are fitted well by the formula (1.1). Thus, for the trivial knot, we have obtained an estimate of the characteristic length

$$N(\mathrm{})=(2.5\pm 0.3)\times 10^5.$$
(1.2)
This result means that trivial knots are dominant among SAPs on the cubic lattice when the step number $`N`$ is less than $`10^5`$. It implies that when $`N>N(\mathrm{})`$, the majority of SAPs on the cubic lattice have some non-trivial knots. The large value of the characteristic length might be a consequence of strong self-avoiding effect of the SAP on the cubic lattice.

This paper is organized as follows. In section 2 we will estimate the universal amplitude ratio and the critical exponent $`\nu _{\mathrm{𝖲𝖠𝖯}}`$ of the gyration radius. In section 3 we will explicitly discuss the knotting probability for the trivial knot, and obtain the estimate of the characteristic length. We will also show that the knotting probability for the trivial knot decays almost linearly since the characteristic length of the trivial knot is so large. Finally, we will discuss a possibility that the characteristic lengths take an unique value without depending on knot types, in the last section.

## 2 The mean-square of the gyration radius

### 2.1 Previous results on SAW and SAP

An $`N`$-step Self-Avoiding Walk (SAW) $`w`$ in $`^3`$ is a sequence $`w_0,w_1,\mathrm{},w_N`$ of $`N+1`$ distinct points in $`^3`$ such that each point $`w_i`$ is one of the nearest neighbors of its predecessor $`w_{i1}`$: $`|w_iw_{i1}|=1\text{for}i=1,\mathrm{},N`$. It is also subject to a constraint that any site never be occupied by two or more points in a sequence $`\{w_i\}`$. The points $`w_0`$ and $`w_N`$ are the endpoints of $`w`$. The components of $`w_i`$ are represented by $`w_i^{(\alpha )}\text{for}\alpha =1,2,3`$. The Self-Avoiding Polygon (SAP) is a special case of the SAW that makes a ring. We consider an ($`N1`$)-step SAW and denote it by $`w=\{w_0,w_1,\mathrm{},w_{N1}\}`$. If the endpoint $`w_0`$ is the nearest neighbor of the endpoint $`w_{N1}`$, this is an $`N`$-step SAP in $`^3`$.

The mean-square of the gyration radius of the SAP $`R_G^2_{SAP}`$ is smaller than that of the SAW $`R_G^2_{SAW}`$. This fact comes from the following difference. The endpoints of SAW are free, while those of SAP are constrained, that is, the endpoints of SAP should meet at the nearest neighbor sites of the lattice. This means that the SAW has a possibility of a longest end-to-end distance whereas the SAP does not.

The RG argument gives several results of the SAW and the SAP. Among these, we are most interested in the amplitude ratio $`R_G^2_{SAP}/R_G^2_{SAW}`$. According to the RG theory, the amplitude ratio should be universal, i.e. it does not depend on the details of the models for the SAW or the SAP. The ratio has been evaluated by using the RG equation and the cluster expansion up to $`O(\epsilon )`$ :

$$\frac{R_G^2_{SAP}}{R_G^2_{SAW}}=0.568.$$
(2.1)
The ratio has also been estimated in Ref.  through numerical simulations with $`N800`$: It is given by $`0.538\pm 0.006`$.

Another interesting result of the RG argument is that the critical exponent for the mean-square of the gyration radius is universal. It is well-known that the exponent $`\nu _{\mathrm{𝖲𝖠𝖶}}`$ of SAW corresponds to the critical exponent of the O$`(n)`$ vector model in the limit of $`n`$ going to zero. The precise estimate of the critical exponent $`\nu _{\mathrm{𝖲𝖠𝖶}}`$ has been made by the $`\epsilon `$-expansion method through this correspondence :

$$\nu _{\mathrm{𝖲𝖠𝖶}}=0.5882\pm 0.0011.$$
(2.2)
The exponent $`\nu _{\mathrm{𝖲𝖠𝖶}}`$ has also been precisely evaluated by the Monte Carlo simulations, where the best estimate is given by the following :

$$\nu _{\mathrm{𝖲𝖠𝖶}}=0.5877\pm 0.0006.$$
(2.3)
### 2.2 The amplitude ratio

In our simulation of SAPs on the cubic lattice, we have employed a length-conserving dynamical algorithm which keeps endpoints fixed. This algorithm was investigated in detail by Madras, Orlitsky, and Shepp . We call this algorithm MOS pivot for short. Some details are described in Appendix A.

We calculate the ratio $`R_G^2_{SAP}/R_G^2_{SAW}`$ and the exponent $`\nu _{\mathrm{𝖲𝖠𝖯}}`$ using the MOS pivot and compare resulting values with theoretical and other simulated values. The mean-square of the gyration radius is calculated by

$$R_G^2_{SAP}=\frac{1}{T}\underset{t=1}{\overset{T}{}}\left[\frac{1}{N}\underset{i=0}{\overset{N}{}}\left(w_i^{[t]}\frac{1}{N}\underset{j=0}{\overset{N}{}}w_j^{[t]}\right)^2\right],$$
(2.4)
where $`T`$ is the number of polygons and $`N`$ the step number, and the symbol $`w_i^{[t]}`$ denotes the $`i`$-th site in the $`t`$-th SAP in the sequence of SAPs given by “200 successful MOS pivot operations”. (See Appendix A for details.)

Let us note that the gyration radius should have the asymptotic behavior

$$R_G^2_{SAP}=A_{\mathrm{𝖲𝖠𝖯}}N^{2\nu _{\mathrm{𝖲𝖠𝖯}}}(1+BN^\mathrm{\Delta })\text{as}N\mathrm{}.$$
(2.5)
Choosing $`\frac{1}{2}`$ for the exponent $`\mathrm{\Delta }`$ of the correction term, the numerical data given by the equation (2.4) are plotted in figure 1 and fitted to $`A_{\mathrm{𝖲𝖠𝖯}}N^{2\nu _{\mathrm{𝖲𝖠𝖯}}}(1+BN^{\frac{1}{2}})`$, where $`A_{\mathrm{𝖲𝖠𝖯}}`$, $`B`$ and $`\nu _{\mathrm{𝖲𝖠𝖯}}`$ are fitting parameters. Since there are not enough data to fit in a four-parameter curve in our simulation, we have assumed the fixed value for the exponent $`\mathrm{\Delta }`$. We plot the mean-square of the gyration radius from $`N=`$500 to $`N=`$30000 in figure 1. Here, we set $`T`$, the number of polygons, to 10000. The error bars denote one standard deviation given by the Poisson distribution to the number $`T`$ of polygons with $`N`$.

The three parameters are obtained by fitting the data in figure 1

$`\nu _{\mathrm{𝖲𝖠𝖯}}`$ $`=`$ $`0.5867\pm 0.0017,`$ (2.6)
$`A_{\mathrm{𝖲𝖠𝖯}}`$ $`=`$ $`0.1101\pm 0.0037,`$ (2.7)
$`B`$ $`=`$ $`0.06\pm 0.04.`$ (2.8)
The errors are subjective $`68.3\%`$ confidence intervals. For the SAW, we recall that Madras et al.  estimated the exponent $`\nu _{\mathrm{𝖲𝖠𝖶}}`$ by the pivot algorithm up to $`N=`$80000. The estimate of $`\nu _{\mathrm{𝖲𝖠𝖶}}`$ is given by (2.3) together with the following estimates of the fitting parameters

$`A_{\mathrm{𝖲𝖠𝖶}}`$ $`=`$ $`0.19455\pm 0.00007,`$ (2.9)
$`B^{}`$ $`=`$ $`0.11432\pm 0.00465,`$ (2.10)
$`\mathrm{\Delta }`$ $`=`$ $`0.56\pm 0.03.`$ (2.11)
Here, $`B^{}`$ corresponds to $`AB`$ in our estimated values. The estimate (2.6) is in good agreement with that of the RG argument (2.2) and that of the simulation for the SAW (2.3).

From (2.7) and (2.9), we obtain the amplitude ratio

$$\frac{R_G^2_{SAP}}{R_G^2_{SAW}}=0.566\pm 0.019,$$
(2.12)
where an error is a subjective $`68.3\%`$ confidence interval. The amplitude ratio (2.12) is consistent with the estimated value of the RG argument (2.1). The estimate of Ref. is a littile smaller than (2.12). However, it should be consistent with the estimate of Ref. , if we consider the fact that the step numbers $`N`$ of the simulations in Ref.  are much less than 30000.

## 3 The characteristic length $`N(\mathrm{})`$

### 3.1 A method of evaluation of the knotting probability

We evaluate the knotting probability in the following way. We generate $`M`$ SAPs of $`N`$ steps and then enumerate the number $`M(K)`$ of those polygons which are equivalent to a given knot type $`K`$. We define the knotting probability $`P_K(N)`$ by the ratio $`M(K)/M`$.

In our simulation, we determine knot types of SAPs using the 2nd-order Vassiliev-type invariant and the Alexander polynomial evaluated at $`t=1`$. The Vassiliev-type invariants have the following advantages: (1) we can calculate them in polynomial time, and (2) we can calculate them without consuming a large memory area . The Vassiliev-type invariants are not complete invariants. However, in the practical sense we can safely say that if the value of the Vassiliev-type invariant computed for an SAP is zero, this SAP is a trivial knot. We will see in section 3.3 that complicated knots are very rare events in our data. Even if non-trivial knots are misidentified as the trivial knot by the Vassiliev-type invariant (this chance is very small), it would not affect the results of this paper.

For calculating the knotting probability, we generate random sequences of SAPs. We make five seeds for each of the six step numbers $`N=`$500, 1000, 1500, 2000, 2500, and 3000, and then apply an operator $`P`$ on them for a large number of times (20000 times). Here, $`P`$ denotes 200 successful MOS pivots. (See Appendix A.) Then, we get effectively random sequences of SAPs. A sequence of SAPs derived from each seed has a set of 20000 SAPs. Here the sequence consists of 20000 $`\times `$ 200 successful MOS pivot operations. For each step number $`N`$, we thus get samples of 100000 SAPs.

In order to analyze the behavior of the knotting probability, we use the fitting formula (1.1). Here, we write it again:

$$P_K(N)=C(K)N^{m(K)}\mathrm{exp}\left[\frac{N}{N(K)}\right],$$

where $`C(N)`$, $`N(K)`$, and $`m(K)`$ are fitting parameters. In particular $`N(K)`$ is called characteristic length with knot type $`K`$. This formula was introduced by Deguchi and Tsurusaki . They pointed out that the formula (1.1) is suitable for the knotting probabilities of the Gaussian and rod-bead models. We will show that this is suitable also for the knotting probability $`P_K(N)`$ of the cubic lattice model.

### 3.2 The random events of knotted polygons

Let us discuss statistics of knotted polygons generated in our simulation. We will see in section 3.3 that almost all polygons are of trivial knots. Therefore, we may assume that non-trivial knots are generated such as the Poisson random events: the number of trivial knots between two knotted SAPs will follow the Poisson distribution, if the SAPs are randomly constructed.

Let us consider the “time interval” of the Poisson random events. We recall that the time $`t`$ is a discrete number. We measure the length $`L`$ of the time interval of SAPs from time $`t_1`$ at which a knotted SAP appears to time $`t_2`$ at which the next knotted SAP appears after $`t_1`$, and set $`L=t_2t_1`$. If non-trivial knots are generated as the Poisson random events, the time interval $`L`$ follows the function:

$`D(N,L)=A(N)P_{\mathrm{}}(N)^{L1}.`$ (3.1)
We call it the discrete distribution function of the time interval $`L`$. Here, $`P_{\mathrm{}}(N)`$ is the knotting probability for the trivial knot. We have introduced $`A(N)`$ for a technical reason and hence $`D(N,L)`$ is not necessarily normalized.

The discrete distribution functions of the time interval $`L`$ are numerically evaluated as follows. The time interval $`L`$ is a discrete random variable. We introduce a sequence of natural numbers $`\{l_0,l_1,\mathrm{}\}`$, where $`l_i=50\times i`$. Then, we count the number of polygons with $`L`$ taking values in $`[l_i,l_{i+1})`$ and plot it at $`l_i`$ for each $`i`$. This is the discrete distribution function of the time interval $`L`$ obtained from numerical evaluation.

We plot the discrete distribution functions obtained by the numerical evaluation and $`D(N,L)`$ as a function of $`L`$, where $`P_{\mathrm{}}(N)`$ is estimated from table 1 and $`A(N)`$ is chosen to fit these distributions (See figure 2a and 2b). Error bars denote one standard deviation. They are estimated by applying the Poisson distribution to the number of samples included in interval $`L`$. For $`N=`$1000, 1500, 2000, 2500, and 3000, these graphs show a fairly good agreement with the function $`D(N,L)`$. In the case of N=500, the data deviate from the function. This is not unexpected since we have too few samples of knotted polygons. We do not plot the graph for $`N=`$500 in this paper.

We expect from figure 2a and 2b that the MOS pivot makes a uniform ensemble for knots. Thus, we can calculate the knotting probability and the characteristic length using the MOS pivot.

### 3.3 The number of unknotted polygons

The table 1 gives the number of each knot type with respect to $`N`$. Here, errors correspond to 68.3% confidence intervals. They are estimated by applying the binomial distribution to the number $`M(K)`$ of polygons for knot $`K`$. In the table 1 we explain the notations: $`\mathrm{}`$,$`3_1`$,$`4_1`$ denote the trivial knot, the trefoil knot and the figure eight knot, respectively. The other knot types are denoted by $`etc`$.

The table 1 tells clearly that almost all the generated SAPs are trivial knots. $`M(3_1)`$ is also much larger than $`M(4_1)`$. The other knots ($`etc.`$) are nearly equal to zero.

Next we focus on the knotting probability for the trivial knot, and plot $`P_{\mathrm{}}(N)`$ as a function of $`N`$ (figure 3). The error bars are one standard deviation. $`P_{\mathrm{}}(N)`$ decays linearly with respect to $`N`$. It is expected that $`P_{\mathrm{}}(N)`$ becomes decays exponentially when $`N`$ goes to infinity. In  it was shown that the knotting probability $`P_{\mathrm{}}(N)`$ tends to zero “exponentially”:

$$P_{\mathrm{}}(N)=A\mathrm{exp}[\kappa N+o(N)],\text{when}N\mathrm{}.$$
(3.2)
Thus, this is a natural situation.

The asymptotic shape (3.2) is realized for the trivial knot in our case when the fitting parameters of the formula (1.1) take the following values: $`|m(\mathrm{})|1`$ and $`C(\mathrm{})1`$. In fact, using the least-squares estimation, we find:

$`C(\mathrm{})`$ $`=`$ $`1.0035\pm 0.0035`$ (3.3)
$`m(\mathrm{})`$ $`=`$ $`(4.7\pm 5.7)\times 10^4`$ (3.4)
$`N(\mathrm{})`$ $`=`$ $`(2.5\pm 0.3)\times 10^5`$ (3.5)
$`\chi ^2`$ $`=`$ $`0.748`$ (3.6)
$`Prob(\chi ^2>0.748)`$ $`=`$ $`0.862`$ (3.7)
(errors are one standard deviation). Here, the $`\chi ^2`$ value (i.e. sum of square of normalized deviations from the regression line) can serve as a criterion of good fit. It should be distributed as $`\chi ^2`$ with $`n3`$ degrees of freedom, where $`n`$ is the number of data points in the fit. $`Prob(\chi ^2>0.748)`$ is the probability that $`\chi ^2`$ would exceed the observed value, in this case $`86.2\%`$. This implies that the formula (1.1) is suitable. It is remarkable that the characteristic length of the trivial knot is much larger than the value expected from the rod-bead model.

There have been a few simulation studies on several lattices. In  the knotting probability for the trivial knot was calculated on a face-centered cubic (FCC) lattice. It was shown that assuming the two-parameter fitting formula: $`P_{\mathrm{}}(N)=C(\mathrm{})e^{\alpha (\mathrm{})N}`$, the parameters were given by

$`\alpha (\mathrm{})`$ $`=`$ $`(7.6\pm 0.9)\times 10^6`$
$`C(\mathrm{})`$ $`=`$ $`1.0011\pm 0.003`$
$`\chi ^2`$ $`=`$ $`2.7`$
$`Prob(\chi ^2>2.7)`$ $`=`$ $`0.44,`$
where errors were one standard deviation. Our interpretation of $`\alpha (\mathrm{})`$ taking $`7.6\times 10^6`$ is that the characteristic length is $`1.3\times 10^5`$. Comparing the cubic lattice with the FCC lattice, the characteristic length is larger on the cubic lattice than on the FCC lattice. In  it was also shown that the exponent $`\alpha (\mathrm{})`$ in the above form was $`(5.7\pm 0.5)\times 10^6`$ on the cubic lattice. This corresponds to the characteristic length of $`1.8\times 10^5`$. We can expect that this is consistent with our estimated value within the error bars. Here we note that the connection of Ref.  shall be discussed in section 4 and also note that the simulation results of Refs.  seem to contain some information on the characteristic length of the SAP on the cubic lattice. However, we are unable to derive any appropriate estimate from them.

Let us return to our data. We estimate not only the characteristic length $`N(\mathrm{})2.5\times 10^5`$ but also the exponent of a correction term $`m(\mathrm{})0`$. These parameters show that the knotting probability for the trivial knot serves a “pure exponential” decay on the cubic lattice. This is a new result for the cubic lattice model. Thus, the simulation in the paper improves that of .

We have several interpretations of the characteristic length $`N(\mathrm{})`$. Since $`N(\mathrm{})`$ is so large, we found that the number of knotted polygons are much smaller than that of unknotted polygons in our simulation. We expect that the knotting probability for the trivial knot decreases to about $`30\%`$ at $`N=N(\mathrm{})`$. Thus, non-trivial knots become majority of SAPs when $`N`$ is larger than $`N(\mathrm{})`$. In addition to this, we will see in section 4 that $`N(\mathrm{})N(K)`$ for any knot type $`K`$ is due to the fact that $`N(\mathrm{})`$ is large.

In figure 4, we plot the knotting probability for the trefoil knot $`P_{3_1}(N)`$ as a function of $`N`$, where error bars are one standard deviation. The data points almost lie on a straight line. Since the data of the trefoil knot are only six points in the step number, we cannot fit the knotting probability $`P_{3_1}(N)`$ to formula (1.1). We leave the selection of fitting parameters for the trefoil knot or more complicated knots as future investigations.

The knotting probability $`P_{3_1}(N)`$ may have a finite-size effect. When we set $`m(3_1)=1`$ as expected from , the formula (1.1) does not match with $`P_{3_1}(N)`$ in our data. Deguchi and Tsurusaki argued that a finite-size effect appeared in the knotting probabilities in . When we plot the straight line fitted to $`P_{3_1}(N)`$, the line intersects with the $`x`$-axis at a positive value. This is the finite-size effect. We introduce offset parameter $`N_{ini}`$ and replace $`N`$ by $`\stackrel{~}{N}=NN_{ini}`$ in the formula (1.1). When we fix $`m(3_1)=1`$, our rough estimation gives as $`N_{ini}140`$ for the trefoil knot. On the other hand, we roughly estimate $`N_{ini}0`$ for the trivial knot. We expect that such finite-size effect would also appear in the knotting probabilities for the SAPs with more complicated knot types. It could be confirmed by generating SAPs with much larger $`N`$.

In section 2.2 we calculated the gyration radius including all possible knots, and then estimated the universal amplitude ratio and the universal exponent $`\nu `$. Without classifying knot types, however, we can effectively consider only trivial knots when $`N<N(\mathrm{})`$. The universal amplitude ratio and the universal exponent are evaluated effectively for the trivial knots, although 2% of the SAPs are non-trivial knots and we neglect their influence.

## 4 A consequence of <br>the large characteristic length $`N(\mathrm{})`$

The trivial knot dominates among SAPs on the cubic lattice when $`N`$ is less than $`2.5\times 10^5`$. Our interpretation is that the cubic lattice has so strongly the excluded-volume effect that it almost prevents appearances of knotted polygons.

We have a conjecture that the appearance of a complicated knot is a rare event on the cubic lattice. According to , non-trivial knots occupy a large number of configurations of SAPs for $`NN(\mathrm{})`$. We believe that the above situation is realized also on the cubic lattice. If the formula (1.1) is a suitable form of the knotting probability and if the trefoil and figure eight knots on the cubic lattice behave like these of continuum models, the ratio $`M(4_1)/M(3_1)`$ for each step number $`N`$ should depend only on the ratio $`C(4_1)/C(3_1)`$. We recall that $`M(3_1)M(4_1)`$ for $`N<`$ 3000 from the table 1. Then, at $`N2.5\times 10^5`$ almost all SAPs are expected to be $`3_1`$ knotted polygons on the cubic lattice unlike SAPs on continuum models.

We see that our estimated value $`N(\mathrm{})`$ is related to the growth constant $`\mu `$ from the viewpoint of Orlandini et al.. The asymptotic behavior of the number of $`N`$-step polygons $`c_N`$ is given by

$$c_N=aN^{\alpha 3}\mu ^N\left(1+bN^\mathrm{\Delta }+o(N^1)\right),$$
(4.1)
where $`\alpha `$ and $`\mathrm{\Delta }`$ are critical exponents. For fixed knot type $`K`$, it is believed that the number $`c_N(K)`$ of polygons with knot $`K`$ should have a similar form:

$$c_N(K)=a(K)N^{\alpha (K)3}\mu (K)^N\left(1+b(K)N^{\mathrm{\Delta }(K)}+o(N^1)\right).$$
(4.2)
Then, the knotting probability $`P_K(N)`$ is given by $`c_N(K)/c_N`$, and this implies that the characteristic length $`N(\mathrm{})`$ relates to the growth constants $`\mu `$ and $`\mu (\mathrm{})`$. We can estimate the ratio $`\mu /\mu (\mathrm{})`$ from the value of $`N(\mathrm{})`$:

$$\frac{\mu }{\mu (\mathrm{})}e^{1/N(\mathrm{})}1+(4\pm 2)\times 10^6,$$
(4.3)
where an error is two standard deviations.

Let us discuss the independence of the characteristic length $`N(K)`$ with respect to knot type $`K`$. Orlandini et al.  calculated the growth constants directly and showed the following equality:

$$\mu (\mathrm{})=\mu (K)=4.6836\pm 0.0038\text{for}anyknottypeK,$$
(4.4)
where an error corresponds to $`95\%`$ confidence interval. In addition Guttmann estimated the growth constant

$$\mu =4.68393\pm 0.00002$$
(4.5)
using exact enumeration and series analysis (an error is one standard deviation). From (4.4) and (4.5), we obtain the ratio

$$\frac{\mu }{\mu (\mathrm{})}=\frac{\mu }{\mu (K)}1+(7\pm 8)\times 10^5,$$
(4.6)
where an error is two standard deviations. From (4.3) and (4.6), we expect that the characteristic length is independent of the knot type: $`N(\mathrm{})N(K)2.5\times 10^5`$ for any knot type $`K`$.

Although the independence of $`\mu (K)`$ with respect to knot type $`K`$ is pointed out by Orlandini et al. , we can also confirm it more precisely through our simulation. While Orlandini et al.  calculated the growth constants, we have no direct calculation of them. However, we can predict that the differnce between $`\mu (\mathrm{})`$ and $`\mu (K)`$ is very small when we use the following inequalities: $`\underset{N\mathrm{}}{lim\; inf}N^1\mathrm{log}c_N(K)\mu (\mathrm{})`$ and $`\underset{N\mathrm{}}{lim\; sup}N^1\mathrm{log}c_N(K)<\mu `$ for any knot type $`K`$, which were proven in . If both $`\underset{N\mathrm{}}{lim\; inf}N^1\mathrm{log}c_N(K)`$ and $`\underset{N\mathrm{}}{lim\; sup}N^1\mathrm{log}c_N(K)`$ exist and take the same value $`\mu (K)`$, then we have $`\mu (\mathrm{})\mu (K)<\mu `$. These inequalities and the estimate (4.3) limit the ratio $`\mu (K)/\mu (\mathrm{})`$ to

$$\left|\frac{\mu (K)}{\mu (\mathrm{})}1\right|(4\pm 2)\times 10^6.$$
(4.7)
This is a strong bound. Thus, we can expect that $`\mu (\mathrm{})\mu (K)`$ for any knot type $`K`$.

## Acknowledgement

We would like to thank Takeo Inami for a careful reading of the manuscript and valuable comments. A.Y. is supported by a Research Assistant Fellowship of Chuo University. This work is partially supported by the Grant-in-Aid for Encouragement of Young Scientists (No. 12740231).

Appendix

## Appendix A The method for generating SAPs

Let us discuss the method for making SAPs in the paper. We first construct a seed SAP $`w`$, and then derive a random sequence of SAPs $`\{w^{[1]},w^{[2]},\mathrm{},w^{[t1]},w^{[t]},\mathrm{}\}`$, where $`w^{[t]}=P(w^{[t1]})`$ and $`w^{[0]}=w(0)`$. Here the operator $`P`$ will be defined later.

We construct a ‘seed’ SAP, by combining two SAWs which have the same endpoints, in the following way: first we make an $`N/2`$-step SAW using the myopic self-avoiding walk (MSAW) algorithm , where $`N`$ is an even integer; secondly we perform the MOS pivot transformations with respect to $`k=0,l=N/2`$. ; finally we concatenate the endpoints of the new and original SAWs respectively (figure 5), and get an SAP with the step number $`N`$ if it has no self-intersections.

Let us discuss how to construct a random sequence of SAPs in our simulations. If one of the MOS pivot transformations changes a given SAP into a different SAP, this operation is called a successful MOS pivot operation. We consider a sequence of successful MOS pivot operations. Then, we define an operator $`P`$ by 200 successful MOS pivot operations in the sequence. We note that the SAP obtained by a single successful MOS pivot is not independent from the original one: they are correlated. However, the correlation decays almost completely after 200 successful MOS pivot operations, which shall be shown in Appendix B. Thus, we may consider that for any given SAP $`w`$, $`P(w)`$ is independent from $`w`$.

For generating SAPs, we use the Mersenne Twister that is a pseudo random number generator . This algorithm has the following properties: (1) we can get many samples because the period is $`2^{19937}1`$, (2) we treat high dimensional space (max $`623`$ dimensions), (3) pseudo random numbers are generated fast and (4) we can use the memory efficiently. Thus, the Mersenne Twister is a high performance generator.

## Appendix B The decay of correlations among the SAPs obtained by the MOS pivot operations

In order to check the validity of the random sequence of SAPs constructed in the paper, we show explicitly how the correlation among the SAPs decays after applying a number of MOS pivot operations.

Let us regard the number $`\tau `$ of successful MOS pivot operations as the time of evolution of the shape of SAP under the sequential MOS pivot operations: the seed SAP $`w`$ of a random sequence corresponds to $`w(0)`$ for $`\tau =0`$; $`w(\tau )`$ is defined by $`p^\tau (w)`$ for $`\tau >0`$. Here $`p`$ have denoted a successful MOS pivot operation. Let us now define the correlation function for the structure of SAP with the step number $`N`$ in the following

$$C(\tau )=\frac{_{\alpha =1}^3_{i=0}^{N1}\left(w_i^{(\alpha )}(0)\frac{1}{N}_{j=0}^{N1}w_j^{(\alpha )}(0)\right)\left(w_i^{(\alpha )}(\tau )\frac{1}{N}_{l=0}^{N1}w_l^{(\alpha )}(\tau )\right)}{_{\alpha =1}^3_{i=0}^{N1}\left(w_i^{(\alpha )}(0)\frac{1}{N}_{j=0}^{N1}w_j^{(\alpha )}(0)\right)^2},$$
(2.8)
where $`w_i^{(\alpha )}(\tau )`$ is the $`\alpha `$-component of the $`i`$-th site of SAP (or SAW) after $`\tau `$ pivot operations and $``$ denotes the statistical average.

In figure 6a and 6b, we plot the correlation $`C(\tau )`$ versus the time $`\tau `$ in the case of $`N=`$ 500 and 1000, respectively. Error bars show one standard deviation which are estimated by assuming that the data follows the Poisson distribution. In our simulation, we evaluated the correlations of SAPs with the step number $`N`$ at 500 and 1000, and generated 10000 sequences (starting from 10000 different seeds) to take the statistical average.

The decay rate of the correlation function $`C(\tau )`$ is slower in the case of the MOS pivot than in the case of normal pivot algorithm (figure 6a and 6b). This is due to the difference in the numbers of independent pivot transformations. Note that the number of independent transformations for the normal pivot algorithm  is $`d!2^d1`$ (the normal pivot algorithm is an algorithm with one endpoint free while the other fixed) , which correspond to the number of all the elements of $`d`$-dimensional orthogonal group on a hypercubic lattice, while for the MOS pivot $`2d(d1)+1`$ in the $`d`$-dimensional hypercubic lattice. In the cubic lattice, the MOS pivot has 13 pivot transformations while the normal pivot algorithm has 47.

## Figure Caption