A study of polymer knots using a simple knot invariant written consisting of multiple contour integrals
aa r X i v : . [ c ond - m a t . s o f t ] J u l A study of polymer knots using a simple knot invariant written consisting of multiplecontour integrals
Yani Zhao ∗ and Franco Ferrari † CASA* and Institute of Physics, University of Szczecin, Szczecin, Poland (Dated: April 11, 2018)In this work the thermodynamic properties of short polymer knots (up to 120 segments) definedon a simple cubic lattice are studied with the help of the Wang-Landau Monte Carlo algorithm. Thesampling 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 consid-ered. The properties of the specific energy, heat capacity and gyration radius of the knots 3 , and 5 are discussed. It is found that the heat capacity exhibits a sharp peak. If the interactionsare attractive, similar peaks have been observed also in single open chains and have been related tothe transition from a frozen crystallite state to an expanded coil state. Some other peculiarities ofthe behavior of the analyzed observables are presented, like for instance the increasing or decreasingof the knot specific energy at high temperatures with increasing polymer lengths depending if theinteractions are attractive or repulsive. Besides the investigation of the thermodynamics of polymerknots, the second goal of this paper is to introduce a method for distinguishing the topology of aknot based on a topological invariant which is in the form of multiple contour integrals and explicitlydepends on the physical trajectory of the knot. The chosen invariant, denoted here ̺ ( C ), is relatedto the second coefficient of the Conway polynomial. It has been first isolated from the amplitudesof a Chern-Simons field theory with gauge group SU ( N ). It is shown that this invariant is veryreliable in distinguishing the topology of polymer knots. One of the advantages of the proposedapproach is that it allows to reduce the number of samples needed by the Wang-Landau algorithm.Some solutions to speed up the calculations of ̺ ( C ) exploiting Monte Carlo integration techniquesare developed. PACS numbers:
I. INTRODUCTION
Long polymers are very likely to be found in the con-figuration of knots or links. The topological properties ofpolymers with closed conformations play indeed an im-portant role in physics, chemistry and biology. For thatreason, they are being actively investigated [1–26]. Aparticularly challenging problem is that of the statisticalmechanics of polymer knots. Up to now, a satisfactoryanalytical model exists only in the case of two polymerrings linked together [27, 28], but there is no analogousmodel for a knot despite many attempts, see for instance[30–33] for a review on this subject. Moreover, the scal-ing laws of the most important observables of polymerknots, like for instance the gyration radius, are still asubject of intense research [34–38].The main difficulty behind the treatment of polymerknots, as well as polymer links, is to distinguish thewealth of different topological configurations of such sys-tems. This problem arises in analytical models becauseit is necessary to impose topological constraints in or-der to avoid that the statistical fluctuations affect theinitial topological state. This is a physical requirement,dictated by the fact that, once a polymer knot or linkhas been formed, its topological state cannot be modified ∗ Electronic address: yanizhao@fermi.fiz.univ.szczecin.pl † Electronic address: ferrari@fermi.fiz.univ.szczecin.pl without breaking the covalent bonds holding together themonomers. Thus, unwanted changes of the topologicalconfiguration must be detected and rejected. In numer-ical simulations, instead, a widely used way to generatepolymer knots or links is to consider self avoiding walks(SAW’s) that, at a certain point, intersect themselvesforming closed conformations [12, 39–41]. The problem inthis case is to sort out the topological configurations thatare of interest from all the other configurations producedwithin this approach. This process of generating knotand links is not very efficient for our purposes, becausethe probability of formation of knots or links of a giventype from SAW’s is very low, see for instance Ref. [42–44].Alternatively, it is possible to start from a polymer ringwith a seed trajectory that is out of equilibrium, but al-ready in the desired topological configuration. Later, thesystem is equilibrated using the so-called pivot transfor-mations [45, 46]. There exist pivot transformations thatare automatically preserving the topological state of thesystem, see for example [47], however they are able tomodify only a very small part of the whole polymer, afact that considerably increases the time for reaching theequilibrium, especially for very long polymers. In [48]another method has been proposed, called the Pivot Al-gorithm and Excluded Area (PAEA) method, in whichmore general pivot transformations are considered, in-cluding those that could potentially change the topology.In this case, large transformations are not easy to be im-plemented.In conclusion, apart from a few exceptions, like for in-stance the already cited works of [47, 48] and notably thedynamic Monte Carlo approach [49], that however leadsto some level of polydispersity, in most numerical com-putations involving polymer knots, topological invariantsare exploited whenever it becomes necessary to distin-guish the topological configuration of the system underinvestigation. The most popular topological invariantsare given in the form of polynomials or of multiple con-tour integrals computed along the physical trajectories ofthe polymers. The latter invariants are easier to be usedin analytical models than the former, because the coeffi-cients of the polynomials are not directly related to thepolymer conformation. The particular simplicity of theGauss linking number, which consists in a double con-tour integral, is the main reason for which it has beenpossible to derive an analytical model of two linked poly-mer rings [27–29]. Moreover, the Gauss linking numberhas already been used in numerical simulations of poly-mer systems, see for instance [50, 51]. Unfortunately,in order to distinguish the topology of knots, there isno simple invariant like the Gauss linking number. Sofar, the statistical properties of polymer knots have beenstudied numerically with the help of topological invari-ants like the Alexander polynomials [12] or the HOMFLYpolynomials [52], see the detailed textbook by Kleinert[33] for an extensive review on that subject. Of course,there is a plenty of other knot invariants that have beenor could be applied, for instance the Conway polynomi-als [53], the Arf-Casson invariant [54] or the Milnor [55]and Vassiliev-Kontsevich invariants [56]. The latter threehave an explicit representation in terms of multiple inte-grals computed along the knot trajectories.The purpose of this paper is to show that knot in-variants that are given in the form of multiple contourintegrals can represent a valid alternative in numericalcalculations to polynomial knot invariants or even tothe PAEA method, which is able to detect the topol-ogy changes exactly and has been proven to allow veryfast computations of the statistical properties of poly-mer knots [48, 57]. In particular, we will concentrateon a topological invariant denoted here ̺ ( C ), where C denotes the trajectory of the knot. ̺ ( C ) has been de-rived from the one-loop amplitudes of non-abelian Chern-Simons field theories with gauge group SU ( N ) [54]. Itis related to the Arf-Casson invariant and to the secondcoefficient of the Conway polynomials [54]. Its value canbe analytically computed for any given knot configura-tion. ̺ ( C ) is the simplest knot invariant represented interms of multiple contour integrals.The knot invariant ̺ ( C ) is applied here in order toderive the average values of the specific energy, heat ca-pacity and gyration radius of several different knot con-figurations by means of the Wang-Landau algorithm [58].We find in this way that the examined knots, correspond-ing namely to the trefoil 3 , the figure-eight 4 and 5 [75],undergo with the temperature a phase transition, whichis probably from a frozen crystallite state to an expandedcoil state similarly to what happens in the case of a single polymer chain proposed in [59]. Other physical proper-ties of polymer knots are discussed. Compared with thePAEA method, the use of ̺ ( C ) allows to reduce the num-ber of samples necessary for the calculations of the aver-ages of the observables with the Wang-Landau algorithm.As it will be seen, this reduction is due to the fact thatwith ̺ ( C ), large pivot transformations can be exploitedwhich are able to change relevant portions of the knot.In this way, the exploration of the whole set of availableconformations becomes faster. Despite the decreasing ofthe number of samples, the computations last in generallonger than those performed with the PAEA method, be-cause the expression of ̺ ( C ) contains quadruple integralsthat should be evaluated numerically and this takes sometime. As a consequence, we present here the results forrelatively short polymer knots up to L = 120, where L is the number of segments composing the knot on asimple cubic lattice. Actually, there is no problem instudying longer polymer knots. The reason is that in-variants given in the form of contour integrals can becomputed for arbitrarily deformed knots, not necessarilydefined on a lattice, provided the topological configura-tion remains the same after the deformation. Thanks tothis fact, we have found that it is possible to shorten thenumber of segments by a factor three, speeding up thecalculation of ̺ ( C ) considerably. Even using this trick,the method requires too long times (exceeding a monthon a modern workstation) if L > ̺ ( C )becomes however competitive in the equilibration of verylong polymers because, as already mentioned, it offers thepossibility of exploring very fast the set of conformationscompatible with the topological constraints.The rest of the current work is organized as follows.Our simulation set-up and the topological invariant ̺ ( C )are introduced in Section II. The Monte Carlo integra-tion method we use to calculate the knot invariant ̺ ( C )is introduced in Section III. An explanation of the Wang-Landau algorithm with particular attention to its appli-cations to the statistical mechanics of polymers is pro-vided in Section IV. The results on the thermal propertiesof different polymer knots are presented in Section V. Wealso compare our calculations with the results obtainedwith the PAEA method in Refs. [48, 57]. Finally, inSection VI we draw our conclusions and possible gener-alizations of this work are briefly discussed. II. METHODOLOGYA. General outline of the used methodology
First of all, a brief digression on the used terminologyis in order. Throughout this work the word configurationrefers to a particular topological state of a polymer knot.The word conformation will instead denote the particularshape in the space of the trajectory of a polymer knot ina given topological configuration.At this point it is possible to go back to the outlineof the methodology. We adopt the strategy of consider-ing at the beginning a seed conformation of the knot tobe analyzed. The knot is a self-avoiding polygon definedon a simple cubic lattice with edges of unit length. Themonomers are located on the vertices of the lattice. Let L denote the total length of the polygon. Since it consistsof edges or segments of unit length, L coincides also withthe number of segments composing the polygon. Thecrossing of the knot trajectory with itself at some pointon the lattice is forbidden. The starting seed configura-tion is equilibrated and then used to compute the ther-modynamic properties of the studied knot. The relevantobservables are calculated by means of the Wang-LandauMonte Carlo algorithm [58], which will be explained insome more details later. Both cases of attractive and re-pulsive interactions are considered. For the equilibrationof the knot and the sampling in the Wang-Landau al-gorithm, the polymer conformations are randomly modi-fied by exploiting the pivot transformations described inRef. [45]. These transformations can involve any number N of segments such that 1 < N ≤ L and are not prevent-ing the crossing of the lines of the knot trajectory C . Forthat reason, the invariant ̺ ( C ) should be applied in or-der to check if the topology of the knot has been alteredafter each pivot transformation. If this is the case, thetransformation is rejected and a new one is considered. B. The topological invariant ̺ ( C ) To avoid topology changes of a polymer knot C po-tentially occurring after the pivot transformations, wewill use in this work a topological invariant that hasbeen derived from the one-loop amplitude of the Wilsonloop in non-abelian SU ( N ) Chern-Simons field theories [54, 60, 61]. The most important characteristic of thisinvariant, that will be denoted ̺ ( C ), is that it can beexpressed in the form of a sum of multiple contour inte-grals: ̺ ( C ) = ̺ ( C ) + ̺ ( C ) (1)where the knot C is represented as an oriented closedpath of length L . The contribution ̺ ( C ) is given by thetriple integral: ̺ ( C ) = − π I C dx µ Z x dy ν Z y dz ρ I µνρ ( ~x, ~y, ~z ) , (2)with I µνρ ( ~x, ~y, ~z ) = ǫ αβγ ǫ µασ ǫ νβλ ǫ ργτ × Z d ~ω ( ω − x ) σ | ~ω − ~x | ( ω − y ) λ | ~ω − ~y | ( ω − z ) τ | ~ω − ~z | (3)while the second part ̺ ( C ) is: ̺ ( C ) = 18 π I C dx µ Z x dy ν Z y dz ρ Z z dw σ ǫ σνα ǫ ρµβ × ( w − y ) α | ~w − ~y | ( z − x ) β | ~z − ~x | (4)In the above formulas greek letters denote space indexes.The variables x µ , y ν , z ρ and w σ , µ, ν, ρ, σ = 1 , ,
3, arethe components of the radius vectors describing the posi-tions of four points on the same curve C . ǫ µνρ representsinstead the completely antisymmetric tensor uniquely de-fined by the condition ǫ = 1. The integrations alongthe path C in (2) and (4) are path ordered. This can beseen explicitly by parametrizing the trajectory C withthe arc-length: ̺ ( C ) = − π Z L ds dx µ ( s ) ds Z s dt dy ν ( t ) dt Z t du dz ρ ( u ) du I µνρ ( ~x, ~y, ~z ) , (5)and ̺ ( C ) = 18 π Z L ds dx µ ( s ) ds Z s dt dy ν ( t ) dt Z t du dz ρ ( u ) du Z u dv dw σ ( v ) dv ǫ σνα ǫ ρµβ ( w ( v ) − y ( t )) α | ~w ( v ) − ~y ( t ) | ( z ( u ) − x ( s )) β | ~z ( u ) − ~x ( s ) | (6)It has been shown that the knot invariant appearingabove is related to the second coefficient a ( C ) of theConway polynomial of a knot C through the followingrelation [54]: a ( C ) = 12 (cid:20) ̺ ( C ) + 112 (cid:21) (7)The Conway polynomials are well known and their coeffi-cients can be computed analytically for every knot topol- ogy, so that thanks to Eq. (7) it is easy to derive also thevalues of ̺ ( C ). In Table I we give a list of the second coef-ficients of the Conway polynomial and the correspondingvalues of ̺ ( C ) for the knot configurations that will bestudied here. ̺ ( C ) is the simplest known knot invariant that can beexpressed in the form of contour integrals. Like any otherknot invariant, ̺ ( C ) is not able to distinguish differentknots unambiguously. For example, the trefoil knot 3 TABLE I: This table provides the values of the second coef-ficients of the Conway polynomials and of the correspondingtopological invariants for the trefoil 3 , the figure-eight 4 andthe knot 5 . knot type a ( C ) ̺ ( C )3 -1 − has ̺ ( C ) = , exactly the same value of the knots6 , , and many others. However, we should keepin mind that the main role of a knot invariant in study-ing the thermal and mechanical properties of polymerknots is not to guess its topological configuration. Infact, the topological configuration is known since the be-ginning. The problem is rather to preserve that configu-ration against thermal fluctuations, because without anyconstraint the polymer trajectories are allowed to crossthemselves, a fact that can potentially alter a knot. Theprobability that due to thermal fluctuations a polymerring jumps from one knot configuration to another withthe same value of ̺ ( C ) is very low, as it has been observedin our simulations. What emerges from them is that it isvery unlikely that a knot passes to another configurationwith the same value of ̺ ( C ) after a pivot transforma-tion. Most probably, one ends up with the trivial knotor, somewhat less frequently, with a conformation withlower number of crossings C ′ such that ̺ ( C ) = ̺ ( C ′ ). Forthe purposes of this work, it is thus possible to affirm that ̺ ( C ) is a powerful knot invariant. To convince oneself,it is sufficient to recall that, if one considers the simplestknots up to ten crossings, there are particular topologicalconfigurations that are uniquely distinguished by ̺ ( C ),like 9 and 10 or are very efficiently distinguished fromall the others because their corresponding values of ̺ ( C )occur rarely. Luckily, if the values of ̺ ( C ) for two topo-logically different knots are not the same, the smallestdifference between them is 2. As an example, the knot5 has ̺ ( C ) = , while for the topologically inequivalentknots 5 and 9 we have ̺ ( C ) = and ̺ ( C ) = re-spectively. This allows to choose the number of samplingpoints in such a way that the variance in the Monte Carloevaluation of ̺ ( C ) is low enough that the probability ofconfusing two different knot topologies due to numericalerrors is negligible.The price to be paid for this efficiency in distinguishingknots is the complicated expression of ̺ ( C ). The mosttime-consuming contribution to ̺ ( C ) is the quadruplecontour integral necessary to compute ̺ ( C ) in Eq. (4).For a knot of length L , the evaluation time of ̺ ( C ) scalesas L . This is approximately one order more than thetime necessary to evaluate the Alexander polynomial ofa knot [12, 62, 63], which scales as ( M − . Here M denotes the number of crossings which is necessary torepresent the knot by projecting it on an arbitrary plane, see [63] for more details. Of course, also the computa-tion of the Alexander polynomial becomes prohibitive forpolymers which are long or have compact conformations,because in these cases the number of crossings M drasti-cally increases [14]. Moreover, the scaling law ( M − of the computational time is true only if the determinantof a M × M matrix that arises in the algorithm for com-puting the Alexander polynomial is evaluated with themethod of Gaussian elimination, which is subjected onround-off error that become important when M is large.One advantage of ̺ ( C ) is that its calculation can beextended without any effort to any kind of trajectory, notnecessarily on a cubic lattice. This fact is very helpfulwhen the polymer is long, so that it is advisable to de-crease the number L of its segments. In a very simpleway it is possible to reduce L by a factor three withoutdestroying the topology by replacing in the knot everygroup of three contiguous segments with a single seg-ment. As well, there is no problem in distinguishing thechanges of topology when the number N involved in thepivot transformations becomes large.Like the Alexander polynomial, also ̺ ( C ) is not ableto distinguish uniquely two different topological config-urations and is subjected to numerical errors. However,we have seen above that, for the goals of this work, theinvariant ̺ ( C ) is powerful enough. III. MONTE CARLO EVALUATION OF PATHORDERED CONTOUR INTEGRALS ON ALATTICEA. Simpson’s rule vs. Monte Carlo method
First of all, we introduce some notation that will beuseful in this Section. The contour C describing the phys-ical trajectory of the knot in space is represented here asa curve ~x ( s ), with 0 ≤ s ≤ L . Of course, C is consistingof a set of discrete segments, see Fig. 1 for an examplewith L = 10, and this fact should be taken into account.To this purpose, let’s denote with ~x i , i = 1 , · · · , L , thelocations of the lattice sites through which the closed con-tour C is passing. The i -th segment of the loop C formsa vector ~x i +1 − ~x i for i = 1 , · · · , L −
1. The L -th segmentis instead associated with the vector ~x − ~x L . Next, let ~x i (˜ s ), with i = 1 , · · · , L , be the restriction of the curve ~x ( s ) to the i − th segment. Here we have introduced thesegment’s arc-length ˜ s such that 0 ≤ ˜ s ≤
1. On the i − thsegment, ˜ s is related to s by the formula˜ s = s − i for i = 1 , . . . , L − s = s for i = L (9)Explicitly, the expression of x i (˜ s ) is given by: ~x i ( s ) = ~x i + ˜ s ( ~x i +1 − ~x i ) i = 1 , . . . , L − ~x L ( s ) = ~x L + ˜ s ( ~x − ~x L ) (11)From the above equations, it is easy to derive also thederivative of ~x i (˜ s ) with respect to ˜ s . Whenever it will be FIG. 1: This figure illustrates the notation used in Section IIIto describe closed contours on a simple cubic lattice with thehelp of the example of a short contour of length L = 10.The coordinate of the point p , represented in the figure by anempty circle, is ~x (˜ s p ) = ( x (˜ s p ) , x (˜ s p ) , x (˜ s p )), where ˜ s p isthe distance of the point p from the lattice site ~x . necessary to specify points on different elements of theloop C , for instance on segments i, j, k, l, . . . , they willbe denoted with the symbols ~x i (˜ s ) , ~y j (˜ t ) , ~z k (˜ u ) , ~w l (˜ v ) , . . . ,where i, j, k, l = 1 , · · · , L and 0 ≤ s, t, u, v ≤ ̺ ( C ) and ̺ ( C ) displayed in Eqs. (5) and (6) respec-tively in a form that is suitable for applying the standardSimpson’s rule: ̺ ( C ) = − π L X i =1 i X j =1 j X k =1 Z d ˜ s dx µi (˜ s ) d ˜ s Z − δ ij (1 − ˜ s )0 d ˜ t dy νj (˜ t ) d ˜ t × Z − δ jk (1 − ˜ t )0 d ˜ u dz ρk ( u ) d ˜ u I µνρ ( ~x i (˜ s ) , ~y j (˜ t ) , ~z k (˜ u )) (12)and ̺ ( C ) = 18 π L X i =1 i X j =1 j X k =1 k X l =1 Z d ˜ s dx µi (˜ s ) d ˜ s Z − δ ij (1 − ˜ s )0 d ˜ t dy νj (˜ t ) d ˜ t Z − δ jk (1 − ˜ t )0 d ˜ u dz ρk (˜ u ) d ˜ u × Z − δ kl (1 − ˜ u )0 d ˜ v dw σl (˜ v ) d ˜ v ǫ σνα ǫ ρµβ ( w l (˜ v ) − y j (˜ t )) α | ~w l (˜ v ) − ~y j (˜ t ) | ( z k (˜ u ) − x i (˜ s )) β | ~z k (˜ u ) − ~x i (˜ s ) | (13)The boundaries in the integrals of ̺ ( C ) and ̺ ( C ) havebeen chosen in such a way that the path ordering of thetrajectories is preserved. For example, if the integralsover ˜ s and ˜ t are performed over two different segments,then the boundaries for both variables are ranging be-tween 0 and 1. Instead, if the integrals are performed onthe same segment i = j , then it must be that 0 ≤ ˜ s ≤ ≤ ˜ t ≤ ˜ s . Let us note that the integrals over thevariable ˜ v in Eq. (13) can be performed exactly. Theremaining integrals over ˜ s, ˜ t and ˜ u should be evaluatednumerically, for instance by means of the Simpson’s rule.It turns out that for the evaluation of these integrals, thezeroth order approximation, obtained by substitutions of the kind: Z − δ ij (1 − ˜ s )0 dy νj (˜ t ) f ν ( ~y j (˜ t )) ∼ ( y νj (1) − y νj (0)) f ν ( ~y j (1)) + f ν ( ~y j (0))2 (14)is sufficient and gives satisfactory results. The prob-lem is that, even exploiting the crude approximationsof Eq. (14) in order to distinguish the topology changesof a given polymer knot, still a sum of L terms shouldbe evaluated in order to derive the value of ̺ ( C ) fromEq. (13). Already in the case of polymers of length L = 100 or more, this number becomes prohibitivelyhigh for practical purposes. In fact, the Wang-Landauprocedure used to compute the density of states requiresseveral millions of samples to be evaluated. For thatreason, it is much better to estimate ̺ ( C ) and ̺ ( C )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 s, t, u and v : ̺ ( C ) = Z L ds Z s dt Z t duF ( s, t, u ) (15)and ̺ ( C ) = Z L ds Z s dt Z t du Z u dvF ( s, t, u, v ) (16)where F ( s, t, u ) = − π dx µ ( s ) ds dy ν ( t ) dt dz ρ ( u ) du × I µ,ν,ρ ( ~x ( s ) , ~y ( t ) , ~z ( u )) (17) and F ( s, t, u, v ) = 18 π dx µ ( s ) ds dy ν ( t ) dt dz ρ ( u ) du dw σ ( v ) dv × ǫ σνα ǫ ρµβ ( w ( v ) − y ( t )) α | ~w ( v ) − ~y ( t ) | ( z ( u ) − x ( s )) β | ~z ( u ) − ~x ( s ) | (18)The variables s, t and u in Eq. (15) span a space of vol-ume V = L , while the variables s, t, u and v in Eq. (16)span a space of volume V = L . To evaluate the righthand sides of Eqs. (15) and (16) via Monte Carlo inte-gration, we can exploit the general formula Z b a dξ Z ξ a dξ · · · Z ξ m − a m dξ m f ( ξ , · · · , ξ m ) ≈ N " N X i =1 f ( ξ ( i )1 , · · · , ξ ( i ) m )( b − a ) m Y σ =2 ( ξ ( i ) σ − a σ ) (19)where the ξ ( i ) σ ’s, i = 1 , · · · , N and σ = 1 , · · · , m denoterandomly chosen variables in the range:[ a , b ] when σ = 1[ a σ , ξ σ ] when σ = 2 , . . . , m (20)In the numerical evaluation of ̺ ( C ), we consider thetrajectory C of the knot, which in principle is a polygonon a simple cubic lattice, as a continuous curve ~x ( s ). Fora given value of s , the segment on which the point ~x ( s )is located is identified by the relation i = [ s ] + 1 (21)where [ s ] denotes the integer part of s . The componentsof the curve ~x ( s ) are obtained using its restriction ~x i (˜ s ) tothe i − th segment, whose components can be computedexploiting Eqs. (10) and (11). The components of ~y ( t ), ~z ( u ) and ~w ( v ) are derived analogously.We stress the fact that the integrands F ( s, t, u ) and F ( s, t, u, v ) are regular, even if for instance F ( s, t, u, v )seems to be divergent when ~w ( v ) = ~y ( t ) or ~z ( u ) = ~x ( s ).Analytically, it is possible to prove that these singular-ities cancel as expected in a topological invariant. Innumerical computations the situation looks however dif-ferent, because one has to cope with terms that are sep-arately diverging, but whose sum is finite. A regulariza-tion is thus necessary in order to eliminate these ambi-guities. Let us first consider the computation of ̺ ( C ).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 situ-ations by choosing randomly the variable s, t, u, v is very low. In fact, if the calculations are performed using dou-ble precision variables, the number of digits after thefloating point is so high, that in practice the conditionsdisplayed in Eq. (22) are never satisfied. More seriousis the case of ̺ ( C ). After the integration over ~ω inEq. (3), we get the explicit expression of F ( s, t, u ) whichis reported in Appendix A. We see that, besides the sin-gularities 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 therandom sampling, the divergences of Eq. (23) will ap-pear due to the same reasons explained in the case of theanalogous divergences in Eq. (22), the condition (24) isvery easy to be realized. It is sufficient for example thatthe two vectors ~y ( t ) − ~x ( s ) and ~z ( u ) − ~x ( s ) have oppositeorientations and two zero components. Depending on thetopology of the knot, this situation may occur rather of-ten. Even if the number of cases in which this happensis negligible with respect to the total number of sampledpoints generated during the Monte Carlo computation of ̺ ( C ), still the result may be spoiled by overflow or un-derflow errors. To cure the singularity in Eq. (24) whenthe quantity A is equal to zero, we use the framing regu-larization introduced in [64]. This consists in performingan almost infinitesimal shift of the curves ~x ( s ), ~y ( t ), ~z ( u )and ~w ( v ) along a direction which is normal to the knot C . By denoting with ~n ( s ) the unit vector that gives thenormal direction to C , this means to replace for instance ~x ( s ) with the quantity ~x ( s )+ ǫ~n ( s ), where ǫ is very small,let’s say of the order ǫ ∼ − . This is sufficient to elim-inate all singularities occurring when the condition (24)is fulfilled while preserving the topological properties of ̺ ( C ) as shown in Ref. [54]. We have checked that theresult of the calculation of ̺ ( C ) is not much sensitive tothe value of ǫ .With the above setup, we have found that a number ofMonte Carlo samples of the order of a few millions is suf-ficient to evaluate both contributions ̺ ( C ) and ̺ ( C ) tothe knot invariant ̺ ( C ) with a satisfactory precision forpolymers of length up to L = 125. Smaller knots requirea smaller amount of samples. To fix the ideas, with aknot of length L = 60, three million samples are enoughto evaluate ̺ ( C ) with a variance of about 0 .
2. For a knotsof length 125, instead, with a number of samples of fivemillions the obtained variance is of the order of 1 .
2. Thisis not a bad result if we consider that for a knot with 125segments, the volume that is needed to be explored bythe Monte Carlo sampling for the computation of ̺ ( C ) isequal to ∼ . × . Supposing that we have a four-dimensional hypercube of such a volume, the length of itssides will be around 56 lattice units. To evaluate ̺ ( C )with five millions samples means that each dimension isexplored in the average only about 47 times. In orderto tune the computational time, the two most relevantparameters are the number of sampling points Υ used inthe Monte Carlo integration and the variance Σ. Υ andΣ are related. If Υ is too small, the error in the estima-tion of ̺ ( C ) becomes large. In that case, the Monte Carloevaluation of the integrals contained in ̺ ( C ) is faster, butthe rejection rate of the pivot transformations increasesdue to the high uncertainty on ̺ ( C ). On the other side,if Υ is too big, the rejection rate of the pivot transforma-tions decreases, but the time needed for the calculationof ̺ ( C ) becomes unacceptably long. A good choice is tofix Υ in such a way that the variance is approximatelyequal to 1. This is a safe estimation of the maximum pos-sible error, because, as mentioned before, the minimumdifference between two nearest non-coinciding values of ̺ ( C ) is 2.The above evaluations of the performance of the cal-culations have been made by assuming no particular ac-tion in order to improve the computational time. Forinstance, without any problem it is possible to reducethe size of the knot by a factor three by replacing everyset of three contiguous segments with a single one. It iseasy to realize that this does not alter the knot topologyon a simple cubic lattice. Essentially, after this procedurea knot is obtained, whose length is one third of the lengthof the original knot. Of course, the new knot is definedoff lattice. Another possibility to speed up the calcula-tions consists in detecting particular elements of the knotthat can be safely replaced by shorter ones. With thesetricks the problem of studying the thermal properties ofpolymer knots with length up to L ∼
400 becomes treat-able.
IV. THE WANG-LANDAU METHOD
The Wang-Landau (WL) method [58] used here tocompute the density of states has been already exten-sively discussed in the physical literature. Its conver-gence has been rigorously proven in [65]. Here we limitourselves to a brief review concerning the application ofthe WL algorithm to polymer knots.Basically, the WL method is a self-adjusting proce-dure to compute the so-called density of states φ i . Forinstance, let us consider the partition function Z = Z D Xe − βH ( X ) (25)of a system with Hamiltonian H ( X ), where X is a pos-sible microstate of the system and β = T denotes theBoltzmann factor in thermodynamic units, in which theBoltzmann constant is set to one: k B = 1. Let us sup-pose that the admitted energy values E i are discrete with i = 0 , , . . . , so that Z may be rewritten in the form Z = X i e − βE i φ i (26)where φ i = Z D Xδ E i ,H ( X ) (27)To compute the φ ′ i s, the WL algorithm proceeds as fol-lows. Let g ( E i ) denote the would be density of statesand M ( E i ) the energy histogram. At the zeroth approx-imation, we put: g (0) ( E i ) = 1 , M ( E i ) = 0 (28)Successively, a Markov chain of microstates X (1) , X (2) , X (3) , . . . is generated. In our case, the X ( i ) ’s differ from each other by transformations inwhich an element of the knot’s trajectory of length N ischanged by using the so-called pivot moves. We use aset of three possible pivot moves, called the inversion,reflection and interchange transformations. They havebeen discussed in Ref. [45], which the interested readermay consult for further details on this subject. Both thelocation of the element of the knot to be transformedand the kind of pivot transformation to be applied,are randomly selected. Also the number N is chosenrandomly between a given interval.The probability of transition from a microstate X i ofenergy E i to a microstate X i ′ of energy E i ′ is given by: p ( i → i ′ ) = min (cid:20) , g (0) ( E i ) g (0) ( E i ′ ) (cid:21) (29)The microstate X i ′ is accepted only if p ( i → i ′ ) ≥ η , η being a randomly generated number in the interval [0 , p ( i → i ′ ) ≥ η is not satisfied, the oldmicrostate X i is accepted once again. In both cases, oncea new microstate with energy E j has been selected with j = i ′ or j = i , the corresponding would be density ofstates g (0) ( E j ) and the energy histogram are updated asshown below: g (0) ( E j ) = f g (0) ( E j ) , (30) M ( E j ) = M ( E j ) + 1 (31)where f >
1. Here we put f = e . We remark thatEq. (30) modifies the probability that microstates of en-ergy E j are accepted. In fact, the next time in which amicrostate of this kind will randomly appear after a pivottransformation, its probability to be selected by the ruleof Eq. (29) will be damped by a factor f . This procedureof sampling new microstates that are chosen or rejectedaccording to the transition probability (29) continues un-til the energy histogram becomes flat. Since in a realsimulation it is nearly impossible to obtain a completelyflat histogram, a deviation of no more than 20% of the M ( E i ) ′ s from their average value is admitted. Clearly,in order to have an almost flat energy histogram, the mi-crostates corresponding to different energies E i should bealmost equiprobable. In other words, after the WL pro-cedure is completed, the probability of the occurrence ofmicrostates with energy E i is a constant independent of i .Let’s call this probability the WL probability and denoteit with the symbol P W L ( E i ). To relate P W L ( E i ) with thedensity of states φ i , we remember that microstates arerandomly generated with the help of pivot transforma-tions. Thus, the WL probability P W L ( E i ) must be equalto the unbiased probability P unbiased ( E i ) of obtaining amicrostate of energy E i by pivot transformations timesthe damping factor ( g (0) ( E i )) − computed using the WLalgorithm explained before. In formulas: P W L ( E i ) = g (0) ( E i )) − × P unbiased ( E i ) (32)At this point, we assume that the unbiased probability P unbiased ( E i ) is proportional to the density of states φ i .More precisely: P unbiased ( E i ) = φ i P j φ j (33)where the sum over all possible energies P j φ j is an irrel-evant constant. The above relation is intuitive, becausethe greater is the density of microstates for a given valueof the energy E i , the greater is the probability to ob-tain one of such microstates by random transformations.Substituting Eq. (33) in Eq. (32), we arrive at the desiredresult: P W L ( E i ) = g (0) ( E i )) − × φ i P j φ j (34)If the energy histogram is flat, also P W L ( E i ) becomes aconstant, so that it is possible to write up to irrelevantconstants φ i = g (0) ( E i ) (35) Since the g (0) ( E i ) ′ s are delivered by the WL algorithm,also the density of states φ i is known.Actually, if f is too big, the statistical errors on the g (0) ( E i )’s may grow large and the above equation is sat-isfied very roughly. On the other side, if f is too small,it is necessary an enormous number of microstates dur-ing the sampling in order to derive the g (0) ( E i )’s. Forthis reason, in the WL procedure the density of statesis computed by successive approximations. Let us intro-duce to this purpose the modification factors f ν , with ν = 0 , , . . . and f = e . At the beginning of the ν − thapproximation, the value of the factor f ν − is decreasedusing the relation: f ν = p f ν − (36)Moreover, the would be density of states g ( ν ) ( E i ) is ini-tialized in such a way that it coincides with the density ofstates g ( ν − ( E i ) obtained from the ( ν − − th approxi-mation: g ( ν ) ( E i ) = g ( ν − ( E i ) (37)Finally, the energy histogram is set to zero. At this point,the would be density of states g ( ν ) ( E i ) is computed at thenext order by generating new microstates and applyingthe same procedure used above to evaluate g (0) ( E i ). Oneshould proceed in this way until, for some integer ¯ ν , themodification factor f ¯ ν becomes sufficiently small, f ¯ ν ∼ · − according to the original article [58] and the changesin the g ( ν )( E i ) ′ s become statistically irrelevant.To conclude this Section, a digression on the ergod-icity of the pivot moves used here is in order. In thework [45], this ergodicity has been proved on a cubic lat-tice for d − dimensional self avoiding walks (SAW) withfixed ends, including those in which the ends are locatedat the distance of one lattice size and thus may be con-sidered as closed. More precisely, it has been verified in[45] that, starting from an arbitrary SAW of length L in d − dimensions, it is possible to reduce it by a finite num-ber of pivot moves to a given canonical SAW of the samelength. Of course, no special restriction has been requiredon these moves concerning their ability to preserve thetopology of a knot. Thus, some of those moves may allowthe crossings of the lines of the SAW, a fact that can po-tentially destroy the topology of the knot. For that rea-son, the proof of the ergodicity of the pivot moves definedin [45] is not sufficient in our case, in which the topologyof the studied polymer knot is fixed since the beginning.A complete proof would require to show that it is possibleto reach, starting from an arbitrary knot configuration, agiven seed configuration of that knot by applying succes-sive pivot moves. Moreover, the knot obtained after eachmove should have the same topology of the initial one.Despite the fact that we did not succeed to obtain such aproof up to now, our empirical investigations, performedon various knot configurations of different lengths, seemto suggest that the pivot moves of Ref. [48] are ergodicalso for real polymer knots, in which the trajectories can-not cross themselves and thus the topology is preserved.Indeed, in all analyzed cases, it has been possible to con-clude that with the pivot moves of Ref. [45], togetherwith the WL algorithm, conformations with every possi-ble number of contacts can be visited after a sufficientlylong run. A precise definition of contacts will be pro-vided in the next Section after introducing the potentialof the short-range forces that will take into account ofthe interactions between the monomers. V. THERMAL PROPERTIES OF KNOTSA. Observables
In this Section, we study the thermal properties of thepolymer knots 3 , 4 and 5 . In particular, the specificenergy h E ( β ) i , the heat capacity C ( β ) and the gyrationradius h R G i ( β ) will be computed. Both cases of attrac-tive and repulsive short-range interactions, will be con-sidered. To this purpose we introduce the following po-tential between two monomers: V IJ = + ∞ if I = Jε if d = | ~R I − ~R J | = 1 and I = J ±
10 otherwise (38)with ε being a small energy scale determining thestrength interactions between pairs of non-bondedmonomers. The condition ε < ε > ~R I denotes the position vector of the I -th seg-ment. Conformations such that I = J are automaticallydiscarded, because no crossing of the trajectories is pos-sible in a real polymer knot. Besides, when a crossingoccurs, the knot is no longer mathematically well defined.The total Hamiltonian of a polymer knot in a givenmicrostate X is given by: H ( X ) = L X I,J =1 V IJ (39)To simplify the expression of H ( X ), it is convenient toclassify the microstates according to their number of con-tacts m . Two non-bonded monomers are said to form acontact if their reciprocal distance on the lattice is equalto 1 (see the second condition in Eq. (38)). m countsthe number of contacts of every pair of noncontiguousmonomers appearing in the conformation. Clearly, for amicrostate X m with a number m of contacts and with nooverlapping monomers, the Hamiltonian reads as follows: H ( X m ) = mε (40)The density of states φ m defined by Eq. (27) is computedby means of the Wang-Landau algorithm illustrated inSection IV. In our settings, h E ( β ) i and C ( β ) are expressed as fol-lows: h E ( β ) i = P m mεe − βmε φ m P m e − βmε φ m (41) C ( β ) = 1 T ( h E ( β ) i − h E ( β ) i ) (42)while the mean square radius of gyration h R G i ( β ) is givenby [57, 66]: h R G i ( β ) = P m h R G i m e − βmε φ m P m e − βmε φ m (43)with h R G i m = L P LI,J =1 h ( ~R I − ~R J ) i m denoting theaverage of the gyration radius computed over states with m contacts.Finally, we wish to discuss how the problem of rareevents has been treated in this work. While there isenough evidence that the pivot algorithm is ergodic, it isalso true that certain states occur very rarely and some-times may require the computations of billions of trialconformations before being obtained. This is the case ofvery compact ( m > L ) or very swollen ( m ∼
0) confor-mations. For polymers of any length L , at least up to L = 300, the maximal length that has been studied withthe present approach, states with a number of contactshigher than L should be considered as rare events. Forshort polymer knots, samples with m = L or more areextremely rare. When L = 50 or 70, the lowest energystates that have been reached after a few millions of sam-ples are characterized respectively by 48 and 69 contacts.Short polymers which are heavily knotted are also diffi-cult to be found in conformations with a small numberof contacts. For a trefoil with L = 50 or L = 70 the statewith m = 0 should be considered as rare, while for 8 m < m ≤ L . For short and topologicallycomplex knots, m has additionally been bounded frombelow to be greater than 0. In the case of some knots,which are not considered here, the restriction m < m affect the calculations at lower and higher temperaturesat most by a relative error of a few percent. The dataabout the average gyration radius show for example thatthe difference between the mean square gyration radiiobtained by cutting or not the five lowest values of m isvery small. On the other side, without imposing boundson m , we have observed that the flattening procedureof the energy histogram in the Wang-Landau algorithmrequires an enormous amount of sampling, especially if0the rare events occur when the modification factor f ν issmall. B. Simulation results based on the knot invariant ̺ ( C ) In Fig. 2, left panel, the results for the specific en-ergy and the heat capacity in the attractive case are dis-played. Let us note that in all figures we have used thesymbol T for the normalized temperature Tε . As it is pos-sible to see, the specific energy increases with increasingtemperatures as it should be, because as the tempera-ture grows, more and more energetic states are excited.Moreover, at high temperatures longer polymers have ahigher specific energy than shorter ones. This behaviormay be explained by the following two observations. Thefirst observation, which is true for both short and longpolymers, is that, when the temperature is very high,the energy gap ε < V IJ defined inEq. (38) is much smaller than the energy related to thethermal fluctuations. As a consequence, the polymer issupposed to be in a more swollen conformation than atlow temperatures, where the interactions dominate overthe thermal fluctuations. The second observation is thatthe effects of knotting are less and less important withincreasing polymer lengths. To convince ourselves thatthis is the case, it is sufficient to imagine a knot which islocalized in a small part of the polymer. If the polymeris long, the localized part will be not relevant in compari-son with the rest of the polymer, which will behave moreor less like a unknotted ring. Also numerical simulationsconfirm this trend, see for example [48]. Moreover, ourcalculations show that the differences in the specific en-ergy and heat capacity between two trefoils of lengths L = 90 and L = 120 are much less marked than thesame differences between two trefoils of lengths L = 70and L = 90. Indeed, the data of the trefoil with L = 120have not been reported in Fig. 2(a) and (b) because itis difficult to distinguish them from those of the trefoilwith L = 90. Combining the two observations above,one can conclude that, at high temperatures, the poly-mer will tend to swell, but the effect of the topologicalconstraints will be to counteract this swelling process.On the other side, the topological constraints and theireffects become less important when polymers are long.As an upshot, if the temperature T and the knot topol-ogy are kept fixed, we expect that the average distancebetween the monomers and thus the specific energy of theknot will increase with increasing polymer length when T >>
1. This explains the behavior of the specific en-ergy of Fig. 2 at higher temperatures. The behavior of h E ( β ) i at lower temperatures should be taken with somecare because, as already mentioned, the number of con-tacts has been limited by the condition m ≤ L , so thatthe lowest energy state cannot be reached. Always forthat reason, the specific energy at zero temperature isslightly different for different lengths. Concerning the heat capacity, we can see in Fig. 2,right panel, that there is only one sharp peak in thewhole temperature region. The interpretation of thispeak is somewhat difficult. It is certainly a pseudo phasetransition caused by the fact that we are working on alattice with a finite size system, as those discussed inRefs. [67, 68]. Similar pseudo phase transitions have beenalready observed in knots, see [47]. In [68] it has beenstressed that, along with the advances in constructinghigh resolution equipment, such pseudo phase transitionin real system become more and more important. It islikely that the peaks in the heat capacity correspond tothe transition of the knot from a frozen crystallite state toan expanded state similar to what happens in the case ofa single polymer chain discussed in [59]. The data con-cerning the gyration radius in Fig. 6, left panel, seemsto confirm this hypothesis, because the gyration radiusstats to grow abruptly more or less in the same range oftemperatures in which the peak is appearing. This hy-pothesis could also explain the absence of a second peakor shoulder connected to the coil-globule transition. Asit was discussed in [59], in fact, if the range of the inter-actions is very short, then an open chain admits just twopossible 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 lowestenergy conformations exhibit some kind of ordering, butprobably a full ordering is forbidden by the topology ofthe knot. In conclusion, further analysis is necessary inorder to identify the nature of these transitions.Analogous considerations can also be made in the re-pulsive case displayed in Fig. 3. One difference is that lowtemperatures correspond now to the swollen state, whileat high temperatures, when the energy barrier ε becomesnegligible with respect to the energy carried by the ther-mal fluctuations, the monomers are allowed to get nearerand the average number of contacts is growing. As a con-sequence, the gyration radius for repulsive interactionsdecreases at higher temperatures as shown in Figs. 6(b)and 7(b). Moreover, while as expected the specific energyincreases with increasing temperatures, it turns out that,contrarily to what happens when the forces are attrac-tive, longer polymers have a lower specific energy thanshorter ones at any temperature. To explain this de-crease of the energy of longer polymers, we recall thatthe swelling of the polymer knot is hindered by the topo-logical constraints, in such a way that more complex knotconfigurations correspond to more compact polymer con-formations. As well, due to the fact that those parts ofthe polymer trajectory which are affected by the topolog-ical constraints will become less and less important whenthe polymer length increases, it is licit to conclude thatlonger polymers will in the average admit conformationsthat are more swollen with respect to those of shorterpolymers as we argued in the attractive case. The differ-ence is that, if the interactions are repulsive, more swollenconformations are less energetic, which explains why thespecific energy of longer polymers is in the average lower1 FIG. 2: In the attractive case, the specific energy (in units of ε ) and heat capacity of the trefoil as functions of the normalizedtemperature T = Tε . The polymer length can take the values L = 50 (circles), 70 (rectangles) and 90 (diamonds).FIG. 3: In the repulsive case, the specific energy (in units of ε ) and heat capacity of the trefoil as functions of the normalizedtemperature T = Tε . The polymer length can take the values L = 50 (circles), 70 (rectangles) and 90 (diamonds). than that of shorter polymers. The presence of sharppeaks in the heat capacities, see Figs. 3(b) and 5(b), hasnot a straightforward interpretation like those occurringwhen the interactions are attractive. Apparently, as ar-gued in [48], at T = 0 the knot is in one of the lowestenergy conformations which are allowed. As the temper-ature increases, at the beginning the system is unableto pass to a more compact configuration unless T ∼ ε .After that threshold, the knot more and more contactsare possible between the monomers unless saturation isreached.To check the effects of the topology on the behavior of the polymer, we have tested different knot configura-tions of the same length. We present here the data ofpolymers with L = 70. We observe that, in the attrac-tive case, the increasing of the knot complexity resultsin the decreasing of the specific energy and heat capac-ity at high temperatures, see Figs. 4(a) and (b). As onemay expect from the previous discussion, in the repulsivecase it is exactly the converse, i. e. the increasing of theknot complexity results in the increasing of the specificenergy and heat capacity when the temperature grows,as shown in Fig 5(a) and (b). Once again, this is due tothe fact that, if the topology of the knot is more com-2 FIG. 4: Specific energy (in units of ε ) and heat capacity for knots 5 , 4 and 3 in the attractive case. The polymers havelength L = 70. Inset: Specific energy and heat capacity for knots 5 , 4 and 3 with length L = 120 in the attractive case. plex, the knot conformation contains more contacts andis more compact. To have more contacts implies that, inthe average, the energy of each monomer and thus thespecific energy is lower if the forces are attractive andhigher if the forces are repulsive. This scenario is con-firmed by the data on the gyration radius of differentknots with length L = 70 in both the attractive and re-pulsive case, see Fig. 7. The gyration radii h R G i of thethree different knots 3 , and 5 satisfy the inequality( h R G i ) > ( h R G i ) > ( h R G i ) independently if the in-teractions are attractive or repulsive. Finally, as we pre-viously mentioned, at higher temperatures the effects ofthe interactions should become negligible. Accordingly,even if the gyration radius exhibits a different behaviordepending if the interactions are repulsive or attractive,in Fig. 8 we see that the values of the radius of gyra-tion computed in these two different cases get closer andcloser when the temperature is increasing. Eventually,they should converge if the temperature is high enough. C. Comparison with the PAEA method
The PAEA method allows to use topology changingpivot transformations in numerical simulations of poly-mer systems and prevents the accidental transition toanother topology without making use of topological in-variants. The idea of the PAEA method is based on thefollowing observation. After a pivot transformation isperformed on a randomly chosen element ∆ k of length N of a knot, ∆ k is transformed into the new element∆ k ′ . Since ∆ k and ∆ k ′ have their two ends in common,because these ends are untouched by the pivot trans-formation, it is easy to realize that, together, ∆ k and ∆ k ′ form a small loop or, if N is large, a set of smallloops of total length 2 N . The topology is preserved bychecking whether the part of the knot unaffected by thepivot transformation crosses an arbitrary surface spannedaround the small loop(s). If the crossing happens, thenthe trial conformation will be rejected and the knot un-dergoes another pivot transformation. Otherwise, thetrial conformation is accepted and a new pivot transfor-mation is applied to it. More details can be found in[48].What is crucial for the success of the PAEA methodis to classify all the possible small loops of 2 N segmentsthat can be produced after a random pivot transforma-tion and to construct suitable surfaces having these smallloops as borders. When the length N of the segmentschanged by the pivot transformations in greater thanfive, the number of all loops of this kind becomes large.With increasing values of N , it becomes more and moredifficult to construct the surfaces mentioned above. Un-fortunately, this is a not problem that can be solved au-tomatically by the computer using some algorithm. Forthis reason, up to now the PAEA algorithm has been de-veloped only in the case N = 4 ,
5. The advantage of thetechnique discussed in this work based on the knot invari-ant ̺ ( C ) is that the number of segments involved in thepivot transformations may arbitrarily range within theinterval 1 < N ≤ L . Larger values of N change biggerportions of the knot and this leads to a faster equilibra-tion process than the PAEA method. Moreover, we haveobserved that, as the length N of the segments to bechanged is increasing, the number of samples N samples needed to get a flat energy histogram is decreasing. Thisfact is shown in Table II, where as an example the ratio( N samples ) P AEA / ( N samples ) ̺ ( C ) is displayed in the case3 FIG. 5: Specific energy (in units of ε ) and heat capacity for knots 5 , 4 and 3 in the repulsive case. The polymers have length L = 70. Inset: Specific energy and heat capacity for knots 5 , 4 and 3 with length L = 120 in the repulsive case.FIG. 6: Mean square gyration radius of trefoil knots of lengths L = 50 , ,
90 as a function of the temperature. (a) Plot of themean square gyration radius in the attractive case; (b) Plot of the mean square gyration radius in the repulsive case. of the knot 3 with length L = 120. In the calcula-tions with the PAEA method the pivot transformationswere limited to N = 4, while in the calculations with the ̺ ( C ) invariant N = 24. All the approximation degrees ν = 0 , . . . ,
16 used in computing the density of states havebeen listed. As we can see from Table II, the number ofsamples ( N samples ) P AEA needed in the PAEA method tomake the energy histogram flat is always larger than thecorresponding number of samples ( N samples ) ̺ ( C ) neces-sary to distinguish the topology with the help of ̺ ( C ).The explanation of this increasing in the efficiency of the method based on the knot invariant ̺ ( C ) with respect tothe PAEA method is that with the invariant ̺ ( C ) largepivot transformations are allowed and they are able tomodify a relevant portion of the polymer. The larger isthe number N of segments affected by a pivot transfor-mation, the greater is the difference between the numbersof contacts of the knot before and after the transforma-tion. As a consequence, with a large pivot transformationit is possible to jump from conformations that have verydifferent values of m , thus accelerating considerably theexploration of the set of all possible conformations of a4 FIG. 7: Mean square gyration radius of knots 5 , 4 and 3 of length L = 70 as a function of the temperature. (a) Plot of themean square gyration radius in the attractive case; (b) Plot of the mean square gyration radius in the repulsive case.FIG. 8: Convergence at high temperatures of the mean squaregyration radii computed in the attractive and repulsive casesfor a trefoil knot of length L = 70. knot compatible with the given topological configuration.Of course, the calculation of ̺ ( C ) implies the evalua-tion of quadruple integrals with a good precision, other-wise there will be not a sufficient control of the topology.To perform the necessary integrations becomes challeng-ing with increasing polymer lengths even using a MonteCarlo integration method as we did here. However, sincewhat we need to calculate is the second Conway coeffi-cient and its value is dependent on the topological con-figuration of the knot, but not on its length, it is alwayspossible to shorten the knot provided its topology is notaffected. A simple way to do that is to group together TABLE II: The ratio ( N samples ) PAEA / ( N samples ) ̺ ( C ) (de-noted with P AEA/̺ ( C ) in this table) obtained from the com-putations of the density of states of the knot 3 with length L = 120. The initial value of the modification factor is f = e . f ν P AEA/̺ ( C ) f ν P AEA/̺ ( C ) f = e f = √ f f = √ f f = √ f f = √ f f = √ f f = √ f f = √ f f = √ f f = √ f f = √ f f = √ f f = √ f f = √ f f = √ f f = √ f f = √ f triplets of contiguous segments of the knot. This allowsto decrease by a factor three the length of the polymerwhich, after this procedure, will be of course no longerdefined on a simple cubic lattice. Such reduction methodcan be combined with other reduction algorithms, like forinstance the KMT reduction scheme proposed in [63, 69].Finally, we wish to come back to the problem of ergod-icity. We have seen that the PAEA method is restrictedup to now to small pivot transformations with N = 4 , ̺ ( C )in order to distinguish the topology. With this knot in-variant, it is in fact possible to consider pivot transfor-5 FIG. 9: Comparison of the densities of states φ m for a trefoilknot with length L = 120 computed using the PAEA method(grey squares) and the knot invariant ̺ ( C ) evaluated with aMonte Carlo (MC) algorithm (black circles). mations involving an arbitrary number of segments as ex-plained before. We have found that the densities of statescomputed with the PAEA method and with the invari-ant ̺ ( C ) are in complete agreement with each other. InFig. 9 it is shown for example the case of the densities ofstates for a trefoil knot of length L = 120 computed withthe PAEA method (grey squares) and with the invariant ̺ ( C ) evaluated using the Monte Carlo integration algo-rithm of Section III (black circles). As it is possible torealize, both densities of states coincide. VI. CONCLUSIONS
In this work we have computed the specific energy, thespecific heat capacity and the radius of gyration of thetrefoil knot 3 , the figure-eight 4 and the knot 5 . Poly-mer knots have been generated on a simple cubic latticeand the conformations used in the Wang-Landau algo-rithm to get the density of states Ω m have been sampledby means of pivot transformations. The topology of poly-mer knots needs to be preserved during the sampling pro-cedure. To this purpose, the knot invariant ̺ ( C ) has beenemployed. We have considered both attractive and repul-sive short-range interactions between monomers. First,we have analyzed the simplest possible knot 3 at variouslengths L = 50 , , , h R G i is displayed in Figs. 6(a) and 6(b) for the attrac-tive and repulsive cases respectively. It turns out thatlonger polymers have bigger mean square gyration ra-dius independently of the fact that the interactions arerepulsive or attractive. Indeed, as it is intuitive, longerpolymers should occupy larger volumes. Our simulationsshow also that, in the attractive case, h R G i increases withthe increasing of the temperature while, on the contrary,it decreases in the repulsive case. This phenomenon isconnected with the analogous increasing and decreasingof the specific energy mentioned before. As a matter offact, shorter radii of gyration imply higher specific ener-gies in the repulsive case and lower specific energies inthe attractive case.The influences of topology on the thermal propertiesof polymer knots have been studied by comparing knotsof different types but of the same length. In Figs. 4, 5and 7 are respectively reported the data of the specificenergy, heat capacity and gyration radius for the knots3 , and 5 with length L = 70. As it can be observedfrom Fig. 4, in the attractive case the specific energy andheat capacity decrease with the knot complexity. This isdue to the fact that with the increasing of the complex-ity of the topological configuration while the knot lengthis kept fixed, the conformation of the polymer becomesmore compact and thus the number of contacts betweenthe monomers becomes larger. The opposite situationis observed in the repulsive case. These results are inagreement with the conclusions of Refs. [48] and [57].Beside studying the thermodynamic properties of poly-6mer knots, the second purpose of this work has been acheck of the feasibility of the use in numerical simulationsof knot invariants which are given in the form of multiplecontour integrals. For this reason, we have consideredone of the simplest knot invariant, namely the quantity ̺ ( C ) which is related to the second coefficient of the Con-way polynomial. The most serious disadvantage of thiskind of invariants is that the evaluation of the multipleintegrals is time consuming even within the Monte Carloapproach adopted here. However, one should keep inmind that these knot invariants can be applied to anyspatial curve C representing the knot under investiga-tion, not necessarily defined on a cubic lattice or withsegments that are straight lines. This allows to reducethe number of segments composing the polymer consid-erably. For example, it is possible to replace in the knottrajectory up to three segments with a single segments,reducing in this way the polymer length by a factor three.Up to now, our method can be efficiently exploited withpolymer knots with lengths L ≤ ̺ ( C ) ofthe knot invariant ̺ ( C ) can be performed analytically.In this way, the time for calculating ̺ ( C ), which scaleswith the polymer length as t ∼ L , may be reduced to t ∼ L . While there are faster algorithms to preserve thetopology, like for instance the PAEA method discussedin [48], which is both exact and very fast for small pivottransformations involving up to N = 5 segments, the use of knot invariants in the form of multiple contour inte-grals has the advantage that it works with any numberof segments N . This allows a faster equilibration of thepolymer starting from a seed configuration and limitsthe number of sampling during the calculations with theWang-Landau Monte Carlo algorithm.In the future, we plan to extend the present approach,which is valid for knots, to the case of three linked poly-mers. This is possible because the triple Milnor linkinginvariant [55], which is able to distinguish the topologyof a link formed by three knots, is composed by four in-tegrals that have the same tensor structure of the twocomponents ̺ ( C ) and ̺ ( C ) of ̺ ( C ), see for example[74]. Appendix A
As shown in [54], in the contribution ̺ ( C ) to theknot invariant ̺ ( C ) the integration over the variable ~ω in Eq. (2) can be performed exactly giving as a result: ̺ ( C ) = − π I dx µ Z x dy ν Z y dz ρ H µ,ν,ρ ( ~y − ~x, ~z − ~x )(A1)After putting a = ~y − ~x and b = ~z − ~x , the tensor H µ,ν,ρ ( a, b ) may be expressed as follows: H µ,ν,ρ ( a, b ) = C C C [ δ νρ ( a − b ) µ + δ µρ b ν − δ µν a ρ ] − C C C ǫ λστ a σ b τ (cid:20) ǫ ρµα δ λν (cid:18) a α + b α | a || b | (cid:19) + ǫ νµα δ λρ (cid:18) b α + a α | b || a | (cid:19)(cid:21) + C C ǫ λστ a σ b τ (cid:26) ǫ ρµα δ λν (cid:20) b α | a − b | − | a | b + ( a − b ) α | a | + | b | ( a − b ) (cid:21) + ǫ νµα δ λρ (cid:20) a α | a − b | − | b | a + ( b − a ) α | a | + | b | ( a − b ) (cid:21)(cid:27) (A2)where C ( a, b ) = 2 π | a || b || a − b | (A3) C ( a, b ) = 1 | a || b | + a µ b µ (A4) C ( a, b ) = | a | + | b | − | a − b | (A5) | a | , | b | and | a − b | mean the module of the variable a, b and ( a − b ), like | a | = p ( y − x ) + ( y − x ) + ( y − x ) . a α , b α , ( a − b ) α mean the α -th component of the corre-sponding variable, α = 1 , , Acknowledgments
The support of the Polish National Center of Science,scientific project No. N N202 326240, is gratefully ac-knowledged. The simulations reported in this work wereperformed in part using the HPC cluster HAL9000 of theComputing Centre of the Faculty of Mathematics andPhysics at the University of Szczecin. [1] A. Yu. Grosberg,
Phys.-Usp. (1997), 12. [2] W. R. Taylor, Nature (London) (2000), 916. [3] V. Katritch, J. Bednar, D. Michoud, R. G. Scharein, J.Dubochet and A. Stasiak, Nature (1996), 142.[4] V. Katritch, W. K. Olson, P. Pieranski, J. Dubochet andA. Stasiak,
Nature (1997), 148.[5] M. A. Krasnow, A. Stasiak, S. J. Spengler, F. Dean, T.Koller and N. R. Cozzarelli,
Nature (1983), 559.[6] B. Laurie, V. Katritch, J. Dubochet and A. Stasiak,
Bio-phys. Jour. (1998), 2815.[7] J. I. Su lkowksa, P. Su lkowksa, P. Szymczak and M.Cieplak, Phys. Rev. Lett. (2008), 058106.[8] J. F. Marko,
Phys. Rev. E (2009), 051905.[9] Z. Liu, E. L. Zechiedrich and H. S. Chan, Biophys. J. (2006), 2344.[10] S. A. Wasserman and N. R. Cozzarelli, Science (1986), 951.[11] D. W. Sumners, Knot theory and DNA,
New ScientificApplications of Geometry and Topology, edited by D. W.Sumners, Proceedings of Symposia in Applied Mathemat-ics, American Mathematical Society, Providence, RI (1992), 39.[12] A. V. Vologodski, A. V. Lukashin, M. D. Frank-Kamenetski and V. V. Anshelevich, Zh. Eksp. Teor. Fiz. (1974), 2153; Sov. Phys. JETP (1975), 1059; M.D. Frank-Kamenetskii, A. V. Lukashin and A. V. Volo-godskii, Nature (London) (1975), 398.[13] E. Orlandini and S. G. Whittington,
Rev. Mod. Phys. (2007), 611.[14] C. Micheletti, D. Marenduzzo and E. Orlandini, Phys.Reports (2011), 1.[15] S. D. Levene, C. Donahue, T. C. Boles and N. R. Coz-zarelli,
Biophys. J. (1995), 1036.[16] T. Vettorel, A. Yu. Grosberg and K. Kremer, Phys. Biol. (2009), 025013.[17] P. Virnau, Y. Kantor and M. Kardar, J. Am. Chem. Soc. (43) (2005), 15102.[18] P. Piera´nski, S. Przyby l and A. Stasiak,
EPJ E (2)(2001), 123.[19] J. Yan, M. O. Magnasco and J. F. Marko, Nature (1999), 932.[20] J. Arsuaga, M. Vazquez, S. Trigueros, D. W. Sumnersand J. Roca,
PNAS (2002), 5373.[21] J. Arsuaga, M. Vazquez, P. McGuirk, S. Trigueros, D.W. Sumners and J. Roca, PNAS (2005), 9165.[22] R. Metzler, A. Hanke, P. G. Dommersnes, Y. Kantor andM. Kardar,
Phys. Rev. Lett. (2002), 188101.[23] P. Pieranski, S. Clausen, G. Helgesen and A. T. Skjeltorp, Phys. Rev. Lett. (1996), 1620.[24] Y. Diao, C. Ernst and E. J. Janse van Rensburg, IdealKnots, A. Stasiak, V. Katritch and L. H. Kauffman(Eds) , (World Scientific, Singapore, 1998), 52.[25] D. W. Sumners,
Notices of the Am. Math. Soc. (5)(1995), 528.[26] L. Faddeev and A. J. Niemi, Nature (1997), 58.[27] F. Ferrari and I. Lazzizzera,
Phys. Lett. B (1998),167.[28] F. Ferrari and I. Lazzizzera,
Nucl. Phys. B (3)(1999), 673.[29] F. Ferrari, H. Kleinert and I. Lazzizzera,
Eur. Phys. J
B18 (2000), 645.[30] A. L. Kholodenko and T. A. Vilgis,
Phys. Rep. (1998), 251.[31] F. Ferrari,
Annalen der Physik (Leipzig) (2002) 4, 255.[32] F. Ferrari, Topological field theories with non-semisimplegauge group of symmetry and engineering of topological invariants , chapter published in
Trends in Field The-ory Research , O. Kovras (Editor), Nova Science Publish-ers (2005), ISBN:1-59454-123-X. See also the reprint ofthis article in
Current Topics in Quantum Field TheoryResearch , O. Kovras (Editor), Nova Science Publishers(2006), ISBN: 1-60021-283-2.[33] H. Kleinert,
Path Integrals in Quantum Mechanics,Statistics, Polymer Physics and Financial Markets ,(World Scientific Publishing, Singapore, 2009).[34] A. Yu. Grosberg, A. Feigel and Y. Rabin,
Phys. Rev. E (1996), 6618.[35] B. Marcone, E. Orlandini, A. L. Stella and F. Zonca, J.Phys. A (2005), L15.[36] S. R. Quake, Phys. Rev. Lett. (1994), 3317.[37] M. L. Mansfield and J. F. Douglas, J. Chem. Phys. (2010), 044903.[38] J. Suzuki, A. Takano and Y. Matsushita,
J. Chem. Phys. (2013), 024902.[39] J. des Cloizeaux and M. L. Mehta,
J. Physique (1979),665.[40] G. Brinke and G. Hadziioannou, Macromol. (1987),480.[41] J. P. J. Michels and F. W. Wiegel, Phys. Lett. A (1984), 381; Proc. R. Soc. A (1986), 269.[42] M. Baiesi and E. Orlandini,
Phys. Rev. E (2012),031805.[43] D. Meluzzi, D. E. Smith and G. Arya, Annual Review ofBiophysics , (2010), 349.[44] A. Yao, H. Marsuda, H. Tsukahara, M. K. Shimamura,and T. Deguchi, J. Phys. A (2001), 7563.[45] N. Madras, A. Orlitsky and L. A. Sepp, J. Stat. Phys. (1990), 159.[46] E. J. J. van Rensburg and S. G. Whittington, J. Phys.A (1990), 3573.[47] A. Swetnam, C. Brett and M. P. Allen, Phys. Rev. E (2012), 031804.[48] Y. Zhao and F. Ferrari, JSTAT J. Stat. Mech. (2012),P11022, arXiv 1209.1717.[49] A. Baumg¨artner and K. Binder,
J. Chem. Phys. (1979), 2541.[50] R. Everaers and K. Kremer, Phys. Rev. E , R37 (1996),1.[51] R. Everaers and K. Kremer, Entanglement effects inmodel polymer networks, Lecture Notes in Physics Vol-ume (1999), 221.[52] W. Michalke, M. Lang, S. Kreitmeier and D. Goritz,
Phys. Rev. E (2001), 012801.[53] J. H. Conway, An enumeration of knots and links, andsome of their algebraic properties , Computational Prob-lems in Abstract Algebra (Proc. Conf., Oxford, 1967),(Pergamon, Oxford, 1970), 329.[54] E. Guadagnini, M. Martellini and M. Mintchev,
Nucl.Phys. B (1990), 581.[55] J. Milnor,
Ann. of Math. (1954),177.[56] M. Kontsevich, Adv. Sov. Math. part 2 (1993), 137.[57] Y. Zhao and F. Ferrari, Acta Physica Polonica B (2013), 1001.[58] F. Wang and D. P. Landau, Phys. Rev. Lett. (2001),2050.[59] M. P. Taylor, W. Paul and K. Binder, J. Chem. Phys. (2009), 114907.[60] J. M. F. Labastida and A.V. Ramallo,
Phys. Lett. B (1990), 92.[61] A. S. Cattaneo, P. Cotta-Ramusino and M. Martellini, Nucl. Phys. B (1-2) (1995), 355.[62] T. Deguchi and K. Tsurusaki,
Numerical application ofknot invariants and universality of the random knotting ,published in
Knot Theory , Banach Center Publications (1998), 77.[63] K. Koniaris and M. Muthukumar, Phys Rev Lett (1991), 2211.[64] E. Witten, Commun. Math. Phys. (1989), 351.[65] C. Zhou and R. N. Bhatt,
Phys. Rev
E72 (2005), 025701(R).[66] P. N. Vorontsov-Velyaminov, N. A. Volkov, A. A.Yurchenko and A. P. Lyubartsev,
Polymer Science, Ser.A (2010), 742.[67] T. W¨ust and D. P. Landau, Phys. Rev. Lett. (2009),178101.[68] T. Vogel, M. Bachmann and W. Janke,
Phys. Rev.
E76 (2007), 061803. [69] W. R. Taylor and A. Asz´odi,
Protein Geometry, Classi-fication, Topology and Symmetry , (New York, Taylor &Francis, 2005).[70] F. Rampf, W. Paul and K. Binder,
Europhys. Lett. (2005), 628634.[71] H. B. Kolli and K. P. N. Murthy, Phase Transition ina Bond Fluctuating Linear Homo Polymer , AIP Conf.Proc. (2010), 117.[72] Z. Wang and X. He,
J. Chem. Phys. (2011), 094902.[73] L. I. Klushin and A. M. Skvortsov,
J. Phys. A: Math.Theor. (2011), 473001.[74] L. Leal and J. Pineda, Mod. Phys. Lett. A23