A study of polymer knots using a simple knot invariant written consisting of multiple contour integrals

# A study of polymer knots using a simple knot invariant written consisting of multiple contour integrals

## Abstract

In this work the thermodynamic properties of short polymer knots (up to 120 segments) defined on a simple cubic lattice are studied with the help of the Wang-Landau Monte Carlo algorithm. The sampling process is performed using pivot transformations starting from a given seed conformation. Both cases of short-range attractive and repulsive interactions acting on the monomers are considered. The properties of the specific energy, heat capacity and gyration radius of the knots and are discussed. It is found that the heat capacity exhibits a sharp peak. If the interactions are attractive, similar peaks have been observed also in single open chains and have been related to the transition from a frozen crystallite state to an expanded coil state. Some other peculiarities of the behavior of the analyzed observables are presented, like for instance the increasing or decreasing of the knot specific energy at high temperatures with increasing polymer lengths depending if the interactions are attractive or repulsive. Besides the investigation of the thermodynamics of polymer knots, the second goal of this paper is to introduce a method for distinguishing the topology of a knot based on a topological invariant which is in the form of multiple contour integrals and explicitly depends on the physical trajectory of the knot. The chosen invariant, denoted here , is related to the second coefficient of the Conway polynomial. It has been first isolated from the amplitudes of a Chern-Simons field theory with gauge group . It is shown that this invariant is very reliable in distinguishing the topology of polymer knots. One of the advantages of the proposed approach is that it allows to reduce the number of samples needed by the Wang-Landau algorithm. Some solutions to speed up the calculations of exploiting Monte Carlo integration techniques are developed.

## I Introduction

Long polymers are very likely to be found in the configuration of knots or links. The topological properties of polymers with closed conformations play indeed an important role in physics, chemistry and biology. For that reason, they are being actively investigated grosberg (); dna (); katritch96 (); katritch97 (); krasnow (); laurie (); cieplak (); marko (); liu (); wasserman (); sumners (); vologodski (); orlandini (); marenduzzoorlandinimichelettiphysrep (); levene (); kurt (); mehran (); pp2 (); yan (); arsuaga (); arsuaga2 (); metzler (); pieranski (); diao (); sumners2 (); faddeev (). A particularly challenging problem is that of the statistical mechanics of polymer knots. Up to now, a satisfactory analytical model exists only in the case of two polymer rings linked together FFIL (); FFIL2 (), but there is no analogous model for a knot despite many attempts, see for instance kholodenkovilgis (); FFAP (); FFNOVA (); kleinert () for a review on this subject. Moreover, the scaling laws of the most important observables of polymer knots, like for instance the gyration radius, are still a subject of intense research grosbergjpa (); stella (); quake (); mansfielddouglas (); suzukietal ().

The main difficulty behind the treatment of polymer knots, as well as polymer links, is to distinguish the wealth of different topological configurations of such systems. This problem arises in analytical models because it is necessary to impose topological constraints in order to avoid that the statistical fluctuations affect the initial topological state. This is a physical requirement, dictated by the fact that, once a polymer knot or link has been formed, its topological state cannot be modified without breaking the covalent bonds holding together the monomers. Thus, unwanted changes of the topological configuration must be detected and rejected. In numerical simulations, instead, a widely used way to generate polymer knots or links is to consider self avoiding walks (SAW’s) that, at a certain point, intersect themselves forming closed conformations vologodski (); descloizeaux (); hazi (); michels (). The problem in this case is to sort out the topological configurations that are of interest from all the other configurations produced within this approach. This process of generating knot and links is not very efficient for our purposes, because the probability of formation of knots or links of a given type from SAW’s is very low, see for instance Ref. probability1 (); probability2 (); probability3 (). Alternatively, it is possible to start from a polymer ring with a seed trajectory that is out of equilibrium, but already in the desired topological configuration. Later, the system is equilibrated using the so-called pivot transformations pivot (); Rensburg (). There exist pivot transformations that are automatically preserving the topological state of the system, see for example swetnam (), however they are able to modify only a very small part of the whole polymer, a fact that considerably increases the time for reaching the equilibrium, especially for very long polymers. In yzff () another method has been proposed, called the Pivot Algorithm and Excluded Area (PAEA) method, in which more general pivot transformations are considered, including those that could potentially change the topology. In this case, large transformations are not easy to be implemented.

In conclusion, apart from a few exceptions, like for instance the already cited works of swetnam (); yzff () and notably the dynamic Monte Carlo approach dmca (), that however leads to some level of polydispersity, in most numerical computations involving polymer knots, topological invariants are exploited whenever it becomes necessary to distinguish the topological configuration of the system under investigation. The most popular topological invariants are given in the form of polynomials or of multiple contour integrals computed along the physical trajectories of the polymers. The latter invariants are easier to be used in analytical models than the former, because the coefficients of the polynomials are not directly related to the polymer conformation. The particular simplicity of the Gauss linking number, which consists in a double contour integral, is the main reason for which it has been possible to derive an analytical model of two linked polymer rings FFIL (); FFIL2 (); FFIL3 (). Moreover, the Gauss linking number has already been used in numerical simulations of polymer systems, see for instance ralf (); Kremer1 (). Unfortunately, in order to distinguish the topology of knots, there is no simple invariant like the Gauss linking number. So far, the statistical properties of polymer knots have been studied numerically with the help of topological invariants like the Alexander polynomials vologodski () or the HOMFLY polynomials michalke (), see the detailed textbook by Kleinert kleinert () for an extensive review on that subject. Of course, there is a plenty of other knot invariants that have been or could be applied, for instance the Conway polynomials conway (), the Arf-Casson invariant GMMknotinvariant () or the Milnor milnor () and Vassiliev-Kontsevich invariants kontsevich (). The latter three have an explicit representation in terms of multiple integrals computed along the knot trajectories.

The purpose of this paper is to show that knot invariants that are given in the form of multiple contour integrals can represent a valid alternative in numerical calculations to polynomial knot invariants or even to the PAEA method, which is able to detect the topology changes exactly and has been proven to allow very fast computations of the statistical properties of polymer knots yzff (); yf (). In particular, we will concentrate on a topological invariant denoted here , where denotes the trajectory of the knot. has been derived from the one-loop amplitudes of non-abelian Chern-Simons field theories with gauge group GMMknotinvariant (). It is related to the Arf-Casson invariant and to the second coefficient of the Conway polynomials GMMknotinvariant (). Its value can be analytically computed for any given knot configuration. is the simplest knot invariant represented in terms of multiple contour integrals.

The knot invariant is applied here in order to derive the average values of the specific energy, heat capacity and gyration radius of several different knot configurations by means of the Wang-Landau algorithm wanglandau (). We find in this way that the examined knots, corresponding namely to the trefoil , the figure-eight and 1, undergo with the temperature a phase transition, which is probably from a frozen crystallite state to an expanded coil state similarly to what happens in the case of a single polymer chain proposed in binder (). Other physical properties of polymer knots are discussed. Compared with the PAEA method, the use of allows to reduce the number of samples necessary for the calculations of the averages of the observables with the Wang-Landau algorithm. As it will be seen, this reduction is due to the fact that with , large pivot transformations can be exploited which are able to change relevant portions of the knot. In this way, the exploration of the whole set of available conformations becomes faster. Despite the decreasing of the number of samples, the computations last in general longer than those performed with the PAEA method, because the expression of contains quadruple integrals that should be evaluated numerically and this takes some time. As a consequence, we present here the results for relatively short polymer knots up to , where is the number of segments composing the knot on a simple cubic lattice. Actually, there is no problem in studying longer polymer knots. The reason is that invariants given in the form of contour integrals can be computed for arbitrarily deformed knots, not necessarily defined on a lattice, provided the topological configuration remains the same after the deformation. Thanks to this fact, we have found that it is possible to shorten the number of segments by a factor three, speeding up the calculation of considerably. Even using this trick, the method requires too long times (exceeding a month on a modern workstation) if . The use of becomes however competitive in the equilibration of very long polymers because, as already mentioned, it offers the possibility of exploring very fast the set of conformations compatible with the topological constraints.

The rest of the current work is organized as follows. Our simulation set-up and the topological invariant are introduced in Section II. The Monte Carlo integration method we use to calculate the knot invariant is introduced in Section III. An explanation of the Wang-Landau algorithm with particular attention to its applications to the statistical mechanics of polymers is provided in Section IV. The results on the thermal properties of different polymer knots are presented in Section V. We also compare our calculations with the results obtained with the PAEA method in Refs. yzff (); yf (). Finally, in Section VI we draw our conclusions and possible generalizations of this work are briefly discussed.

## Ii Methodology

### ii.1 General outline of the used methodology

First of all, a brief digression on the used terminology is in order. Throughout this work the word configuration refers to a particular topological state of a polymer knot. The word conformation will instead denote the particular shape in the space of the trajectory of a polymer knot in a given topological configuration.

At this point it is possible to go back to the outline of the methodology. We adopt the strategy of considering at the beginning a seed conformation of the knot to be analyzed. The knot is a self-avoiding polygon defined on a simple cubic lattice with edges of unit length. The monomers are located on the vertices of the lattice. Let denote the total length of the polygon. Since it consists of edges or segments of unit length, coincides also with the number of segments composing the polygon. The crossing of the knot trajectory with itself at some point on the lattice is forbidden. The starting seed configuration is equilibrated and then used to compute the thermodynamic properties of the studied knot. The relevant observables are calculated by means of the Wang-Landau Monte Carlo algorithm wanglandau (), which will be explained in some more details later. Both cases of attractive and repulsive interactions are considered. For the equilibration of the knot and the sampling in the Wang-Landau algorithm, the polymer conformations are randomly modified by exploiting the pivot transformations described in Ref. pivot (). These transformations can involve any number of segments such that and are not preventing the crossing of the lines of the knot trajectory . For that reason, the invariant should be applied in order to check if the topology of the knot has been altered after each pivot transformation. If this is the case, the transformation is rejected and a new one is considered.

### ii.2 The topological invariant ϱ(C)

To avoid topology changes of a polymer knot potentially occurring after the pivot transformations, we will use in this work a topological invariant that has been derived from the one-loop amplitude of the Wilson loop in non-abelian Chern-Simons field theories GMMknotinvariant (); LR (); cotta (). The most important characteristic of this invariant, that will be denoted , is that it can be expressed in the form of a sum of multiple contour integrals:

 ϱ(C)=ϱ1(C)+ϱ2(C) (1)

where the knot is represented as an oriented closed path of length . The contribution is given by the triple integral:

 ϱ1(C)=−132π3∮Cdxμ∫xdyν∫ydzρIμνρ(→x,→y,→z), (2)

with

 Iμνρ(→x,→y,→z) =ϵαβγϵμασϵνβλϵργτ (3) ×∫d3→ω(ω−x)σ|→ω−→x|3(ω−y)λ|→ω−→y|3(ω−z)τ|→ω−→z|3

while the second part is:

 ϱ2(C) = 18π2∮Cdxμ∫xdyν∫ydzρ∫zdwσϵσναϵρμβ (4) ×(w−y)α|→w−→y|3(z−x)β|→z−→x|3

In the above formulas greek letters denote space indexes. The variables and , , are the components of the radius vectors describing the positions of four points on the same curve . represents instead the completely antisymmetric tensor uniquely defined by the condition . The integrations along the path in (2) and (4) are path ordered. This can be seen explicitly by parametrizing the trajectory with the arc-length:

 ϱ1(C) = −132π3∫L0dsdxμ(s)ds∫s0dtdyν(t)dt∫t0dudzρ(u)duIμνρ(→x,→y,→z), (5)

and

 ϱ2(C) = 18π2∫L0dsdxμ(s)ds∫s0dtdyν(t)dt∫t0dudzρ(u)du∫u0dvdwσ(v)dvϵσναϵρμβ(w(v)−y(t))α|→w(v)−→y(t)|3(z(u)−x(s))β|→z(u)−→x(s)|3 (6)

It has been shown that the knot invariant appearing above is related to the second coefficient of the Conway polynomial of a knot through the following relation GMMknotinvariant ():

 a2(C)=12[ϱ(C)+112] (7)

The Conway polynomials are well known and their coefficients can be computed analytically for every knot topology, so that thanks to Eq. (7) it is easy to derive also the values of . In Table 1 we give a list of the second coefficients of the Conway polynomial and the corresponding values of for the knot configurations that will be studied here.

is the simplest known knot invariant that can be expressed in the form of contour integrals. Like any other knot invariant, is not able to distinguish different knots unambiguously. For example, the trefoil knot has , exactly the same value of the knots and many others. However, we should keep in mind that the main role of a knot invariant in studying the thermal and mechanical properties of polymer knots is not to guess its topological configuration. In fact, the topological configuration is known since the beginning. The problem is rather to preserve that configuration against thermal fluctuations, because without any constraint the polymer trajectories are allowed to cross themselves, a fact that can potentially alter a knot. The probability that due to thermal fluctuations a polymer ring jumps from one knot configuration to another with the same value of is very low, as it has been observed in our simulations. What emerges from them is that it is very unlikely that a knot passes to another configuration with the same value of after a pivot transformation. Most probably, one ends up with the trivial knot or, somewhat less frequently, with a conformation with lower number of crossings such that . For the purposes of this work, it is thus possible to affirm that is a powerful knot invariant. To convince oneself, it is sufficient to recall that, if one considers the simplest knots up to ten crossings, there are particular topological configurations that are uniquely distinguished by , like and or are very efficiently distinguished from all the others because their corresponding values of occur rarely. Luckily, if the values of for two topologically different knots are not the same, the smallest difference between them is . As an example, the knot has , while for the topologically inequivalent knots and we have and respectively. This allows to choose the number of sampling points in such a way that the variance in the Monte Carlo evaluation of is low enough that the probability of confusing two different knot topologies due to numerical errors is negligible.

The price to be paid for this efficiency in distinguishing knots is the complicated expression of . The most time-consuming contribution to is the quadruple contour integral necessary to compute in Eq. (4). For a knot of length , the evaluation time of scales as . This is approximately one order more than the time necessary to evaluate the Alexander polynomial of a knot vologodski (); deguchitsurusaki (); Muthukumar (), which scales as . Here denotes the number of crossings which is necessary to represent the knot by projecting it on an arbitrary plane, see Muthukumar () for more details. Of course, also the computation of the Alexander polynomial becomes prohibitive for polymers which are long or have compact conformations, because in these cases the number of crossings drastically increases marenduzzoorlandinimichelettiphysrep (). Moreover, the scaling law of the computational time is true only if the determinant of a matrix that arises in the algorithm for computing the Alexander polynomial is evaluated with the method of Gaussian elimination, which is subjected on round-off error that become important when is large.

One advantage of is that its calculation can be extended without any effort to any kind of trajectory, not necessarily on a cubic lattice. This fact is very helpful when the polymer is long, so that it is advisable to decrease the number of its segments. In a very simple way it is possible to reduce by a factor three without destroying the topology by replacing in the knot every group of three contiguous segments with a single segment. As well, there is no problem in distinguishing the changes of topology when the number involved in the pivot transformations becomes large.

Like the Alexander polynomial, also is not able to distinguish uniquely two different topological configurations and is subjected to numerical errors. However, we have seen above that, for the goals of this work, the invariant is powerful enough.

## Iii Monte Carlo evaluation of path ordered contour integrals on a lattice

### iii.1 Simpson’s rule vs. Monte Carlo method

First of all, we introduce some notation that will be useful in this Section. The contour describing the physical trajectory of the knot in space is represented here as a curve , with . Of course, is consisting of a set of discrete segments, see Fig. 1 for an example with , and this fact should be taken into account. To this purpose, let’s denote with , , the locations of the lattice sites through which the closed contour is passing. The -th segment of the loop forms a vector for . The -th segment is instead associated with the vector . Next, let , with , be the restriction of the curve to the th segment. Here we have introduced the segment’s arc-length such that . On the th segment, is related to by the formula

 ~s=s−i for i=1,…,L−1 (8) ~s=s for i=L (9)

Explicitly, the expression of is given by:

 →xi(s) = →xi+~s(→xi+1−→xi)i=1,…,L−1 (10) →xL(s) = →xL+~s(→x1−→xL) (11)

From the above equations, it is easy to derive also the derivative of with respect to . Whenever it will be necessary to specify points on different elements of the loop , for instance on segments , they will be denoted with the symbols , where and .

At this point we are ready to rewrite the quantities and displayed in Eqs. (5) and (6) respectively in a form that is suitable for applying the standard Simpson’s rule:

 ϱ1(C) = −132π3L∑i=1i∑j=1j∑k=1∫10d~sdxμi(~s)d~s∫1−δij(1−~s)0d~tdyνj(~t)d~t (12) × ∫1−δjk(1−~t)0d~udzρk(u)d~uIμνρ(→xi(~s),→yj(~t),→zk(~u))

and

 ϱ2(C) = 18π2L∑i=1i∑j=1j∑k=1k∑l=1∫10d~sdxμi(~s)d~s∫1−δij(1−~s)0d~tdyνj(~t)d~t∫1−δjk(1−~t)0d~udzρk(~u)d~u (13) × ∫1−δkl(1−~u)0d~vdwσl(~v)d~vϵσναϵρμβ(wl(~v)−yj(~t))α|→wl(~v)−→yj(~t)|3(zk(~u)−xi(~s))β|→zk(~u)−→xi(~s)|3

The boundaries in the integrals of and have been chosen in such a way that the path ordering of the trajectories is preserved. For example, if the integrals over and are performed over two different segments, then the boundaries for both variables are ranging between and . Instead, if the integrals are performed on the same segment , then it must be that and . Let us note that the integrals over the variable in Eq. (13) can be performed exactly. The remaining integrals over and should be evaluated numerically, for instance by means of the Simpson’s rule. It turns out that for the evaluation of these integrals, the zeroth order approximation, obtained by substitutions of the kind:

 ∫1−δij(1−~s)0dyνj(~t)fν(→yj(~t)) ∼(yνj(1)−yνj(0))fν(→yj(1))+fν(→yj(0))2 (14)

is sufficient and gives satisfactory results. The problem is that, even exploiting the crude approximations of Eq. (14) in order to distinguish the topology changes of a given polymer knot, still a sum of terms should be evaluated in order to derive the value of from Eq. (13). Already in the case of polymers of length or more, this number becomes prohibitively high for practical purposes. In fact, the Wang-Landau procedure used to compute the density of states requires several millions of samples to be evaluated. For that reason, it is much better to estimate and performing the integration with Monte Carlo techniques.

The idea is to regard the contour integrals in Eqs. (5) and (6) as usual multiple integrals over the variables and :

 ϱ1(C)=∫L0ds∫s0dt∫t0duF1(s,t,u) (15)

and

 ϱ2(C)=∫L0ds∫s0dt∫t0du∫u0dvF2(s,t,u,v) (16)

where

 F1(s,t,u) = −132π3dxμ(s)dsdyν(t)dtdzρ(u)du (17) ×Iμ,ν,ρ(→x(s),→y(t),→z(u))

and

 F2(s,t,u,v)=18π2dxμ(s)dsdyν(t)dtdzρ(u)dudwσ(v)dv ×ϵσναϵρμβ(w(v)−y(t))α|→w(v)−→y(t)|3(z(u)−x(s))β|→z(u)−→x(s)|3 (18)

The variables and in Eq. (15) span a space of volume , while the variables and in Eq. (16) span a space of volume . To evaluate the right hand sides of Eqs. (15) and (16) via Monte Carlo integration, we can exploit the general formula

 ∫b1a1dξ1∫ξ1a2dξ2⋯∫ξm−1amdξmf(ξ1,⋯,ξm)≈1N[N∑i=1f(ξ(i)1,⋯,ξ(i)m)(b1−a1)m∏σ=2(ξ(i)σ−aσ)] (19)

where the ’s, and denote randomly chosen variables in the range:

 [a1,b1]whenσ=1[aσ,ξσ]whenσ=2,…,m (20)

In the numerical evaluation of , we consider the trajectory of the knot, which in principle is a polygon on a simple cubic lattice, as a continuous curve . For a given value of , the segment on which the point is located is identified by the relation

 i = [s]+1 (21)

where denotes the integer part of . The components of the curve are obtained using its restriction to the th segment, whose components can be computed exploiting Eqs. (10) and (11). The components of , and are derived analogously.

We stress the fact that the integrands and are regular, even if for instance seems to be divergent when or . Analytically, it is possible to prove that these singularities cancel as expected in a topological invariant. In numerical computations the situation looks however different, because one has to cope with terms that are separately diverging, but whose sum is finite. A regularization is thus necessary in order to eliminate these ambiguities. Let us first consider the computation of . Here there are potential problems whenever

 →w(v)−→y(t)=0 or →z(u)−→x(s)=0. (22)

However, the probability of the occurrence of such situations by choosing randomly the variable is very low. In fact, if the calculations are performed using double precision variables, the number of digits after the floating point is so high, that in practice the conditions displayed in Eq. (22) are never satisfied. More serious is the case of . After the integration over in Eq. (3), we get the explicit expression of which is reported in Appendix A. We see that, besides the singularities at the points satisfying the conditions:

 →y(t)=→x(s)→z(u)=→x(s)→y(t)=→z(u) (23)

it appears also a pole whenever the equation

 A=|→y(t)−→x(s)||→z(u)−→x(s)|+(→y(t)−→x(s))⋅(→z(u)−→x(s))=0 (24)

is verified. While it is very unlikely that, during the random sampling, the divergences of Eq. (23) will appear due to the same reasons explained in the case of the analogous divergences in Eq. (22), the condition (24) is very easy to be realized. It is sufficient for example that the two vectors and have opposite orientations and two zero components. Depending on the topology of the knot, this situation may occur rather often. Even if the number of cases in which this happens is negligible with respect to the total number of sampled points generated during the Monte Carlo computation of , still the result may be spoiled by overflow or underflow errors. To cure the singularity in Eq. (24) when the quantity is equal to zero, we use the framing regularization introduced in Witten (). This consists in performing an almost infinitesimal shift of the curves , , and along a direction which is normal to the knot . By denoting with the unit vector that gives the normal direction to , this means to replace for instance with the quantity , where is very small, let’s say of the order . This is sufficient to eliminate all singularities occurring when the condition (24) is fulfilled while preserving the topological properties of as shown in Ref. GMMknotinvariant (). We have checked that the result of the calculation of is not much sensitive to the value of .

With the above setup, we have found that a number of Monte Carlo samples of the order of a few millions is sufficient to evaluate both contributions and to the knot invariant with a satisfactory precision for polymers of length up to . Smaller knots require a smaller amount of samples. To fix the ideas, with a knot of length , three million samples are enough to evaluate with a variance of about . For a knots of length , instead, with a number of samples of five millions the obtained variance is of the order of . This is not a bad result if we consider that for a knot with segments, the volume that is needed to be explored by the Monte Carlo sampling for the computation of is equal to . Supposing that we have a four-dimensional hypercube of such a volume, the length of its sides will be around lattice units. To evaluate with five millions samples means that each dimension is explored in the average only about times. In order to tune the computational time, the two most relevant parameters are the number of sampling points used in the Monte Carlo integration and the variance . and are related. If is too small, the error in the estimation of becomes large. In that case, the Monte Carlo evaluation of the integrals contained in is faster, but the rejection rate of the pivot transformations increases due to the high uncertainty on . On the other side, if is too big, the rejection rate of the pivot transformations decreases, but the time needed for the calculation of becomes unacceptably long. A good choice is to fix in such a way that the variance is approximately equal to . This is a safe estimation of the maximum possible error, because, as mentioned before, the minimum difference between two nearest non-coinciding values of is .

The above evaluations of the performance of the calculations have been made by assuming no particular action in order to improve the computational time. For instance, without any problem it is possible to reduce the size of the knot by a factor three by replacing every set of three contiguous segments with a single one. It is easy to realize that this does not alter the knot topology on a simple cubic lattice. Essentially, after this procedure a knot is obtained, whose length is one third of the length of the original knot. Of course, the new knot is defined off lattice. Another possibility to speed up the calculations consists in detecting particular elements of the knot that can be safely replaced by shorter ones. With these tricks the problem of studying the thermal properties of polymer knots with length up to becomes treatable.

## Iv The Wang-Landau method

The Wang-Landau (WL) method wanglandau () used here to compute the density of states has been already extensively discussed in the physical literature. Its convergence has been rigorously proven in zhoubhatt (). Here we limit ourselves to a brief review concerning the application of the WL algorithm to polymer knots.

Basically, the WL method is a self-adjusting procedure to compute the so-called density of states . For instance, let us consider the partition function

 Z=∫DXe−βH(X) (25)

of a system with Hamiltonian , where is a possible microstate of the system and denotes the Boltzmann factor in thermodynamic units, in which the Boltzmann constant is set to one: . Let us suppose that the admitted energy values are discrete with , so that may be rewritten in the form

 Z=∑ie−βEiϕi (26)

where

 ϕi=∫DXδEi,H(X) (27)

To compute the s, the WL algorithm proceeds as follows. Let denote the would be density of states and the energy histogram. At the zeroth approximation, we put:

 g(0)(Ei)=1, M(Ei)=0 (28)

Successively, a Markov chain of microstates is generated. In our case, the ’s differ from each other by transformations in which an element of the knot’s trajectory of length is changed by using the so-called pivot moves. We use a set of three possible pivot moves, called the inversion, reflection and interchange transformations. They have been discussed in Ref. pivot (), which the interested reader may consult for further details on this subject. Both the location of the element of the knot to be transformed and the kind of pivot transformation to be applied, are randomly selected. Also the number is chosen randomly between a given interval.

The probability of transition from a microstate of energy to a microstate of energy is given by:

 p(i→i′)=min[1,g(0)(Ei)g(0)(Ei′)] (29)

The microstate is accepted only if , being a randomly generated number in the interval . If the condition is not satisfied, the old microstate is accepted once again. In both cases, once a new microstate with energy has been selected with or , the corresponding would be density of states and the energy histogram are updated as shown below:

 g(0)(Ej) = f0g(0)(Ej), (30) M(Ej) = M(Ej)+1 (31)

where . Here we put . We remark that Eq. (30) modifies the probability that microstates of energy are accepted. In fact, the next time in which a microstate of this kind will randomly appear after a pivot transformation, its probability to be selected by the rule of Eq. (29) will be damped by a factor . This procedure of sampling new microstates that are chosen or rejected according to the transition probability (29) continues until the energy histogram becomes flat. Since in a real simulation it is nearly impossible to obtain a completely flat histogram, a deviation of no more than 20% of the s from their average value is admitted. Clearly, in order to have an almost flat energy histogram, the microstates corresponding to different energies should be almost equiprobable. In other words, after the WL procedure is completed, the probability of the occurrence of microstates with energy is a constant independent of . Let’s call this probability the WL probability and denote it with the symbol . To relate with the density of states , we remember that microstates are randomly generated with the help of pivot transformations. Thus, the WL probability must be equal to the unbiased probability of obtaining a microstate of energy by pivot transformations times the damping factor computed using the WL algorithm explained before. In formulas:

 PWL(Ei)=g(0)(Ei))−1×Punbiased(Ei) (32)

At this point, we assume that the unbiased probability is proportional to the density of states . More precisely:

 Punbiased(Ei)=ϕi∑jϕj (33)

where the sum over all possible energies is an irrelevant constant. The above relation is intuitive, because the greater is the density of microstates for a given value of the energy , the greater is the probability to obtain one of such microstates by random transformations. Substituting Eq. (33) in Eq. (32), we arrive at the desired result:

 PWL(Ei)=g(0)(Ei))−1×ϕi∑jϕj (34)

If the energy histogram is flat, also becomes a constant, so that it is possible to write up to irrelevant constants

 ϕi=g(0)(Ei) (35)

Since the s are delivered by the WL algorithm, also the density of states is known.

Actually, if is too big, the statistical errors on the ’s may grow large and the above equation is satisfied very roughly. On the other side, if is too small, it is necessary an enormous number of microstates during the sampling in order to derive the ’s. For this reason, in the WL procedure the density of states is computed by successive approximations. Let us introduce to this purpose the modification factors , with and . At the beginning of the th approximation, the value of the factor is decreased using the relation:

 fν=√fν−1 (36)

Moreover, the would be density of states is initialized in such a way that it coincides with the density of states obtained from the th approximation:

 g(ν)(Ei)=g(ν−1)(Ei) (37)

Finally, the energy histogram is set to zero. At this point, the would be density of states is computed at the next order by generating new microstates and applying the same procedure used above to evaluate . One should proceed in this way until, for some integer , the modification factor becomes sufficiently small, according to the original article wanglandau () and the changes in the s become statistically irrelevant.

To conclude this Section, a digression on the ergodicity of the pivot moves used here is in order. In the work pivot (), this ergodicity has been proved on a cubic lattice for dimensional self avoiding walks (SAW) with fixed ends, including those in which the ends are located at the distance of one lattice size and thus may be considered as closed. More precisely, it has been verified in pivot () that, starting from an arbitrary SAW of length in dimensions, it is possible to reduce it by a finite number of pivot moves to a given canonical SAW of the same length. Of course, no special restriction has been required on these moves concerning their ability to preserve the topology of a knot. Thus, some of those moves may allow the crossings of the lines of the SAW, a fact that can potentially destroy the topology of the knot. For that reason, the proof of the ergodicity of the pivot moves defined in pivot () is not sufficient in our case, in which the topology of the studied polymer knot is fixed since the beginning. A complete proof would require to show that it is possible to reach, starting from an arbitrary knot configuration, a given seed configuration of that knot by applying successive pivot moves. Moreover, the knot obtained after each move should have the same topology of the initial one. Despite the fact that we did not succeed to obtain such a proof up to now, our empirical investigations, performed on various knot configurations of different lengths, seem to suggest that the pivot moves of Ref. yzff () are ergodic also for real polymer knots, in which the trajectories cannot cross themselves and thus the topology is preserved. Indeed, in all analyzed cases, it has been possible to conclude that with the pivot moves of Ref. pivot (), together with the WL algorithm, conformations with every possible number of contacts can be visited after a sufficiently long run. A precise definition of contacts will be provided in the next Section after introducing the potential of the short-range forces that will take into account of the interactions between the monomers.

## V Thermal properties of knots

### v.1 Observables

In this Section, we study the thermal properties of the polymer knots , and . In particular, the specific energy , the heat capacity and the gyration radius will be computed. Both cases of attractive and repulsive short-range interactions, will be considered. To this purpose we introduce the following potential between two monomers:

 VIJ=⎧⎪⎨⎪⎩+∞if I=Jεif d=|→RI−→RJ|=1 and I≠J±10otherwise (38)

with being a small energy scale determining the strength interactions between pairs of non-bonded monomers. The condition characterizes the attractive case, while characterizes the repulsive case. Moreover, denotes the position vector of the -th segment. Conformations such that are automatically discarded, because no crossing of the trajectories is possible in a real polymer knot. Besides, when a crossing occurs, the knot is no longer mathematically well defined.

The total Hamiltonian of a polymer knot in a given microstate is given by:

 H(X)=L∑I,J=1VIJ (39)

To simplify the expression of , it is convenient to classify the microstates according to their number of contacts . Two non-bonded monomers are said to form a contact if their reciprocal distance on the lattice is equal to (see the second condition in Eq. (38)). counts the number of contacts of every pair of noncontiguous monomers appearing in the conformation. Clearly, for a microstate with a number of contacts and with no overlapping monomers, the Hamiltonian reads as follows:

 H(Xm)=mε (40)

The density of states defined by Eq. (27) is computed by means of the Wang-Landau algorithm illustrated in Section IV.

In our settings, and are expressed as follows:

 ⟨E(β)⟩ = ∑mmεe−βmεϕm∑me−βmεϕm (41) C(β) = 1T2(⟨E(β)2⟩−⟨E(β)⟩2) (42)

while the mean square radius of gyration is given by yf (); pnv ():

 ⟨R2G⟩(β)=∑m⟨R2G⟩me−βmεϕm∑me−βmεϕm (43)

with denoting the average of the gyration radius computed over states with contacts.

Finally, we wish to discuss how the problem of rare events has been treated in this work. While there is enough evidence that the pivot algorithm is ergodic, it is also true that certain states occur very rarely and sometimes may require the computations of billions of trial conformations before being obtained. This is the case of very compact () or very swollen () conformations. For polymers of any length , at least up to , the maximal length that has been studied with the present approach, states with a number of contacts higher than should be considered as rare events. For short polymer knots, samples with or more are extremely rare. When or , the lowest energy states that have been reached after a few millions of samples are characterized respectively by and contacts. Short polymers which are heavily knotted are also difficult to be found in conformations with a small number of contacts. For a trefoil with or the state with should be considered as rare, while for is rare. The problem is that rare events act like bottlenecks in the WL algorithm, which prevent the energy histogram to become flat unless a prohibitively high number of samples is used. Such bottlenecks appearing in polymer simulations have been discussed in binder (). In that reference, the introduction of a cutoff in the energy and a slightly different criterion of flatness has been proposed. Accordingly, our simulations have been restricted to conformations in which the number of contacts is limited by the condition . For short and topologically complex knots, has additionally been bounded from below to be greater than . In the case of some knots, which are not considered here, the restriction has been required. These cuts in the number of contacts affect the calculations at lower and higher temperatures at most by a relative error of a few percent. The data about the average gyration radius show for example that the difference between the mean square gyration radii obtained by cutting or not the five lowest values of is very small. On the other side, without imposing bounds on , we have observed that the flattening procedure of the energy histogram in the Wang-Landau algorithm requires an enormous amount of sampling, especially if the rare events occur when the modification factor is small.

### v.2 Simulation results based on the knot invariant ϱ(C)

In Fig. 2, left panel, the results for the specific energy and the heat capacity in the attractive case are displayed. Let us note that in all figures we have used the symbol for the normalized temperature . As it is possible to see, the specific energy increases with increasing temperatures as it should be, because as the temperature grows, more and more energetic states are excited. Moreover, at high temperatures longer polymers have a higher specific energy than shorter ones. This behavior may be explained by the following two observations. The first observation, which is true for both short and long polymers, is that, when the temperature is very high, the energy gap of the potential defined in Eq. (38) is much smaller than the energy related to the thermal fluctuations. As a consequence, the polymer is supposed to be in a more swollen conformation than at low temperatures, where the interactions dominate over the thermal fluctuations. The second observation is that the effects of knotting are less and less important with increasing polymer lengths. To convince ourselves that this is the case, it is sufficient to imagine a knot which is localized in a small part of the polymer. If the polymer is long, the localized part will be not relevant in comparison with the rest of the polymer, which will behave more or less like a unknotted ring. Also numerical simulations confirm this trend, see for example yzff (). Moreover, our calculations show that the differences in the specific energy and heat capacity between two trefoils of lengths and are much less marked than the same differences between two trefoils of lengths and . Indeed, the data of the trefoil with have not been reported in Fig. 2(a) and (b) because it is difficult to distinguish them from those of the trefoil with . Combining the two observations above, one can conclude that, at high temperatures, the polymer will tend to swell, but the effect of the topological constraints will be to counteract this swelling process. On the other side, the topological constraints and their effects become less important when polymers are long. As an upshot, if the temperature and the knot topology are kept fixed, we expect that the average distance between the monomers and thus the specific energy of the knot will increase with increasing polymer length when . This explains the behavior of the specific energy of Fig. 2 at higher temperatures. The behavior of at lower temperatures should be taken with some care because, as already mentioned, the number of contacts has been limited by the condition , so that the lowest energy state cannot be reached. Always for that reason, the specific energy at zero temperature is slightly different for different lengths.

Concerning the heat capacity, we can see in Fig. 2, right panel, that there is only one sharp peak in the whole temperature region. The interpretation of this peak is somewhat difficult. It is certainly a pseudo phase transition caused by the fact that we are working on a lattice with a finite size system, as those discussed in Refs. wuest (); janke (). Similar pseudo phase transitions have been already observed in knots, see swetnam (). In janke () it has been stressed that, along with the advances in constructing high resolution equipment, such pseudo phase transition in real system become more and more important. It is likely that the peaks in the heat capacity correspond to the transition of the knot from a frozen crystallite state to an expanded state similar to what happens in the case of a single polymer chain discussed in binder (). The data concerning the gyration radius in Fig. 6, left panel, seems to confirm this hypothesis, because the gyration radius stats to grow abruptly more or less in the same range of temperatures in which the peak is appearing. This hypothesis could also explain the absence of a second peak or shoulder connected to the coil-globule transition. As it was discussed in binder (), in fact, if the range of the interactions is very short, then an open chain admits just two possible states, the crystallite and the expanded coil ones. For a knot, the situation is however more complicated. The visual analysis of the samples shows that the lowest energy conformations exhibit some kind of ordering, but probably a full ordering is forbidden by the topology of the knot. In conclusion, further analysis is necessary in order to identify the nature of these transitions.

Analogous considerations can also be made in the repulsive case displayed in Fig. 3. One difference is that low temperatures correspond now to the swollen state, while at high temperatures, when the energy barrier becomes negligible with respect to the energy carried by the thermal fluctuations, the monomers are allowed to get nearer and the average number of contacts is growing. As a consequence, the gyration radius for repulsive interactions decreases at higher temperatures as shown in Figs. 6(b) and 7(b). Moreover, while as expected the specific energy increases with increasing temperatures, it turns out that, contrarily to what happens when the forces are attractive, longer polymers have a lower specific energy than shorter ones at any temperature. To explain this decrease of the energy of longer polymers, we recall that the swelling of the polymer knot is hindered by the topological constraints, in such a way that more complex knot configurations correspond to more compact polymer conformations. As well, due to the fact that those parts of the polymer trajectory which are affected by the topological constraints will become less and less important when the polymer length increases, it is licit to conclude that longer polymers will in the average admit conformations that are more swollen with respect to those of shorter polymers as we argued in the attractive case. The difference is that, if the interactions are repulsive, more swollen conformations are less energetic, which explains why the specific energy of longer polymers is in the average lower than that of shorter polymers. The presence of sharp peaks in the heat capacities, see Figs. 3(b) and 5(b), has not a straightforward interpretation like those occurring when the interactions are attractive. Apparently, as argued in yzff (), at the knot is in one of the lowest energy conformations which are allowed. As the temperature increases, at the beginning the system is unable to pass to a more compact configuration unless . After that threshold, the knot more and more contacts are possible between the monomers unless saturation is reached.

To check the effects of the topology on the behavior of the polymer, we have tested different knot configurations of the same length. We present here the data of polymers with . We observe that, in the attractive case, the increasing of the knot complexity results in the decreasing of the specific energy and heat capacity at high temperatures, see Figs. 4(a) and (b). As one may expect from the previous discussion, in the repulsive case it is exactly the converse, i. e. the increasing of the knot complexity results in the increasing of the specific energy and heat capacity when the temperature grows, as shown in Fig 5(a) and (b).

Once again, this is due to the fact that, if the topology of the knot is more complex, the knot conformation contains more contacts and is more compact. To have more contacts implies that, in the average, the energy of each monomer and thus the specific energy is lower if the forces are attractive and higher if the forces are repulsive.

This scenario is confirmed by the data on the gyration radius of different knots with length in both the attractive and repulsive case, see Fig. 7. The gyration radii of the three different knots and satisfy the inequality independently if the interactions are attractive or repulsive.

Finally, as we previously mentioned, at higher temperatures the effects of the interactions should become negligible. Accordingly, even if the gyration radius exhibits a different behavior depending if the interactions are repulsive or attractive, in Fig. 8 we see that the values of the radius of gyration computed in these two different cases get closer and closer when the temperature is increasing. Eventually, they should converge if the temperature is high enough.

### v.3 Comparison with the PAEA method

The PAEA method allows to use topology changing pivot transformations in numerical simulations of polymer systems and prevents the accidental transition to another topology without making use of topological invariants. The idea of the PAEA method is based on the following observation. After a pivot transformation is performed on a randomly chosen element of length of a knot, is transformed into the new element . Since and have their two ends in common, because these ends are untouched by the pivot transformation, it is easy to realize that, together, and form a small loop or, if is large, a set of small loops of total length . The topology is preserved by checking whether the part of the knot unaffected by the pivot transformation crosses an arbitrary surface spanned around the small loop(s). If the crossing happens, then the trial conformation will be rejected and the knot undergoes another pivot transformation. Otherwise, the trial conformation is accepted and a new pivot transformation is applied to it. More details can be found in yzff ().

What is crucial for the success of the PAEA method is to classify all the possible small loops of segments that can be produced after a random pivot transformation and to construct suitable surfaces having these small loops as borders. When the length of the segments changed by the pivot transformations in greater than five, the number of all loops of this kind becomes large. With increasing values of , it becomes more and more difficult to construct the surfaces mentioned above. Unfortunately, this is a not problem that can be solved automatically by the computer using some algorithm. For this reason, up to now the PAEA algorithm has been developed only in the case . The advantage of the technique discussed in this work based on the knot invariant is that the number of segments involved in the pivot transformations may arbitrarily range within the interval . Larger values of change bigger portions of the knot and this leads to a faster equilibration process than the PAEA method. Moreover, we have observed that, as the length of the segments to be changed is increasing, the number of samples needed to get a flat energy histogram is decreasing. This fact is shown in Table 2, where as an example the ratio is displayed in the case of the knot with length . In the calculations with the PAEA method the pivot transformations were limited to , while in the calculations with the invariant . All the approximation degrees used in computing the density of states have been listed. As we can see from Table 2, the number of samples needed in the PAEA method to make the energy histogram flat is always larger than the corresponding number of samples necessary to distinguish the topology with the help of . The explanation of this increasing in the efficiency of the method based on the knot invariant with respect to the PAEA method is that with the invariant large pivot transformations are allowed and they are able to modify a relevant portion of the polymer. The larger is the number of segments affected by a pivot transformation, the greater is the difference between the numbers of contacts of the knot before and after the transformation. As a consequence, with a large pivot transformation it is possible to jump from conformations that have very different values of , thus accelerating considerably the exploration of the set of all possible conformations of a knot compatible with the given topological configuration.