text:
# R-GMA: First results after deployment

## I Introduction

The Grid Monitoring Architecture (GMA)ref:perf-arch  of the GGF, as shown in Figure 1, consists of three components: Consumers, Producers and a directory service, which we prefer to call a Registry).

In the GMA Producers register themselves with the Registry and describe the type and structure of information they want to make available to the Grid. Consumers can query the Registry to find out what type of information is available and locate Producers that provide such information. Once this information is known the Consumer can contact the Producer directly to obtain the relevant data. By specifying the Consumer/Producer protocol and the interfaces to the Registry one can build inter-operable services. The Registry communication is shown on Figure 1 by a dotted line and the main flow of data by a solid line.

The current GMA definition also describes the registration of Consumers, so that a Producer can find a Consumer. The main reason to register the existence of Consumers is so that the Registry can notify them about changes in the set of Producers that interests them.

The GMA architecture was of course devised for monitoring but we think it makes an excellent basis for a *combined* information and monitoring system. We have argued beforeref:time-bncod-2001  that the only thing which characterises monitoring information is a time stamp, so we insist upon a time stamp on all measurements - saying that this is the time when the measurement was made, or equivalently the time when the statement represented by the tuple was true.

The GMA does not constrain any of the protocols nor the underlying data model, so we were free when producing our implementation to adopt a data model which would allow the formulation of powerful queries over the data.

R-GMA is a relational implementation of the GMA which brings the power and flexibility of that model. R-GMA creates the impression that you have one RDBMS per VO. However it is important to appreciate that what our system provides is a way of using the relational model in a Grid environment and that we have *not* produced a general distributed RDBMS. All the producers of information are quite independent. It is relational in the sense that Producers announce what they have to publish via an SQL CREATE TABLE statement and publish with an SQL INSERT and that Consumers use an SQL SELECT to collect the information they need.

R-GMA is built using servlet technology and is being migrated rapidly to web services and specifically to fit into an OGSAref:ogsa-spec  framework.

## II Query types and Producer Types

We have so far defined not just a single Producer but five different types: a DataBaseProducer, a StreamProducer, a ResilientProducer, a LatestProducer and a CanonicalProducer. All appear to be Producers as seen by a Consumer - but they have different characteristics. The CanonicalProducer, though in some respects the most general, is somewhat different as there is no user interface to publish data via an SQL INSERT statement. Instead it triggers user code to answer an SQL query. The other Producers are all Insertable; this means that they all have an interface accepting an SQL INSERT statement.

The other producers are instantiated and given the description of the information they have to offer by an SQL CREATE TABLE statement and a WHERE clause expressing a predicate that is true for the table. Currently this is of the form WHERE (column\_1=value\_1 AND column\_2=value\_2 AND ...). To publish data, a method is invoked which takes the form of a normal SQL INSERT statement.

Three kinds of query are supported: History, Latest and Continuous. The history query might be seen as the more traditional one, where you want to make a query over some time period - including.“all time”. The latest query is used to find the current value of something and a continuous query provides the client with all results matching the query as they are published. A continuous query is therefore acting as a filter on published data.

The DataBaseProducer supports history queries. It writes each record to an RDBMS. This is slow (compared to a StreamProducer) but it can handle joins. The StreamProducer supports continuous queries and writes information to a memory structure where it can be picked up by a Consumer. The ResilientStreamProducer is similar to the StreamProducer but information is backed up to disk so that no information is lost in the event of a system crash. The LatestProducer supports latest queries by holding only the latest records in an RDBMS.

Latest records are defined in terms of something similar to a primary key. Each record has a time stamp, one or more fields which define what is being measured (e.g. a hostname) and one or more fields which are the measurement (e.g. the 1 minute CPU load average). The time stamp and the defining fields are close to being a primary key - but as there is no way of knowing who is publishing what across the Grid, the concept of primary key (as something globally unique) makes no sense. The LatestProducer will replace an earlier record having the same defining fields, as long as the time stamp on the new record is more recent, or the same as the old one.

Producers, especially those using an RDBMS, may need cleaning from time to time. We provide a mechanism to specify those records of a table to delete by means of a user specified SQL WHERE clause which is executed at intervals which are also specified by the user. For example it might delete records more than a week old from some table or it may only hold the newest one hundred rows, or it might just keep one record from each day.

Another valuable component is the Archiver which is a combined Consumer-Producer. You just have to tell it what to collect and it does so on your behalf. An Archiver works by taking over control of an existing Producer and instantiating a Consumer for each table it is asked to archive. This Consumer then connects via the mediator to all suitable Producers and data starts streaming from those Producers, through the Archiver and into the new Producer. The inputs to an Archiver are always streams from a StreamProducer or a ResilientStreamProducer. It will re-publish to any kind of “Insertable”. This allows useful topologies of components to be constructed such as the one shown in Figure 2.

This shows a number of StreamProducers (labelled SP) which is normally the entry point to R-GMA. There is then a layer of Archivers (A) publishing to another StreamProducer. Finally there is an Archiver to a LatestProducer (LP) and an Archiver to a DataBaseProducer (DP) to answer both Latest and History queries.

We intend to allow some kinds of producer to answer more than one kind of query - but for now we are keeping it simple.

## III Applications of R-GMA

R-GMA has applications right across the Grid. First it can be used as a replacement for MDS. A small tool (GIN) has been written to invoke the MDS-like EDG info-providers and publish the information via R-GMA. The info-provider is a small script which can be invoked to produce information in LDIF format. All our information providers conform to the GLUE (http://www.cnaf.infn.it/ sergio/datatag/glue/) schema. Another tool (GOUT) is available to republish R-GMA data to an LDAP server for the benefit of legacy applications. However we expect that most applications will wish to benefit from the power of relational queries. GOUT is an Archiver with a Consumer which periodically publishes to an LDAP database. The GIN-GOUT combination is not efficient - but it works. Both GIN and GOUT are driven by configuration files which define the mapping between the LDAP schema and the relational schema.

R-GMA is also being used for network monitoring where the flexibility of the relational model offers a more natural description of the problem.

It is also being used to locate replica catalogs, and to publish information to two tables, the Service table and the ServiceStatus table. A service publishes its existence when it starts up with an entry in the Service table. It does this using a StreamProducer. An Archiver to a LatestProducer is instantiated to collect all Service information in one place. There are also processes which check the functioning of a service and publish the status frequently to the ServiceStatus table. This is also published via a StreamProducer and is collected by the same Archiver that is used for archiving the Service table. So the Service table says what should exist and the ServiceStatus gives the current state Grid wide.

GRM was written for monitoring parallel applicationsref:sztaki-cluster-monitoring-2000  where it writes logging information to a local file. This has recently been modified to make use of R-GMA for transport.

In addition CMS have adapted their BOSS system which previously wrote to a well known RDBMS to simply publish the job status information via R-GMA. This BOSS work is reported at this conference.

## IV Tools

There are a number of tools available to query R-GMA Producers. There is a command line tool, a Java application: Pulse, and the R-GMA Browser, which is accessible from a Web browser without any R-GMA installation. The Browser offers a few custom queries, and makes it easy for you to write your own.

The command line tool, which is written in Python, is the most powerful. It is designed to do simple things very easily - but if you want to do more complex things you must code them yourself. It supports one instance of each kind of producer and one archiver at any one time. You can also find what tables exist, find details of a table and issue any kind of query.

## V The registry and the mediator

The registry stores information about all producers currently available. Currently there is only one physical Registry per VO. This bottleneck and single point of failure is being eliminated. Code is being written to allow multiple copies of the registry to be maintained. Each one acts as master of the information which was originally stored in that Registry instance and has copies of the information from other Registry instances. Synchronisation is carried out frequently. Currently VOs are disjoint, we plan to allow information to be published to a set of VOs.

The mediator (which is hidden behind the Consumer interface) is the component which makes R-GMA easy to use. Producers are associated with views on a virtual data base. Currently views have the form:

> SELECT * FROM $`<`$table$`>`$ WHERE $`<`$predicate$`>`$

This view definition is stored in the Registry. When queries are posed, the Mediator uses the Registry to find the right Producers and then combines information from them.

## VI Architecture

R-GMA is currently based on Servlet technology. Each component has the bulk of its implementation in a Servlet. Multiple APIs in Java, C++, C, Python and Perl are available to communicate with the servlets. The basic ones are the Java and C++ APIs which are completely written by hand. The C API calls the C++ and the Python and Perl are generated by SWIG. We make use of the Tomcat Servlet container. Most of the code is written in Java and is therefore highly portable. The only dependency on other EDG software components is in the security area.

Figure 3 shows the communication between the APIs and the Servlets. When a Producer is created its registration details are sent via the Producer Servlet to the Registry (Figure 3a). The Registry records details about the Producer, which include the description and view of the data published, *but not the data itself*. The description of the data is actually stored as a reference to a table in the Schema. In practise the Schema is co-located with the Registry. Then when the Producer publishes data, the data are transferred to a local Producer Servlet (Figure 3b).

When a Consumer is created its registration details are also sent to the Registry although this time via a Consumer Servlet (Figure 3c). The Registry records details about the type of data that the Consumer is interested in. The Registry then returns a list of Producers back to the Consumer Servlet that match the Consumers selection criteria.

The Consumer Servlet then contacts the relevant Producer Servlets to initiate transfer of data from the Producer Servlets to the Consumer Servlet as shown in Figures 3d-e.

The data are then available to the Consumer on the Consumer Servlet, which should be close in terms of the network to the Consumer (Figure 3f).

As details of the Consumers and their selection criteria are stored in the Registry, the Consumer Servlets are automatically notified when new Producers are registered that meet their selection criteria.

The system makes use of soft state registration to make it robust. Producers and Consumers both commit to communicate with their servlet within a certain time. A time stamp is stored in the Registry, and if nothing is heard by that time, the Producer or Consumer is unregistered. The Producer and Consumer servlets keep track of the last time they heard from their client, and ensure that the Registry time stamp is updated in good time.

## VII Results so far

Unfortunately we have few results to offer at this stage. It has taken some time to get from the state of having something which passes all its unit tests (about 400 for the Java API) to a stable distributed system - which we think we now have. We have just started running performance tests to understand the behaviour of the code. We have so far tested with many Producers, and one Archiver feeding into a LatestProducer which is then queried to make sure that the Archiver is keeping up with the total flow of data. The Producers are publishing data following the pattern expected of a “typical” site having an SE (Storage Element) and 3 CEs (Computing Elements). We found that we were able to support around 40 to 50 sites publishing data every 30 seconds.

A producer is able to publish individual tuples or a vector of tuples. Changing this buffering had little effect on the maximum number of sites the system is able to support. We will use tools to analyse the code to understand precisely where the bottlenecks are occurring - in particular to find out why buffering did not have the beneficial effect we expected.

In general we expect to be able to do better by modifying tomcat settings, the virtual machine settings and getting more physical memory. To achieve better performance we may need a layer of Archivers combining streams into bigger streams so as to limit the fan-in to any one node. The other way to obtain significantly better performance is not to attempt to get all the information into one place. As the mediator becomes more powerful, it will be able to make use of multiple LatestProducer archives, and carry out a distributed query over them.

The effort involved in making meaningful measurements on such a system as R-GMA should not be underestimated!

## VIII Conclusion

We have a useful architecture and an effective implementation with a number of components which work well together. We hope that R-GMA will have a long, happy and useful life, both in its current form and when reincarnated within an OGSA framework. For more details of R-GMA, please see: http://hepunx.rl.ac.uk/edg/wp3/

###### Acknowledgements.

The authors wish to thank our patient users, the EU and our national funding agencies.


format: text

text:
# 1 Introduction

## 1 Introduction

The attempts to unify fundamental interactions have produced a lot of various models (see, for example,  and references therein). All of them are of rich mathematical structure and most of them are to some degree based on the Riemannian geometry. That’s why we believe that the quantization of the latter is rather important. Even if it has nothing to do with the real gravity, it may play an important role in the further construction of the realistic unified theory.

Recently it has been paid much attention to the dynamical triangulation (DT) approach to quantization of Riemannian geometry  in $`2`$, $`3`$ and $`4`$ dimensions111$`5`$ \- dimensional Euclidean quantum gravity was considered in .. For $`D=2`$ the DT model has a well - defined continuum limit consistent with the predictions of the continuum theory . At $`D=3,4`$ the Euclidean DT models with the discretized Einstein action have two phases: the crumpled phase with infinite fractal dimension and the elongated one, which resembles branched polymer model with the fractal dimension close to $`2`$. It appears that the introducing of the causal structure to these models (correspondent to the transition from Euclidean to Lorentzian quantum gravity) or coupling them to matter changes their behavior and make them more realistic . Nevertheless, the pure Euclidean gravity is still of interest as an area of developing the investigation methods. Moreover, mechanisms observed within this model can be related to more realistic models.

The mathematical structures related to the unification of fundamental forces may include the concept of high dimensional ($`D>4`$) space. In particular, all known superstring models become unambiguous only in the $`10`$ \- dimensional space - time. Unfortunately, most of information about their structure comes from the perturbative methods and the (incomplete) investigation of certain special excitations (such as $`D`$ \- branes). It seems that the nonperturbative investigation of those models would become exhaustive only if the numerical lattice methods are used. The modern lattice theory has an experience of dealing with quantum gravity models in lower dimensions. In particular, the dynamical triangulation method was applied to the systems of $`D=2,3,4`$ dimensions. Therefore we guess that it is reasonable to apply this method to the higher dimensional problems. The investigation of the pure Einstein Euclidean gravity can become the first step. The work related to this step is the content of the present paper.

Namely, we consider the ten - dimensional Euclidean dynamical triangulations of spherical topology. Their behavior was expected to be similar to that of $`3,4`$ and $`5`$ \- dimensional systems. This supposition was confirmed partially. However, it turns out that there are a few features that are not present in the lower dimensional models. For example, the $`10D`$ model does have two phases that are similar to the crumpled and elongated phases of $`3`$, $`4`$ and $`5`$ \- dimensional models. But in our case the phase transition (at least for the observed volumes $`V=8000`$ and $`V=32000`$) is at negative gravitational coupling constant. This means that at ”physical” positive gravitational constant $`G`$ the considered $`10`$ \- dimensional system cannot exist in the crumpled phase unlike $`3`$, $`4`$ and $`5`$ \- dimensional systems. Also this means that for the observed lattice sizes the natural state of the $`10`$ D system (when the Einstein - Hilbert action is turned off) resembles branched polymer while in the low dimensional systems the natural state belongs to the crumpled phase.

It is worth mentioning that the considered volumes are not consistent with the system of ”physical” dimension $`D=10`$. It follows from the observation that already at linear size $`3`$ the analogous rectangular $`10D`$ lattice has the volume $`V=3^{10}=59049`$. However, the observed effective (fractal) dimension of the triangulated space, say, at $`V=32000`$ is $`4`$ in the crumpled phase and $`2`$ in the elongated one. So, the linear size of the system is expected to be of the order of $`(32000)^{\frac{1}{4}}10`$ and $`(32000)^{\frac{1}{2}}100`$ respectively. This is confirmed by the direct measurement of the linear extent.

## 2 The model.

In this section we shall remind briefly the definition of the considered model and the description of the numerical algorithm. For the complete review of the dynamical triangulations see  and references therein. For the full description of the numerical algorithm for arbitrary dimension see .

In the dynamical triangulation approach the Riemannian manifold (of Euclidian signature) is approximated by the simplicial complex obtained by gluing the $`D`$ \- dimensional simplices. Each simplex has $`D+1`$ vertices. All links are assumed to have the same length $`a`$. The metrics inside each simplex is supposed to be flat. Therefore the deviation from the flatness is concentrated on the boundaries of the simplices. The scalar curvature $`R`$ is zero everywhere except the bones ($`D2`$ \- dimensional subsimplices of the triangulation simplices). The Einstein - Hilbert action can be expressed through the number of bones and the number of simplices of the given triangulated manyfold:

$`S`$ $`=`$ $`{\displaystyle \frac{1}{16\pi G}}{\displaystyle R(x)\sqrt{g}d^Dx}={\displaystyle \frac{\mathrm{Vol}_{D2}}{16\pi G}}{\displaystyle \underset{\mathrm{bones}}{}}(2\pi O(\mathrm{bone})cos^1({\displaystyle \frac{1}{D}}))`$ (1)
$`=`$ $`{\displaystyle \frac{\mathrm{Vol}_{D2}}{8G}}(N_{\mathrm{bones}}{\displaystyle \frac{D(D+1)}{4\pi }}N_{\mathrm{simplices}}cos^1({\displaystyle \frac{1}{D}})),`$
where $`O(\mathrm{bone})`$ is the number of simplices sharing the given bone, $`\mathrm{Vol}_j=\frac{a^j\sqrt{j+1}}{j!\sqrt{2^j}}`$ is the volume of a $`j`$ \- dimensional simplex with the edges of length $`a`$, $`N_{\mathrm{bones}}`$ is the total number of bones and $`N_{\mathrm{simplices}}`$ is the total number of simplices.

The metric of the triangulated manyfold is completely defined by the way the simplices are glued together. Therefore, the functional integral over $`Dg`$ in this approach is changed by the summation over the different triangulations (we restrict ourselves with the spherical topology only):

$$Dg\underset{T}{}\frac{1}{C_T}$$
(2)
Where the sum is over the triangulations $`T`$ that approximate different Riemannian manyfolds222Two formally different triangulations may approximate the same Riemannian manyfold. Therefore the correspondent factorization should be implemented. and $`C_T`$ is the symmetry factor of the triangulation itself (the order of its automorphism group).

We consider the model, in which the fluctuations of the global invariant $`D`$ \- volume are suppressed. The partition function has the form:

$$Z_V=\underset{T}{}\frac{1}{C_T}\mathrm{exp}(S(T))=\underset{T}{}\frac{1}{C_T}\mathrm{exp}(\kappa _{D2}N_{D2}\kappa _DN_D\gamma (N_DV)^2)$$
(3)
where we denoted $`N_D=N_{\mathrm{simplices}}`$ , $`N_{D2}=N_{\mathrm{bones}}`$ , and $`\kappa _{D2}=\frac{\mathrm{Vol}_{D2}}{8G}`$. Unfortunately, it is not possible to construct an algorithm that generates the sum over the triangulations of the same volume. So the constant $`\gamma `$ is kept finite. $`\kappa _D(V,\kappa _{D2})`$ is chosen in such a way that

$$<N_D>=\underset{T}{}\frac{1}{C_T}\mathrm{exp}(\kappa _{D2}N_{D2}\kappa _DN_D\gamma (N_DV)^2)N_D=V$$
(4)
This provides that the volume fluctuates around the required value $`V`$. The fluctuations are of the order of $`\delta V\frac{1}{\sqrt{\gamma }}`$. In order to approach the model with constant $`N_D`$ we must keep $`\frac{\delta V}{V}<<1`$. In practice we use $`\gamma =0.005`$ for $`V=16000,32000`$. So $`\frac{\delta V}{V}10^3`$.

For the numerical investigations we used Metropolis algorithm in its form described in . It is based on the following Markov chain. Each step of the chain is the proposal of a deformation $`T_iT_f`$ of the given triangulation $`T_i`$, which is accepted or rejected with the probability $`𝒫(T_iT_f)`$ that satisfies the detailed balance condition

$$\mathrm{exp}(S(T_i))𝒫(T_iT_f)=\mathrm{exp}(S(T_f))𝒫(T_fT_i)$$
(5)
The definition of the proposed deformations is based on the following idea. Let us consider some closed $`D`$ \- dimensional simplicial manyfold of the topology of a $`D`$ \- dimensional sphere. Then, if a connected piece of our original triangulation is equal to a piece of this manyfold, we can replace it by the remaining part of the given manyfold. Thus we obtain the deformed triangulation with the same topology as the original one. Further, let us choose the boundary $`s_{D+1}`$ of a $`D+1`$ \- dimensional simplex $`s_{D+1}`$ as the mentioned above manyfold. There are $`D+1`$ opportunities to distinguish a piece of $`s_{D+1}`$ correspondent to $`p`$ \- dimensional subsimplices of $`s_{D+1}`$ ($`p=0,\mathrm{},D`$). The resulting deformation is called $`(p,Dp)`$ move .

It has been shown that via such moves it is possible starting from any triangulation to reach a triangulation that is combinatorially equivalent 333Two triangulations are combinatorially equivalent if they have subdivisions that are equivalent to each other up to the relabeling of (sub)simplices. The subdivision $`T_s`$ of the given triangulation $`T`$ is another triangulation such that the set of its vertices contains all vertices of $`T`$ and any simplex $`sT_s`$ either belongs to $`T`$ or belongs to some simplex of $`T`$. to an arbitrarily chosen other triangulation . This property is called ergodicity. Due to this property starting from the triangulation that approximates  the given Riemannian manyfold it is possible to reach a triangulation that approximates almost any other Riemannian manyfold. The exceptional cases, in which this is impossible, are commonly believed not to affect physical results.

In practise we choose randomly the type of the move ($`p\{0,\mathrm{},D\}`$), the simplex of a triangulation and its $`p`$ \- dimensional subsimplex. After that we check is there a vicinity of this $`p`$ \- dimensional simplex that is equivalent to the required piece of $`s_{D+1}`$. If so, the suggested move (and the correspondent subsimplex) is called legal and we proceed with checking the possibility to perform the move 444The move is not geometrically allowed if the resulting new simplex already exists in the given triangulation.. If the move is allowed we accept or reject it with the following probability:

$$p(T_iT_f)=\frac{1}{1+(1+\frac{N_D(T_f)N_D(T_i)}{N_D(T_i)})\mathrm{exp}(S(T_f)S(T_i))}$$
(6)
In our calculations we start from the triangulation of minimal size that is $`s_{D+1}`$. Then we allow it to grow randomly until the volume reaches the given value $`V`$. After that we proceed with the normal Metropolis process. During this process the coupling $`\kappa _D`$ is self - tuned automatically in order to get its required value satisfying (4). This is done via redefining it as $`\kappa \kappa +2\gamma (<N_D>V)`$ after each $`10`$ sweeps 555One sweep is $`V`$ suggestions of legal moves..

In order to make a check of our results, we have made two independent program codes. The main program was written in $`C^{++}`$ using modern methods of object \- oriented programming. The algorithm is partially based on the one described in  and on the ideas suggested in . The second program was written in Fortran $`77`$ and was used for checking the main program at small lattice sizes. The calculations were made within the parallel programming environment using the computation facilities of Joint Supercomputer Center (Moscow).

## 3 Numerical results.

We investigate the behavior of $`10D`$ model for $`V=8000`$ and $`V=32000`$. We considered the values of $`\kappa _D`$ varying from $`0.1`$ to $`0.1`$, where the system is found to exist in two phases. The phase transition point is at $`\kappa _D0.03`$ at $`V=8000`$ and at $`\kappa _D0.01`$ at $`V=32000`$.

The self - tuned value of $`\kappa _{10}`$ appears to be independent of $`V`$. In accordance with large volume asymptotics obtained in  its dependence upon $`\kappa _8`$ is linear with a good accuracy. The best fit is:

$$\kappa _{10}=15.57(1)\kappa _8+0.42(1)$$
(7)
We investigate the following variables, which reflect the properties of the triangulated manyfold.

1. The mean curvature carried by a simplex. The normalization is chosen in such a way that it is defined as

$$R=\frac{4\pi N_{\mathrm{bones}}}{D(D+1)N_{\mathrm{simplices}}cos^1(\frac{1}{D})}1$$
(8)
The results for $`V=8000`$ and $`V=32000`$ coincide with a good accuracy. They are represented in the figure 1. In contrast with the four dimensional case (at the observed values of $`\kappa _8`$) $`R>0`$ (in the four dimensional case it becomes negative at the value of $`\kappa _{D2}`$ close to $`0`$). The best fit to the curve of the figure 1 is:

$$R=0.218(1)+0.654(2)\kappa _81.97(3)\kappa _8^2$$
(9)
2. The geodesic distance between two simplices is the length of the shortest path that connects them. The points of this path correspond to the simplices. The links of the path correspond to the pairs of neighboring simplices. We denote the geodesic distance between the simplices $`u`$ and $`v`$ by $`\rho (u,v)`$. Let us fix a simplex $`s`$. Then the ball $`B_R(s)`$ of radius $`RZ`$ is consisted of the simplices $`u`$ such that $`\rho (u,s)`$ is less than or equal to $`R`$. The volume of the ball is defined as the number of simplices contained in it.

One of the most informative characteristics of the triangulated manyfold is the average volume $`𝒱(R)`$ of the balls (of the constant radius $`R`$) as a function of this radius. In practise during each measurement we calculate $`V_R(s)`$ for arbitrary chosen simplex $`s`$. Then we perform its averaging over the measurements (separately for each $`R`$). $`𝒱(R)`$ becomes constant at some value of $`R`$. This value is the averaged largest distance between two simplices of the manyfold, which is also called the diameter $`d`$. We found that for $`V=8000`$ at $`1<R<\frac{d}{2}`$ the dependence of $`\mathrm{log}𝒱`$ on $`\mathrm{log}R`$ is linear. For $`V=32000`$ the same takes place at $`4<R<\frac{d}{2}`$. The slopes give us the definition of the fractal (Hausdorff) dimension of the manyfold: $`\mathrm{log}𝒱=\mathrm{const}+𝒟\mathrm{log}R`$.

3. The mean fractal dimension $`𝒟(V)`$ of the manyfold for the observed volumes $`V`$ as a function of $`\kappa _{D2}`$ is represented on the Fig. 2. The figure indicates that there is a phase transition at critical $`\kappa _{D2}=\kappa _c(V)`$. For the observed volumes $`\kappa _c`$ appears to be small and negative ($`\kappa _c(8000)0.3`$ and $`\kappa _c(32000)0.1`$). For $`\kappa >\kappa _c`$ the fractal dimension is close to two in accordance with the supposition that similar to the $`4`$ \- dimensional case the correspondent phase resembles branched polymers. For $`\kappa <\kappa _c`$ $`𝒟(8000)3.2`$ and $`𝒟(32000)4`$. This is in accordance with the expectations that this phase has a singular nature and corresponds to $`𝒟(\mathrm{})=\mathrm{}`$. Further we shall call these two phases crumpled and elongated respectively. Our results give us arguments in favor of view that the nature of these phases is similar to that of the phases of the lower dimensional models. We must notice, however, that the complete proof has not yet been obtained.

4. The linear size of the system can be evaluated as $`V^{\frac{1}{\mathrm{D}}}`$. So, we expect that in the crumpled phase the diameter $`d10`$ while in the elongated phase $`d100`$. These expectations are in accordance with the direct measurements of the diameter as well as another parameter called linear extent (see ). The linear extent $``$ is defined as the average distance between two simplices of the triangulation:

$$=\frac{1}{V^2}\underset{u,v}{}<\rho (u,v)>$$
(10)
Due to its construction $``$ should be close to half a diameter. In our measurements we calculate the linear extent using the slightly different definition (that, anyway, should lead to the same result after averaging over the measurements). Namely, we calculate

$$(s)=\frac{1}{V}\underset{v}{}<\rho (s,v)>=\frac{1}{V}\underset{R}{}R(V_R(s)V_{R1}(s)).$$
(11)
Performing the averaging over the measurements we obtain the required value of the mean linear extent. Our results on $``$ are represented in the figure Fig. 3. One can see that the linear size of the manyfold is increasing very fast in the elongated phase, while in the crumpled phase it remains almost constant (and close to that of $`3`$, $`4`$ and $`5`$ dimensional models ). We found that the fluctuations of the linear extent almost absent in the crumpled phase and are of the order of the mean size of the manyfold in the elongated phase.

## 4 Discussion.

The force between two particles (of masses $`m`$) in the classical mechanics is equal to

$$F=\frac{Gm^2}{R^2},$$
(12)
where $`R`$ is the distance between them. Therefore negative gravitation coupling $`G`$ and, correspondingly, negative $`\kappa _2`$ in the four - dimensional model would correspond to repulsion instead of attraction. The same picture is, of course, valid for the negative $`\kappa _{D2}`$ in higher - dimensional models. This means that the description of gravity in the dynamical triangulation approach may appear only for the positive $`\kappa _{D2}`$. The phase transition in the pure Euclidean four - dimensional gravity is at $`\kappa _21`$ and the model exists at physical couplings both in the crumpled and in the elongated phases. In the ten dimensional model for the observed volumes the phase transition point is $`0.01`$. So at physical gravitational couplings the system cannot exist in the crumpled phase.

The fractal dimension at $`G>0`$ is found to be close to $`2`$. This means that in order to construct a dynamical triangulation model that corresponds to the dynamics of Riemannian manyfold (of Euclidian metrics) of the dimension $`𝒟>2`$, it is necessary to add something to the pure Euclidean gravity. However, the investigation of the latter seems sensible, because mechanisms that appear in this model may play a role in more realistic models. From this point of view the most interesting is the mechanism of compactification. We started from the model, which has $`\frac{D(D1)}{2}=45`$ local degrees of freedom correspondent to the independent components of the metric tensor $`g`$. Finally we arrive (in the elongated phase) at the system, which has the fractal dimension close to two. We expect that similar to the lower dimensional cases this system can be effectively approximated by a tree graph. In the correspondent branched polymer model  the lengths of the linear pieces of the graph and the way of their gluing together are the dynamical variables. What does then happen with the remaining degrees of freedom? They correspond to the small fluctuations of the pieces of the original manyfold represented by the links of the graph. Those fluctuations, obviously, become the internal degrees of freedom living on the polymer. This dimensional reduction should, in principal, resemble the dimensional reduction of the Kaluza - Klein models . However, the concrete realization of the mechanism has to be investigated. We expect that more realistic realizations of the dynamical triangulation approach should possess the dimensional reduction, which may be similar to that of the pure Euclidean gravity.

The investigation of the $`5`$ \- dimensional model (see ) indicates that it has a complicated phase structure. Namely, at $`\kappa _{D2}<5`$ several different vacua appear. And the system possesses tunnelling between them. In our research of the $`10`$ \- dimensional gravity we restricted ourselves by the ranges of $`\kappa _8(0.1,0.1)`$, where the phase transition between crumpled and elongated phases is observed. However, taking into account the above mentioned property of the $`5D`$ model we expect that the phase structure of the considered model is, probably, not limited by the observed two phases. Another open question is the order of the observed phase transition. Although we have some indications that it is of the first order, the complete investigation of the subject has not been performed.

## 5 Conclusions.

In this paper we report our results on the numerical investigation of the $`10`$ \- dimensional Euclidean quantum gravity in the framework of the dynamical triangulation approach. Our summary is as follows.

1. The considered model contains the phases that we call crumpled and elongated. The arguments in favor of view that those phases resemble the correspondent phases of lower dimensional models were obtained.

2. For the observed volumes the phase transition between the mentioned phases corresponds to small negative values of $`\kappa _8`$. Therefore at the ”physical” positive gravitational coupling $`G\frac{1}{\kappa _8}`$ the model does not contain the crumpled phase contrary to the lower dimensional cases, where critical value of $`\kappa _{D2}`$ is positive. This also means that the natural state of the collection of $`10`$ D simplices (that appears when the Einstein - Hilbert action is turned off and $`\kappa _8=0`$) belongs to the elongated phase (at least for the observed volumes). This is in the contrast with the behavior of the low dimensional models: at $`D=3,4`$ and $`\kappa _{D2}=0`$ the systems exist in the crumpled phase.

3. The average scalar curvature appears to be positive at the observed gravitational couplings ($`\kappa _8(0.1,0.1)`$) in contrast with, say, the $`4`$ \- dimensional case, where it becomes negative at $`\kappa _2`$ close to $`0`$.

4. The fractal dimension of the model is close to $`2`$ in the elongated phase for both considered volumes ($`V=8000`$ and $`V=32000`$). In the crumpled phase the fractal dimension appears to be $`3.2`$ for $`V=8000`$ and $`𝒟4`$ for $`V=32000`$.

5. The linear size of the system in the crumpled phase appears to be close to that of the lower dimensional models (the linear extent is $`10`$). In the elongated phase the mean linear extent and its fluctuations appear to be sufficiently larger than in the crumpled phase.

6. It is expected that the considered phase transition is of the first order. However, the investigation of this subject has not been performed. So the question about the order of the phase transition remains open.

## 6 Acknowledgments.

We are grateful to E.T Akhmedov, B.L.G. Bakker, M.I.Polikarpov and J.Smit for useful discussions. We are very obliged to J.Ambjorn and S.Catterall who have sent us the dynamical triangulation (Fortran and C) program codes. We consulted these codes several times during the elaboration of our own program codes. We are very obliged to the staff of the Joint Supercomputer Center at Moscow and especially to A.V. Zabrodin for the help in computations on supercomputer MVS 1000M. This work was partly supported by Federal Program of Russian Ministry of Industry, Science and Technology N 40.052.1.1.1112 and by the grants RFBR 03-02-16941, RFBR 04-02-16079, RFBR 02-02-17308, RFBR 01-02-117456, RFBR-DFG 03-02-04016, INTAS 00-00111 and CRDF award RP1-2364-MO-02.


text:
# Bose-Einstein Condensation of Molecular Hydrogen in Nanotube Bundles

## Abstract

We evaluate the effects of heterogeneity on the density of states of H2 molecules inside interstitial channels within bundles of carbon nanotubes. As temperature (T) falls, the density increases within those tubes having the greatest binding energy. At T $``$ 10 mK, the molecules undergo Bose-Einstein condensation, exhibiting a singular heat capacity.

The subject of adsorption within bundles of carbon nanotubes has received considerable attention recently, owing to both its fundamental interest and potential applications (e.g. gas storage and isotope separation). One focus of the research is one-dimensional (1D) and quasi-1D phases of matter. These include condensing and crystallizing phases of buckyballs within tubes, He and H2 within the interstitial channels (IC’s) between tubes and various gases within grooves on the outside of the nanotube bundles qsiev ; bob ; cro ; oscCv ; rmp ; aldo ; boro ; c60 ; bienf ; carlo ; hecond . Most analyses (including those of our group) have assumed that the tubes are identical and parallel, forming an ordered lattice. Real nanotube bundles, in contrast, consist of a disordered array of tubes with a distribution of radii. A logical question arises: how reliable are predictions that ignore such a variable environment? Shi and Johnson have recently shown that predictions incorporating heterogeneity agree better with adsorption data than do the idealized models hete . Stimulated by their work, we have explored the behavior of quantum fluids (He and H2) in such a nonuniform environment at low temperature (T). In this paper, we describe an intriguing result: Bose-Einstein condensation (BEC) of H2 molecules occurs as a consequence of the heterogeneity. This paper makes predictions about this phenomenon that are testable experimentally. Similar behavior is expected for 4He atoms. We note that a superfluid phase of para-H2 in confined geometries has been proposed, but experimental evidence of it is lacking thus far h2bec .

The occurrence and properties of the BEC transition are determined by the density of states $`N(E)`$ of the constituent particles. Thus, an initial focus is the determination of that function, derived from the single particle energy spectrum. We assume that particles do not interact (based on the assumption of low density and weak inter-IC interactions). A typical IC presents a highly confining geometry for the molecules if the bundle is close-packed. Because of this confinement (in the x-y plane, perpendicular to the tubes), the molecules have a large zero-point energy of motion. Theoretical values interc  of this energy are of the order of 500 K. This result is consistent with an enormous isotopic heat difference at low coverage (which equals the binding energy difference) between D2 and H2 in a nanotube bundle observed experimentally oscar . That difference was found to be about 250 K, a factor $``$ four larger than the difference found on the graphite surface vid  (due to less localization in that case).

Our analysis assumes that molecules can move within the IC’s in order to achieve chemical equilibrium with a coexisting vapor phase. At low T, because of the confinement, only the lowest energy state of transverse motion is excited. We call this state’s energy $`E_t(𝐑)`$, where $`𝐑=(R_1,R_2,R_3)`$ is a vector whose components are the radii of the tubes surrounding a particular IC. We have evaluated $`E_t(𝐑)`$, with interesting results. The calculation assumes that motion parallel to the IC is that of free molecules, so the total energy of a particle with (z component of) momentum p is given by $`E(p,𝐑)=E_t(𝐑)+p^2/(2m)`$. If the potential is corrugated, one must replace the particle’s mass $`m`$ with its band mass hecond ; band . The density of states for H2 is obtained by summing over the IC’s present in the given sample and integrating over p:

$`N(E)`$ $`=`$ $`{\displaystyle \underset{𝐑,p}{}}\delta [EE(p,𝐑)]`$ (1)
$`=`$ $`{\displaystyle \frac{L}{\mathrm{}\pi }}\left({\displaystyle \frac{m}{2}}\right)^{1/2}{\displaystyle _0^E}𝑑E_t{\displaystyle \frac{g(E_t)}{\sqrt{EE_t}}}`$
Here $`g(E_t)`$ is the density of states for the transverse oscillation problem and $`L`$ is the length of the tubes. If the length is not constant, the variation can be included in $`g(E_t)`$.

The energy $`E_t(𝐑)`$ is evaluated from the potential energy $`V(𝐫,𝐑)`$ of the molecule at position r. To compute $`V`$, we add contributions from the three neighboring tubes, ignoring corrections from more distant tubes and many-body effects associated with the screening of van der Waals interactions by the adjacent tubes milen . Although these approximations introduce some errors in the values of $`E_t`$ they do not affect the key predictions of this work, which are sensitive to the variation of $`E_t(𝐑)`$. The potential from each neighboring tube was derived with the method of Stan et al uptake . Because $`V(𝐫,𝐑)`$ varies rapidly with r, values of $`E_t(𝐑)`$ include large anharmonic and small anisotropic contributions.

The form of $`N(E)`$ for a given collection of nanotube bundles depends on sample preparation. A sample is represented by an ensemble of points (one for each IC) in R space. The density of points in R space, a function $`f(𝐑)`$, is defined so that $`f(𝐑)d𝐑`$ is the number of IC’s within an infinitesimal volume $`d𝐑=dR_1dR_2dR_3`$, centered at R. The function $`f`$ enters the transverse density of states through this expression:

$$g(E)=𝑑𝐑f(𝐑)\delta [EE_t(𝐑)]$$
(2)
At very low T, we need the value of $`f(𝐑)`$ only in the immediate vicinity of the global minimum of the energy ($`E_m`$ at $`𝐑_m`$), but at higher T the specific experimental distribution affects the results quantitatively. Here, we have assumed that the IC’s are uniformly distributed in R space within a radius spread of width 3 Å near $`𝐑_m`$. More general results will be described in a complete report of this work . We find, as one might expect, that $`𝐑_m`$ occurs along the diagonal, symmetry line ($`R_1=R_2=R_3`$). Along this line, which we call the (1,1,1) line, there occurs a global minimum $`E_m=1052.97`$ K at $`R_i`$=9.95 Å. To derive $`g(E)`$, it is important to know the variation of $`E_t(𝐑)`$ near this minimum. Consider a particular displacement (in the $`R_2R_3`$ plane) from $`𝐑_m`$ to a neighboring point for which $`R_1`$ has the same value, while $`R_2`$ and $`R_3`$ are slightly different: $`R_2=R_1+\delta `$ and $`R_3=R_1\delta `$. This change, parallel to the (0,1,-1) direction, yields an extremely small increase in the H2 energy, indicative of a very slow variation of the function $`E_t(𝐑)`$ near $`𝐑_m`$. This behavior, shown in Fig.1 reveals a long valley of low energy states in this direction; the same behavior occurs along the five equivalent directions, e.g. (1,0,-1).

The transverse density of states $`g(E)`$, from Eq. 2, is shown in Fig. 2. Note that $`g(E)`$ is proportional to $`\sqrt{EE_m}`$ for small $`EE_m`$; the prefactor is determined by the principal axes of curvature of the function $`E_t(𝐑)`$. This square root behavior is identical to that found near a van Hove singularity in the phonon density of states of a 3D system near a minimum in the Brillouin zone, for the same reason- phase space topology ziman . At higher energy, instead, the behavior of $`g(E)`$ switches to $`1/\sqrt{EE_m}`$. This (-1/2) power law follows from the fact that the higher energy displacements from $`𝐑_m`$ are quasi-1D. As seen in Fig.1, the iso-energy contours are perpendicular to the diagonal, so that the energy gradient in R space is along the diagonal, with essentially constant transverse variation, a 1D situation.

Fig. 2 also presents $`N(E)`$, derived with Eq. 1, which convolutes the transverse spectrum with the 1D motion along the axis. The resulting power law behavior can be understood from realizing that if $`g(E)`$ is proportional to $`(EE_m)^n`$, for some $`n`$, then $`N(E)`$ is proportional to $`(EE_m)^{n+1/2}`$ . Hence, we find that $`N(E)`$ is proportional to $`(EE_m)`$ near threshold. This linear behavior is that characteristic of a 4D gas in free space; the result implies that this system exhibits 4D gas behavior at low T. For $`E>20`$ mK above threshold, instead, $`N(E)`$ becomes approximately constant, corresponding to the density of states of a 2D gas. Thus, there arises a dimensionality crossover originating from the anomalous transverse density of states. We emphasize that the 4D regime is a direct consequence of the existence of a minimum in the function $`E_t(𝐑)`$, a result that is not sensitive to the details of the calculation groove .

The thermodynamic behavior of the system is derived with the usual bose gas theory. For a given number of molecules, $`N`$, the chemical potential $`\mu `$ is determined from the relation

$$N=𝑑E\frac{N(E)}{e^{\beta (E\mu )}1}$$
(3)
At a specified $`T=1/(k_B\beta )`$, this relation yields a maximum value $`N_{max}`$ when $`\mu `$ equals the lowest energy of the system, $`E_m`$, at which point BEC begins. That is, a macroscopic fraction of the molecules fall into the lowest energy state when $`N>N_{max}`$. Equivalently, at fixed $`N`$, BEC occurs when T falls below the inverse function $`T_c=T(N_{max})`$. The resulting dependence on $`N`$ of $`T_c`$ is shown in Fig. 3. As seen there, $`T_c`$ is of order 10 mK, which is experimentally accessible. Fig. 4 shows the specific heat $`C_N(T)`$, calculated from the energy as a function of $`\mu `$ and T. The novel behavior observed in the figure is a result of the unusual form of $`N(E)`$. At low T, $`C_N`$ is proportional to $`T^2`$, a consequence of the 4D (linear) variation of $`N(E)`$ at low $`E`$. Note the presence of singular behavior of $`C_N(T)`$ as one approaches the transition from above and a cusp at $`T_c`$ itself. At relatively high T, $`C_N(T)/(Nk_B)`$ is essentially unity because $`N(E)`$ is 2D-like at high $`E`$; a nondegenerate 2D gas has $`C_N(T)/(Nk_B)=1`$. Behavior for $`T>0.2`$ K, not shown, is very sensitive to the distribution of nanotube radii.

Discussions of these results with colleagues have led to several questions concerning the transition. One is this: since each IC contributes a density of states characteristic of a 1D system, why does non-1D behavior arise here? The answer is that particles can exchange between different IC’s because of their common access to the vapor region. In practice, this may be a very slow process, leading to nonequilibrium behavior. The presence of breaks or holes in the tubes may alleviate this kinetic problem without invalidating the model, since heterogeneity is an essential aspect of the model. Another question is whether similar behavior occurs for other bose gases. Indeed, BEC of 4He is predicted by a similar analysis to occur if the sample’s distribution of tubes includes those with some near its energy minima in R space, which occur near 8.5 Å. Finally, one might wonder about the effects of interparticle interactions, which have been ignored up to this point he4 . Indeed, some previous studies of H2 in IC’s have found that a liquid-vapor condensation occurs at a higher temperature ($``$ K) (in the absence of heterogeneity) hecond . A very recent study, however, found that nanotubes’ screening of the intermolecular interaction reduces $`T_c`$ to about 10 mK milen . However, that calculation omitted the role of heterogeneity, which is relevant, according to the Harris criterion, since the 3D specific heat critical exponent is positive harris . We expect that this condensation temperature is further reduced by disorder, enabling the BEC transition to occur.

We summarize our results as follows. Heterogeneity alters the qualitative behavior of the low energy spectrum of H2 molecules. The lowest-lying states of the system are those of particles in that channel. As T falls, particles aggregate in the (essentially 4D) space of quantum states, (R,p), with bose statistics having a dramatic effect, i.e. BEC, below a transition temperature of order 20 mK. Anomalous behavior is predicted for the specific heat, a consequence of the unusual density of states, which is 4D-like at very low energy and 2D-like at somewhat higher energy. An experimental probe of the real-space molecular density should reveal the needle-like concentration, below $`T_c`$, of a macroscopic fraction of the particles within the lowest energy channel.

Most intriguing to us is that this transition is a direct consequence of disorder, since the perfectly uniform system of identical nanotubes yields strictly 1D, nonsingular behavior. Such a dramatic effect of heterogeneity occurs elsewhere in low temperature physics. Examples include the spin-glass transition schiffer , the effect of tunneling states on thermal behavior of glasses glasses  and the effects of disorder on monolayer films hedisorder .

This research has benefited from discussions with Moses Chan, Vin Crespi, Susana Hernández, Jainendra Jain, Peter Schiffer, Paul Sokol, Flavio Toigo and David Weiss and support from NSF and the Hydrogen Storage program at Penn State. F.A. acknowledges funding from MIUR-COFIN 2001.


text:
# Thermodynamics of DNA packaging inside a viral capsid: the role of DNA intrinsic thickness

## I Introduction

The goal of the present study is to characterize the equilibrium thermodynamics of a thick polymer chain confined in a spherical cavity. In particular we aim at establishing what force is necessary to apply, and how much energy is to be spent, to achieve a given packing density for the polymer. This problem is intimately related to the problem of DNA packaging inside a viral capsid phi29review ; bloomfield ; reviewcapsid , which is under an increasing theoretical and experimental attention phi29 ; phi29motor ; gelbart ; odjik ; pais ; muthukumar ; podgornik . At variance with previous studies we model explicitly the double-stranded (ds) DNA as a cylindrical tube with a finite thickness (as opposed to a beads-and-springs model) and study how it impacts on the force and work required to package it. We shall further discuss how a finite thickness induces a non-negligible persistence length in a polymer. The experimental reference for the present study is provided by a series of experiments phi29  where the loading process of DNA inside the $`\varphi `$29 virus head was followed under controlled conditions. In this situation, the 6.6 $`\mu `$m long double-stranded DNA polymer was fed by a molecular motor inside the virus head, a prolate $`54`$ nm by $`42`$ nm icosahedronphi29 ; cell ; bloomfield . By using suitable force-feedback measurements it has been possible to carefully measure the force exerted by the motor at various stages of the loading process, and hence the total energy required for the whole packaging process. It has been found phi29  that the internal force exerted by the packaged DNA against the $`\varphi 29`$ capsid walls (which has to be overcome by the portal motor), displays a dramatic increase during the final loading stages. In particular, the maximum value attained by this internal force is $`50`$ pN. This is a strikingly large number when compared to typical forces encountered at the molecular and cellular scale.

In a recent studygelbart , a molecular dynamics simulation has been used to characterize the dynamical packaging of DNA in a spherical cavity. In that approach, the double stranded DNA was coarse-grained in a beads-and-strings chain subject to self-interaction. Interestingly, it was found that, at high packing densities, the resistance to packaging was not dominated by the DNA self-interaction (whether attractive or repulsive). An insight into this result can be obtained by a soluble mean-field model, presented later, which highlights how the divergence of the packing force at high density is not controlled by self-interactions.

These results stimulate the search for an alternative physical mechanism responsible for the onset of the large packaging pressure encountered e.g. in $`\varphi `$29. For this reason, we adopted a minimal model where DNA is treated as a thick chain of equispaced beads. Since the DNA self-attraction or repulsion does not appear to be crucial at high packing densities, we have intentionally kept the model at the simplest level and thus avoided the introduction of any explicit treatment of the known self-interaction terms strey . The only two parameters entering in our model are the spacing between the beads and the chain thickness; we do not attempt any tuning of these parameters but, instead, adopt the corresponding values desumed from the experimental literature. In particular, the natural base spacing is set to 0.34 nm while the radius of hydrated double-stranded DNA is taken as 1.25 nm bloomfield . Note that in this work the thickness of a tube is defined as the radius of its circular section taken normally with respect to the tube centerline.

As we shall show below, the appropriate treatment of the finite thickness of the biopolymer allows to account for a variety of elastic properties of DNA in a natural way,even in the absence of additional phenomenological parameters, such as the bending rigidity, that are used in standard DNA coarse-grained models. For example, it will be shown that the mere difference in thickness of single-stranded and double-stranded DNA can account for the large change in the respective persistence lengths. More importantly, by considering the force necessary to compactify thick chains of increasing number of beads, we estimate the packaging force (and work) which turns out to be in good quantitative agreement with experimental measures for $`\varphi `$29 under the same density conditions. The present analysis involves finite-size scaling techniques to extrapolate results obtained in equilibrium stochastic simulations of chains with up to 200 beads. The most compact conformations obtained for such lengths reveal the clear tendency of thick polymers to occupy the spherical cavity by adopting a spool-like conformation, consistently with previous expectations bloomfield ; podgornik ; odjik  and with the recent numerical study of Kindt et al. with a different model gelbart .

## II Theory and Methods

Due to the large number of atomic constituents of both inorganic and biological polymers, it is necessary to resort to suitably simplified models in order to characterize their physical behaviour by means of analytical or numerical studies. Customarily the polymeric chain is coarse-grained into a succession of discrete beads. Besides capturing the effect of chain connectivity through the introduction of suitable strings or springs between consecutive beads, it is usually necessary to introduce appropriate interactions between the beads in order to capture the salient physical features of the given system. For the case of ds DNA the latter would include effects such as strands self-interaction (either attractive or repulsive) and bending rigidity. As a useful term of comparison for the present and previous numerical studies of DNA packaging, it is useful to consider a mean-field model (inspired by the treatment in Refs.garel ; lise ) which, while incorporating the above mentioned features, is still amenable to an analytic treatment. This model consists of a polymer chain of $`N`$ beads embedded in a discrete three-dimensional space and subject to a contact interaction, $`ϵ`$, and a bending rigidity (corner penalty), $`h`$. Within a mean-field picture garel ; lise  it is possible to obtain explicit expressions for the free energy, $`F_N`$ and loading force, $`f`$, required to package the chain inside a cube of side $`L`$:

$`{\displaystyle \frac{\beta F_N(x)}{N}}={\displaystyle \frac{1x}{x}}\mathrm{log}(1x)+\beta ϵ(13x)`$ (1)
$`\mathrm{log}\left[{\displaystyle \frac{2+4\mathrm{exp}(\beta h)}{e}}\right],`$
$`a\beta f(x)=\mathrm{log}(1x)+\beta ϵ(16x)`$ (2)
$`\mathrm{log}\left[2+4\mathrm{exp}(\beta h)\right]`$
where $`x`$ is the density $`x=N/L^3`$, $`\beta =(k_BT)^1`$ with $`k_B`$ the Boltzmann constant and $`a`$ is the unit length in the system. It is important to notice that the free-energy per bead, $`F_N/N`$ depends on $`N`$ and $`L`$ only through the density $`x`$; this fact will be exploited to analyze the results for our confined off-lattice chains with finite thickness.

Despite the simplified nature of the mean field model, the approximate force of expression (II) provides a useful starting point for both understanding recent packaging studies as well as suggesting new improvements of traditional models. The crucial observation is that in the tight-packing limit, $`x1`$, the divergence of $`f(x)`$ is entirely controlled by the logarithmic term which is independent of both $`ϵ`$ and $`h`$. This result provides a theoretical framework for explaining the numerical findings of Ref.gelbart  who ascertained that the force loading curve was not too sensitive to the potential strength and sign at high packing densities. In other words, from this simple mean-field treatment we have an indication that a resistance to dense confinement comes from entropic effects rather then polymer self-interactions or flexibility. This is also in agreement with the picture proposed in Ref. strey .

To the best of our knowledge the fact that also the finite thickness of DNA is responsible for its elastic properties has not been fully investigated before, and hence is the focus of the present study. Usually, in beads-and-springs models, excluded volume effects are incorporated at the level of hard-core repulsion between pairs of coarse-grained beads. However, the most satisfactory way to account for steric effects is through the modelling of the polymer chain as a tube with finite uniform thickness, $`\mathrm{\Delta }`$, rather then a series of beads connected by one-dimensional springs. The appropriate treatment of the chain thickness usually leads to remarkably different physics than with traditional beads-and-springs models. A striking example is provided by the natural emergence of protein-like secondary motifs in marginally-compact conformations of thick chainsopthelix .

The finite thickness impacts on two distinct features of polymer conformations. On one hand it will constrain the local radius of curvature to be not less than $`\mathrm{\Delta }`$ to avoid singularities. This effect induces naturally a bending rigidity on the polymer. As will be shown later, the resulting persistence length depends crucially on $`\mathrm{\Delta }`$, so that the latter appears to be a decisive factor for the huge difference of the persistence length for single-stranded and double stranded DNA, as we shall show later. On the other hand there is also a non-local effect due to the fact that any two portions of the tube, at a finite arclength separation, cannot interpenetrate. Equivalently, the centerlines of the two portions need to be at a distance greater than 2$`\mathrm{\Delta }`$ (see Fig. 1).

In traditional beads-and-springs models it is only this second effect that is taken into account through pairwise hard-core repulsion. Interestingly, one needs to go beyond pairwise interactions to account for the above mentioned effects in discretised polymer chains. This is best illustrated considering the discrete succession of equispaced ”beads”, ($`\stackrel{}{r}_1,\stackrel{}{r}_2,\mathrm{}\stackrel{}{r}_N`$), as the centerline of a tube with thickness $`\mathrm{\Delta }`$. The requirement on the local radius of curvature is simply enforced by finding the radii of the circles going through any consecutive triplet of points and ensuring that each of them is greater than $`\mathrm{\Delta }`$. The non-local effect can be addressed within the same framework by considering the minimum radius among circles going through any non-consecutive triplet of points. This radius is precisely the distance of minimum approach and, again, has to be greater than $`\mathrm{\Delta }`$. In summary, the appropriate treatment for the finite thickness, $`\mathrm{\Delta }`$, of the tube associated to a discretized centerline is reflected in the fact that the radii of the circles going through any triplet of distinct points has to be greater than $`\mathrm{\Delta }`$Thickness1 .

One is thus naturally led to consider the following Hamiltonian for the unconstrained, non self-interacting, double-stranded DNA:

$$(\mathrm{\Gamma })=\underset{ijk}{}V_3(r_{ijk})$$
(3)
where $`V_3`$ is the three-body potential used to enforce the thickness $`\mathrm{\Delta }`$ of the chain Thickness1 ; Thickness2 ; jstatphys ; opthelix ; Thickness3 ; secstr  and $`\mathrm{\Gamma }`$ is the configuration of the ds DNA. The argument of $`V_3`$ is the radius of the circle going through the triplet of distinct points $`i,j,k`$ and has the form

$$V_3(r)=\{\begin{array}{cc}0\hfill & \text{if }r>\mathrm{\Delta }\text{,}\hfill \\ +\mathrm{}\hfill & \text{otherwise.}\hfill \end{array}$$
(4)
The tube ends were modelled as hemispheres by introducing two “phantom” beads essentially coinciding with the terminal beads of the tube.

Since our experimental reference case is the packaging of genome inside the $`\varphi `$29 capsid, we take $`\mathrm{\Delta }`$ as the thickness appropriate for hydrated double-stranded DNA: $`1.25`$ nm while the spacing between the beads of our discretized tube is naturally taken as the average base spacing: $`a=`$0.34 nm.

As stated earlier our goal is to describe the equilibrium thermodynamics of a self-avoiding thick chain of $`N`$ beads confined in a sphere of radius $`R`$. In particular, we aim at calculating the free energy $`F_N(R)`$, for which it is straightforward to calculate the work done to package it under conditions of no energy dissipation, and also the quantity:

$$\mathrm{\Delta }G_N(R)=F_N(R)F_N(\mathrm{}).$$
(5)
The force that is necessary to apply in order to feed additional chain beads inside the sphere can be calculated by differentiating $`\mathrm{\Delta }G_N(R)`$ with respect to $`N`$:

$$f\frac{\mathrm{\Delta }G_N(R)}{N}.$$
(6)
In principle, the free energy, $`F_N(R)`$ could be obtained through a Monte Carlo simulation where the conformation space is restricted to only those structures that can be confined in a sphere of radius not larger than $`R`$. This is, however, impractical due to the fact that in situations of high density (small $`R`$) it will be very difficult to move through distinct compact structures while respecting the confining constraints.

A useful alternative is to work in the conjugated ensemble where, instead of keeping fixed the radius of the confining sphere, $`R`$, one applies a uniform “hydrostatic” pressure, $`P`$ that compactifies the chain. In this situation, the Boltzmann weight of a configuration, $`\mathrm{\Gamma }`$ is given by:

$$e^{(E+PV)/K_BT}$$
(7)
where the temperature $`T`$ is taken as the room temperature, $`T=300K`$, $`E`$ is the structure energy calculated via Eq. (3) and $`V`$ is the volume of the spherical hull associated to $`\mathrm{\Gamma }`$. To find the radius of the spherical hull enclosing $`\mathrm{\Gamma }`$ we first determine the maximum distance, $`\overline{r}`$ from the centre of mass of any of the beads in $`\mathrm{\Gamma }`$. Since these beads lie on the centerline of a tube of thickness $`\mathrm{\Delta }`$, the required radius is $`R=\overline{r}+\mathrm{\Delta }`$. It should be noted that, in principle, the smallest confining sphere is not necessarily centred on the centre of mass of $`\mathrm{\Gamma }`$. However, we have verified by an explicit numerical search of the location of the centre of the smallest sphere that this approximation is very good when dealing with structures with good overall compactness, which are the focus of this investigation. We stress that $`E`$ can take only two values, 0, and $`\mathrm{}`$ for configurations that are respectively compatible or incompatible with the thickness requirement. This effectively restricts structure space to the ensemble of viable conformations of the thick tube.

By working at fixed temperature, one can sample the configuration space for various values of $`P`$ by using a standard Metropolis criterion. From the distribution (histogram) of the various hull volumes encountered at various pressures one can use an ordinary multiple histogram technique ferren  to recover the density of states, $`W(R)`$ corresponding to the number of viable structures (i.e. compatible with the preassigned thickness $`\mathrm{\Delta }`$) and hull radius $`R`$.

The required free energy, $`F(R)`$ is then:

$$e^{\beta F(R)}\underset{R^{}<R}{}W(R^{})$$
(8)
The reconstructed free energy, $`F(R)`$, is defined up to an additive constant, which reflects the fact that the density of states can be obtained only up to a multiplicative constant. This additive constant is immaterial for the calculation of the packaging work and force in a system with a given number of beads, though it has to be set appropriately in order to simplify the scaling analysis of $`F_N(R)`$ as a function of $`N`$ at fixed $`\beta `$. For this reason we have taken the unconstrained situation, $`R=\mathrm{}`$, as the reference case associated with zero free energy: $`F_N(R=\mathrm{})=0`$ (see Eq. 5).

It is worth to point out that the force estimate in (II) is for a chain that is fully confined in the cube, i.e. we have not considered the case of a chain that is partially contained in the cube while the remainder constitutes an external tail. The latter case may be closer to the in vivo condition; however, the contribution of the floating tail to the free energy is expected to become less and less important compared to that of the confined portion as $`x1`$ in Eqs. 1 and 2, which is the limit case of interest. In addition, in experiments where the packaging force is probed, the protruding tail is kept under considerable mechanical tension and therefore is expected not to make significant contribution to the system entropy.

## III Results

One of the paradigm models fruitfully used to describe DNA chains is provided by the freely-rotating chain Daune . Within this scheme a DNA chain is described as a succession of equally-long bond vectors, $`\stackrel{}{t}_i`$, subject to the constraint that consecutive bonds must be at a given angle, $`\alpha `$. This constraint directly impacts on the chain persistence length (which measures how correlations between two bond vectors decrease with their arclength separation).

For a given length of the bond vectors, the bonding angle $`\alpha `$ has to be set a posteriori so to reproduce e.g. the experimental persistence length, $`\xi `$. It is appealing to notice, as pointed out before, that the notion of intrinsic thickness of the DNA chain, naturally impacts on the minimum attainable radius of local curvature and hence on the allowed range of the bonding angle, $`\alpha `$. This leads to the existence of a non-trivial persistence length, $`\xi `$, induced by the chain thickness. Such dependence can be estimated as follows. Consider a discrete chain of thickness $`\mathrm{\Delta }`$, with the consecutive unitary bond vectors $`\stackrel{}{t}_i\stackrel{}{r}_{i+1}\stackrel{}{r}_i`$ and $`\stackrel{}{t}_{i+1}`$. The constraint over the bonding angle is give by

$`\stackrel{}{t}_i\stackrel{}{t}_{i+1}`$ $``$ $`1{\displaystyle \frac{1}{2\mathrm{\Delta }^2}}.`$ (9)
Neglecting the effects of the chain self-interactions it is possible to calculate explicitly the quantity $`\stackrel{}{t}_n\stackrel{}{t}_m`$, with $`n<m`$, which defines in turn a persistence length, $`\xi `$ via $`\stackrel{}{t}_n\stackrel{}{t}_m\mathrm{exp}\left(\frac{|nm|}{\xi }\right)`$. One thus obtains that $`\xi (\mathrm{\Delta })`$ is given by

$$\xi (\mathrm{\Delta })=\frac{1}{\mathrm{log}\left(1(2\mathrm{\Delta })^2\right)}.$$
(10)
This approximate relationship can be specialized to the case of double-stranded DNA, where the unit bonds are taken as the base spacing, 0.34 nm and the thickness equal to the radius of hydrated dsDNA, 1.25 nm. One thus observes $`\xi =19`$nm, which is about one half of the correct experimental value $`\xi _{exp}50`$ nm. This theoretical estimate, which takes into account the mere effect of steric hindrance of the chain does not include the contribution of electrostatic repulsion of the phosphate groups to the apparent persistence lengthcozzarelli  but appears to capture the correct order of magnitude of the observed persistence length.

As a further validation of this approach it is worth considering the case of single-stranded DNA. In this case the vector length is equal to $``$ 5 Å, which is the typical nucleotide separation projected on a line interpolating the chain, while the thickness is about $``$ 4 Å. The near equality of these two parameters imply that the thickness induced persistence length of single-stranded DNA is comparable with the nucleotide spacing, and in fact one obtains, $`\xi `$ $``$ 1 nm. Again, this is consistent with the experimental data for ssDNA in solutions of high ionic strength (which screen the phosphate charges) for which values for $`\xi `$ of the order of 1 nm are observed plengthssdna1 .

Given the simplified derivation of the parametric dependence of $`\xi `$ on $`\mathrm{\Delta }`$ and the total absence of adjustable parameters in the theory, this agreement is remarkable and testifies the need to model explicitly the DNA intrinsic diameter.

While the soluble models leading to e.g. the results of eqns (II) and (10) are useful for gaining insight in the basic physics of thick polymers, they cannot be used to fully characterize more complicated instances such as the packaging process. To address this issue it is, therefore, necessary to resort to stochastic numerical simulations. In the present study, the equilibrium properties of the thick polymer was characterised by a multiple Markov chain simulationTROW  where several copies of the system each at a different pressure were evolved simultaneously. The evolution was controlled by the Metropolis acceptance of elementary chain distortion involving crankshaft, pivot and slithering movessokal . A preliminary measure of correlation times allowed to collect statistically independent measures for reconstructing the density of states as described in the methods section. The full equilibrium thermodynamics was thus obtained for tubes with a discrete number of beads, $`N`$. Values of $`N`$ up to 200 were considered here. As was argued before, in the limit of large $`N`$ the free energy should be extensive in $`N`$ at fixed filling fraction, $`N/R^3`$. Therefore we expect $`F_N(R)`$ to obey the following scaling relations:

$$\underset{N\mathrm{}}{lim}F_N(R)Ng\left(\frac{N}{R^3}\right),$$
(11)
which enables, for large $`N`$, to express the confining force as a function of $`g`$:

$$af_N(R)=\frac{F}{N}g\left(\frac{N}{R^3}\right)+\frac{N}{R^3}g^{}\left(x\right)|_{x=N/R^3}.$$
(12)
where $`a=0.34`$ nm is the unit length in our model. Note that the same scaling theory applies also to expressions (1) and (II). These relations will be obeyed only approximately for any finite chain, though such corrections decay for increasing values of $`N`$, as shown in the following.

In Fig. 3 we have plotted the reduced free energy, $`g_N=F_N(R)/N`$ as a function of $`N/R^3`$ for different values of $`N`$. It is apparent that, as the size of the system, $`N`$, increases, the curves for $`g_N`$ deviate less from each other, indicating the decrease of finite-size corrections from the limit curve of expression (11). Indeed, the systematic decay of such corrections for increasing $`N`$ is visible in Fig. 2 where we have reported the trend of $`g_N`$ as a function of $`N^1`$ at a fixed density. As visible in Fig. 3, the data are well interpolated by a linear regression with correlation coefficient $`r0.98`$ . Through such regression one can estimate the error on the extrapolated value for $`g`$ as the error on the intercept. The observed correlation coefficient is statistically significant, since it pertains to 16 data points and is suggestive that finite-size corrections decay approximately as $`N^1`$. However, from the systematic decrease of the linear regression slope observed when fitting data for the larger tube lengths only, it can be argued that the convergence is faster than $`1/N`$ and may also be compatible with decays such as $`N^{1.5}`$.

The same rate of convergence of finite-size corrections is inherited by the force obtained from (12) by using $`g_N`$ instead of $`g`$. By considering several values of $`N`$ one obtains a succession of curves for the loading force that, again, appear to approach a limiting curve from above. This is in agreement with the physical intuition that more force is required to pack a shorter (but equally thick) chain to the same target density.

As for the corrections to $`g`$, also the force curves appear to approach the limit curve of interest at a rate somewhat faster than $`1/N`$. Nevertheless, since we cannot give a precise account of the correction exponent we calculated the limit force by using the conservative $`1/N`$ decay estimate and by fitting only the data of the largest tube lengths, $`100<N200`$. The results are provided in Fig. 5 and can be compared with the experimental curve of ref. phi29  as discussed below.

If one considers the prolate $`\varphi `$29 capsid as an ellipsoid the reference density at full genome packaging as $`x=0.063`$. The largest densities reached in our simulation correspond therefore to about 85 % of this value. Neglecting the contribution of the protruding tail in the $`\varphi 29`$ experiment one can then compare the experimental force at 85 % genome packaging, 34 pN, with our peak force of 16 pN, which appear to be correct within only a factor ot two. We stress that this finding exploits only the severe restrictions in configuration space operated by the finite thickness of DNA itself; it is striking that this entropic effect appears to account for half the force observed in the experiment. Furthermore, the limit curve for the force has been obtained with the conservative $`1/N`$ estimate of the decay of finite-size corrections; this arguably leads to an underestimation of the force required to package the DNA chain.

Another quantitative comparison against the experiment can be done considering the total packaging work. The experimental total packaging work in the experiment of ref phi29  was estimated as 2 $`10^4K_BT`$, corresponding to $`1.04K_BT`$ per base pair (or bead in our model). In our case, the work per bead necessary to package the tube to a preassigned density, $`x`$, (starting from an unconstrained situation) is given simply by $`g(x)`$ (see Fig. 2). As before, we can connect the peak density reached in our simulation with the experimental situation where 85 % genome is packaged in the capsid. The value of $`g`$ observed in our model at this density corresponds to 0.55 $`K_BT`$, which differs by less than a factor of two from the experimental work (but the latter refers to 100 % packaging). In summary, even without accounting for the different final densities in the experiment and in our model, we can nevertheless conclude that also a substantial fraction (more than 50%) of the work required to package DNA could result from the entropic cost of compactifying a polymer of finite thickness. These results are also consistent, a posteriori, with the findings of Ref. gelbart , where it was pointed out that the detailed form (and even sign) of the DNA self-interaction was not the determinant factor in the packaging process. The important role played by the loss of configurational entropy in DNA packaging or condensation is also supported by recent numerical and analytical work harvey ; odijk93 .

As we argued before, the finite DNA thickness appears to account for a significant fraction of the physical properties (such as persistence length and resistance to confinement) that are usually ascribed entirely to DNA self-interaction and bending rigidity. As an extension of the model discussed so far, one is lead to consider the case where, on top of the finite-thickness treatment one adds bending rigidity effects. The latter are customarily taken into account by the following energy term:

$$E_{b.r.}=\frac{K_BT}{2}\xi 𝑑s\frac{1}{R(s)^2}$$
(13)
where the integral runs over the DNA contour, $`R`$ denotes the local curvature and $`\xi `$ is the persistence length. In a model which already incorporates DNA thickness the use of the full persistence length in expression (13) certainly overemphasizes the elastic effects (since thickness alone can satisfactorily account for a significant fraction of the DNA persistence length, $`\xi `$). Even in this situation it will be apparent that the additional (overestimated) force increment ascribable to the presence of bending rigidity is of the same order of the one induced by the thickness.

As a first step we provide an upper bound for the confining force increment of the due to the added elastic term. To this purpose let us consider a thick chain that has been compactified so to nearly fill a given cavity. It is assumed that packing conditions are such that there is almost no structural freedom left for how to add an additional bead inside the hull. In particular, the bead has to be added so to achieve the tightest local curvature of the chain. Since the lower bound for the curvature is the thickness itself, $`\mathrm{\Delta }`$, the additional packaging force is estimated as:

$$\mathrm{\Delta }f_{b.r.}=\frac{K_BT\xi }{2\mathrm{\Delta }^2}.$$
(14)
Using $`\xi =50nm`$ as in Ref. gelbart  and $`T=300K`$ one obtains that an upper bound for the force increment due to bending rigidity is $`\mathrm{\Delta }f_{b.r.}66`$ pN. To obtain a more accurate estimate of $`\mathrm{\Delta }f`$ we have undertaken additional studies where the numerical scheme previously discussed for the confined thick chain was generalised to include explicitly the bending rigidity term of eq. (13).

The augmented difficulty of packaging the chain in this situation resulted in the fact that the maximum density achieved was $`N/R^3=0.050`$, corresponding to 80 % of the estimated packing fraction of $`\varphi 29`$. In this case, the confining force for a chain of 200 beads resulted equal to 39 pN. This corresponds to an extra force of around $`12`$ pN, (the pure thickness-induced force at the same filling fraction is $``$ $`27`$ pN, see Fig. 4). This value overestimates the bending rigidity contribution in the thermodynamic limit due to the systematic decrease of the force as a function of chain length. However, even allowing for the overemphasized contribution of bending rigidity due to both finite-size effects and the use of the full pesistence length $`\xi `$ in eq. (13) one can already conclude that the finite-thickness effects are comparable with those associated with bending rigidity (see. Figs. 4 and 5).

This conclusion remains true when finite-size extrapolation is carried out on this second model using chains of 40, 80,100,120,160 and 200 beads. From the thermodynamic data a total force of $`24\pm 7`$ pN can be extrapolated at the highest densities $`N/R^3=0.050`$. Despite the still simplistic nature of the second model considered here, the extrapolated force appears to be fully consistent with the experimental results for phi29 at the corresponding filling fraction of 80 %, where a total force of 27 pN was recorded. Thus, the different calculations presented here all indicate that the contribution to the resistance to packaging arising from the finite DNA thickness is comparable to that of other terms usually considered in coarse-grained models of DNA. This indicates that various phenomenological properties of DNA can be accounted in a more accurate and complete way by including an explicit modelling of DNA thickness besides the traditional terms in the effective Hamiltonian.

Finally, it is interesting that, consistently with the study of ref. gelbart  and with other theoretical studiesodjik ; harvey  the most compact conformations obtained in our stochastic simulations (in thermodynamic equilibrium) show the marked tendency to occupy the spherical cavity by performing a succession of turns in a spool-like fashion with a portion of the tube, while the remainder of the chain fills the space around the spool axis, as visible in Fig. 6.

## IV Conclusions

We have characterized the equilibrium thermodynamics of a polymer chain with finite thickness confined inside a spherical cavity. The model aims at schematizing the packaging of the genome (double-stranded DNA) inside a bacteriophage capsid. Stochastic simulations and reweighting techniques were used to calculate the force (and work) required to compactify to a given density a discrete tube of up to 200 beads. The sole inputs of the model are the experimental values for the spacing of the base pairs in double-stranded DNA, and the thickness of the latter; no tunable parameters or specific forms of self-interaction of the polymer are introduced. By using finite-size scaling techniques we have estimated the force necessary to package the genome inside the $`\varphi `$29 virus. The predicted force (and packaging work) agrees with the experimental one to within a factor of two, thus testifying the importance of an appropriate modelling of the DNA thickness in the packaging process. This fact is further corroborated by a theoretical calculation which shows that the different thickness of single- and double-stranded DNA can satisfactorily account for the large difference of the respective persistence lengths.

ACKNOWLEDGMENTS We are grateful to S. C. Harvey, T. Odijk and A. Maritan, for illuminating discussions and suggestions. We acknowledge financial support from INFM and Murst Cofin 2001.


text:
# Reply to “Comment on ‘Properties of the massive Thirring model from the 𝐗𝐘𝐙 spin chain’ ”

(February 12, 2003)

## Abstract

We elaborate in more details why lattice calculation in \[Kolanovic et al, Phys. Rev. D 62, 025021 (2000)\] was done correctly and argue that incresing the number of sites is not expected to change our conclusions on the mass spectrum.

In the comment on our paper PP00  Fujita, Kobayashi and Takahashi claim that our results for the mass spectrum of the massive Thirring model (MTM) are not reliable FuKoTa03 . In particular, they claim that if one uses spin chain regularisation it is *necessary* to diagonalize the spin chain Hamiltonian with a number of sites $`N`$ larger than 1000. We now explain why we disagree with this criticism.

The criticism of Fujita et al. is based on the following argument. If one makes the standard spin chain regularisation of MTM (which is $`XYZ`$ spin 1/2 chain Lut ; Lusch ) to obtain some reasonable results on continuum extrapolation, one has to satisfy the condition

$$\frac{2\pi }{N}am_02\pi .$$
(1)
Here $`N`$ is the number of sites, $`a`$ is the lattice spacing, and $`m_0`$ is a bare mass parameter. Using (1) Fujita et al. claim that if one wants to obtain any reliable information on the bound state of the MTM one has to take $`N>1000`$, which is much larger than ours $`N16`$. Moreover, they claim that for values of the parameters used in PP00  the left inequality in (1) is even completely violated, i.e., $`2\pi /N>am_0`$. In addition they claim that the mass gap in our calculation is approximately equal to the “resolution” $`2\pi /L`$. From all this, Fujita et al. conclude that our results PP00  cannot be very reliable for the bound state spectrum of the massive Thirring model.

Let us now review standard lattice philosophy. Continuum regime in lattice calculations is obtained when the correlation length $`\xi `$ is much larger than the lattice spacing $`a`$ and at the same time much smaller than spatial extension $`L=Na`$, i.e.,

$$a\xi Na.$$
(2)
As $`\xi =1/M`$, where $`M`$ is the mass gap (mass of the lightest particle), (2) can be equivalently written as

$$\frac{1}{N}Ma1$$
(3)
Now, hardware limitations make restrictictions on $`N`$. For example, in lattice (quenched) QCD maximum lattices which are presently calculable have $`N64`$ (see e.g., aoki ). More importantly here, for exact diagonalisation in two dimensions $`N<30`$. It follows that “$``$” in (3) at best means 5–8 times smaller. So in practical calculations one effectively imposes the condition

$$\frac{1}{N}<Ma<1$$
(4)
and from the quality of scaling and the accuracy of the continuum extrapolation (usually using different methods as a check) one decides weather the extrapolated results are a good approximation of the continuum theory. Of course, one cannot exclude the possibility that for values of scaling parameter larger than accesed there is a complete change of scaling behavior so that the obtained extrapolated results are wrong111For nice discusion see Sec. 9.5 in Ref. ChrHen ..

Let us now apply the above textbook analysis to our lattice calculation ($`XYZ`$ spin chain regularisation of MTM) PP00 . For the sake of clarity we restrict ourself to a interval of coupling constant where elementary fermion is the particle with the lowest mass. First, it is easy to see that (3) is *not* equal to (1) used by Fujita et al. It was shown in Lusch  that (using notation from PP00 )

$$m_0a=\frac{8\gamma }{\pi }\mathrm{sin}\gamma \left(\frac{Ma}{4}\right)^{2\gamma /\pi }.$$

In particular, let us analyze the left inequality in (4). In our extrapolation we had $`Ma>0.2`$ which is for $`N=16`$ reasonably larger than $`1/N=0.06`$. It follows that in our analysis continuum condition is fairly satisfied when *proper* condition is used, contrary to claim of Fujita et al.

As we mentioned above, there is always the possibility that for larger $`N`$ something dramatic happens with the scaling law and that our extrapolations made with $`N16`$ are incorrect. There are several reasons why we believe that this should not be expected in the case of the $`XYZ`$ chain:

* global properties of the energy spectrum are as expected for QFT with calculated mass spectrum. In particular, there are states with energies corresponding to two-fermion and two-(first)breather states
* the same technique was previously succesfully applied for similar perturbed CFT’s, e.g. ordinary and tricritical Ising model in a magnetic field Henkel . It was shown there that passing from $`N14`$ to $`N21`$ only slightly improved precision and conclusions about the spectra remained the same. In fact it was shown that very accurate results can be obtained already for lattices with $`N24`$ CasHas .
* our results agree with the DHN mass formula dhn  obtained by a number of different methods. It is very hard to imagine how bad extrapolation can agree with exact analytic result. We should also mention that a similar analysis for lattice regularization of the sine-Gordon model (periodic $`XXZ`$ chain in transverse magnetic field) also gave results consistent with the DHN formula PP99 .
* Beside the mass ratios in $`L\mathrm{}`$ limit we obtained anomalous dimensions of corresponding states, in particular we had for the first time calculated dimension of the second breather. Our result was subsequently confirmed by analytical calculation FRT99 .

We conclude that our numerical analysis made in PP00  was done correctly and that all results should be trusted, in a sense that enlarging lattice would just increase numerical precision and leave all conclusions unchanged.

Finally, we would like to comment on the Bethe ansatz solution of the MTM. Spectrum of the MTM can be found by solving the Bethe ansatz equations in the continuum approximation bt . The spectrum found there agrees with the results of dhn . Recently, it was argued fu  that there exist no complex solution to the Bethe equations and that there is only one bound state. In bt  precisely the complex solutions (so called Bethe strings) are describing bound states of MTM. We performed numerical analysis of Bethe equations for MTM and tried to reproduce the result of fu . Our results indicate that the powers $`\alpha `$ that determine the scaling behavior of energies of the few lowest states on the density $`\rho `$ differ and even vary with the $`\rho `$ (for definitions see section 4 in fu ). Therefore we were unable to extract conclusive results by letting $`\rho \mathrm{}`$. An indication that the numerical iterations might easily miss the complex solutions (that nevertheless exist) is coming from our study of the Bethe ansatz equations for different spin chains strings . There we studied how the complex (string) solution in two and three particle sectors emerge and disappear when one changes the parameters of the model. We found that iterative numerical methods fail to converge on string solutions, although they exist and can be found analytically.

###### Acknowledgements.

MK would like to thank to J.H. Lowenstein for useful and stimulating discussion.