text:
# Fragment Charge Correlations and Spinodal Decomposition in Finite Nuclear Systems

## 1 Introduction

When enough energy is brought into a nucleus, it breaks up into smaller pieces: this is called multifragmentation MoWo93 ; Gro90 . The knowledge of this process, experimentally observed for many years, was recently strongly improved with the analysis of experiments performed with powerful detection arrays MDA00 ; Hau00 ; Elli00 ; Xu00 ; I29-Fra01 ; Bot01 ; MDA02 ; Radu02 ; Kle02 ; Ell02 . As a theoretical framework for describing this phenomenon, the analysis of the bulk dynamics of nuclear matter, under various conditions of density and internal energy, has found arguments that multifragmentation occurs when nuclear matter has expanded through the spinodal region of negative compressibility Ber83 . In this area of mechanical instabilities, covering a large part of a liquid-gas type coexistence domain, the irreversible growth of local density fluctuations leads the system to the separation into two phases: the spinodal decomposition. For the liquid part primary fragments of nearly equal sizes should be favoured, in relation to the wavelengths of the most unstable modes Ayi95 . Effects like beating of different modes, coalescence of nascent fragments, secondary decay of the excited fragments and, above all, finite size effects are expected to deeply blur this simple picture Jac96 ; Colo97 . However, evidence for an enhancement of events with nearly equal-sized fragments was for the first time observed in 32 AMeV 129Xe + natSn central collisions leading to a fused system which undergoes multifragmentation, by using a charge correlation method I31-Bor01 . By means of this model independent method, recently improved Des02 , we are investigating here the evolution of the equal-sized fragment partitions with the increase of the incident energy up to 50 AMeV, in selected samples of experimental fused events concerning the same system. In the same framework, for collisions at the lowest energy, the predictions of the 3D stochastic mean-field simulations are successfully compared to the experimental results; these simulations take into account the dynamics of the most unstable modes in the spinodal region.

The paper is organized as follows. In section 2 we firstly present the experimental set-up, including the detector array specifications and operating conditions during experiment. Secondly, we recall the criteria allowing to select the experimental events corresponding to fused systems in central collisions. The general features of the related experimental data are evidenced in section 3. In section 4, the charge correlation method with some improvements and developments is described; results of stochastic mean field simulations are used to exemplify and to compare the potentiality of the up-dated version to the original method. The experimental results, obtained at different incident energies, are then shown; they are discussed and interpreted in section 5. Conclusions of this study are drawn in section 6.

## 2 Experimental selection of a finite piece of nuclear matter

Heavy ion collisions at relative velocities comparable to those of the nucleons in the nucleus may provide test pieces of nuclear matter at moderate temperatures which are good candidates to undergo a nuclear liquid-gas type phase transition. Indeed, “fused” systems reaching sizes up to a few hundreds of nucleons may be obtained in the laboratory, by colliding accelerated heavy projectiles with heavy targets. By means of $`4\pi `$ detection arrays of high granularity, it is possible to completely study their disassembly. Such events are expected to appear with high probability for central collisions, but one major problem is to select them among reactions with a dominant binary character I28-Fra01 .

### 2.1 Experimental procedure

The 129Xe + natSn system was studied with the $`4\pi `$ multidetector INDRA, operating at the GANIL accelerator. A thin target of natural tin (350 $`\mu `$g/cm2) was bombarded by 129Xe projectiles at five incident energies: 25, 32, 39, 45 and 50 AMeV. Low intensity beams ($``$ 3 $`\times `$ 107 pps) were used to keep the random coincidence rate below 10-4. The low target thickness allows slow fragments to escape the target. A trigger based on the multiplicity was chosen, requiring at least four modules firing. In the off-line analysis, events having a multiplicity of correctly identified charged particles inferior to the experimental trigger condition were rejected for reasons of coherency. The most peripheral collisions are thus eliminated.

INDRA, which is described in detail in I3-Pou95 ; I5-Pou96 , can be viewed as an ensemble of 336 telescopes covering about 90% of the $`4\pi `$ solid angle. The detection cells are distributed amongst 17 rings centred on the beam axis. Low-energy identification thresholds (from $`0.7A`$MeV for $`Z=3`$ to $`1.7A`$MeV for $`Z=50`$) and large energy ranges were obtained through the design of three layer telescopes, composed of an axial-field ionization chamber operated at 30 mbar of $`C_3F_8`$, a 300 $`\mu `$m silicon detector and a CsI(Tl) scintillator, thick enough to stop all emitted particles, coupled to a phototube. Such a telescope can detect and identify from protons between 1 and 200 MeV to uranium ions of 4 GeV. Past 45, where fast projectile-like fragments are no longer expected, the telescopes comprise only two stages, the ionization chamber operated at 20 mbar and the scintillator. Finally the very forward angles (2 - 3) are occupied by NE102 - NE115 phoswiches.

A charge resolution of one unit was obtained for the whole range of atomic number of fragments identified through $`\mathrm{\Delta }EE`$ method in the Si - CsI(Tl) couple. For CsI(Tl) scintillators a better understanding of the light response was obtained and the contribution of $`\delta rays`$ generated by the incoming heavy ion was taken into account I33-Par02 . A direct consequence was the identification of fragment with a resolution of one charge unit up to $`Z=20`$, and a few charge units for the heaviest fragments, in ionization chamber - CsI(Tl) telescopes I34-Par02 . The exact identification of fragments up to at least $`Z=20`$ reveals essential for the charge correlation studies which are the aim of this paper.

### 2.2 Selection of single source events

A two step procedure has been used to isolate fused systems. The first step was to keep the events for which a quasi-complete detection of the reaction products has been achieved. Significant fractions: $`77\%`$ of the charge of the system, $`Z_{sys}=Z_{proj}+Z_{targ}`$, and $`75\%`$ of the beam momentum, $`P_{proj}`$, in the exit channel were required to be measured for every event. INDRA does not permit isotopic identification for fragments and does not detect the neutrons. For this reason, the momentum used here is calculated from the product of atomic number $`Z`$ and velocity component in the beam direction $`v_z`$: $`P_{tot}=\mathrm{\Sigma }Zv_z`$ and normalized to the incident (projectile) momentum: $`Z_{proj}v_{proj}`$. In the second step we used the flow angle ($`\mathrm{\Theta }_{flow}`$) selection Cug83 ; Lec94 ; I28-Fra01 . This global variable is defined as the angle between the beam axis and the preferred direction of emitted matter in each event. It is determined by the energy tensor calculated from fragment ($`Z5`$) momenta in the reaction centre of mass. Fused events have no memory of the entrance channel and should be isotropic while binary dissipative collisions are focused at small $`\mathrm{\Theta }_{flow}`$. Thus, by selecting only large flow angles, fused events can be well isolated. In I28-Fra01  the minimum flow angle chosen for 32 AMeV collisions was 70. It was 60 for 50 AMeV collisions I9-Mar97 . In this paper we chose 60 for all energies, to get enough statistics (at least 30000 events) without degrading the properties of a compact fused source.

The present selection corresponds to measured cross sections decreasing from 56 to 19 mb when the incident energy goes from 25 to 50 AMeV (see table 1). By taking into account detection efficiency and biases due to the selection (quasi-complete events and flow angle selection) the total cross section for the formation of compact fused systems is estimated to decrease from 250 to 85 mb between 25 and 50 AMeV I39-Hud02 .

## 3 Global properties of selected fused events

Before presenting and discussing charge correlations, which deal with information inside events, it is useful to have an overview of the related inclusive properties of the selected events, like the multiplicity and charge distributions. Entrance channel conditions for the Xe+Sn system of total mass $`A_{sys}=248`$ and charge $`Z_{sys}=104`$ are summarized in Table 1.

### 3.1 Multiplicity and bound charge

Fig. 1 shows multiplicity and total charge ($`Z_{\mathrm{bound}}`$ = $`_{ZZ_{min}}Z`$) distributions normalized to the number of events, at the five incident energies. The total multiplicity distributions, when all the charged reaction products are considered, $`Z1`$, are shown in fig. 1a. They all display a gaussian shape; the most probable multiplicity values and the standard deviations increase with the available energy, as shown in Table 1. At 50 AMeV the average multiplicity reaches one third of the total charge. The distributions of the sum of the charges of all the reaction products, plotted in fig. 1e), for different energies, are practically identical, nicely testifying about the negligible variation of the detection efficiency of INDRA when the incident energy significantly varies. The common limitation at 80 charge units is demanded by software as a criterion for the completeness of event detection. These two plots, considered together, speak about the increasing number of reaction products from central collisions with the incident energy.

The evolution of the reaction product multiplicity distribution, at all energies, when another lower charge limit is imposed, is presented afterwards in the left column of the same figure. As compared to fig. 1a ($`Z1`$), the splitting of the energy dependent distributions is diminishing if $`Z3`$ (fig. 1b); the curves become completely superimposable, except at 25 AMeV, for $`Z5`$ (fig. 1c) and split again, but in the opposite sense, if $`Z8`$ (fig. 1d). Deexcitation simulations for primary fragments having the same N/Z ratio as the total system, performed by means of the statistical code SIMON T18ADN98  at 3 AMeV excitation energies, have indicated small evaporation rates of Li and Be isotopes, but negligible rates of heavier nuclei T27Tab00 . This average value of 3 AMeV for the fragment excitation energy was deduced from fragment-particle correlations for the same system between 32 and 50 AMeV incident energy I11-Mar98 ; I39-Hud02 . Moreover the daughter nucleus atomic number is, on average, one charge unit smaller than that of the parent fragment, with a standard deviation of about one unit. It may be thus inferred that the distribution of the number of primary fragments with $`Z5`$ is not dramatically modified by the secondary decay. More precisely, as indicated in Table 1, the average number of fragments with $`Z5`$ reaches a maximum for an incident energy as low as 39 AMeV and starts decreasing above 45 AMeV.

If the lower limit is set to $`Z_{min}`$=3, for comparison with other works, the maximum value, 7.1 fragments, is reached at 50 AMeV. The same maximum value is obtained for the asymmetric Ni+Au collisions measured with INDRA I37-Bel02 . These results can be compared with a series of measurements of fragment multiplicities published in Sis01 . Using for the present data a similar selection method would increase the INDRA value by $``$10%. These INDRA data, combined with those related to Ni+Ni and Kr+Nb of Sis01  show that the maximum average number of fragments produced in central collisions (multifragmentation) is proportional to the total system mass Asys up to at least A$`{}_{sys}{}^{}`$250; there is no indication of a saturation, as inferred from the Kr+Au data of Wil97  and from percolation calculations. The maximum number of fragments for A$`{}_{sys}{}^{}`$250 is moreover reached for a c.m. available energy of $``$ 12 AMeV in the INDRA measurements, much lower than those quoted in Sis01 .

Pictures f), g), h) in the right column of fig. 1 show the distributions of the sum of the charges of the reaction products having atomic numbers $`Z`$ 3, 5, and 8, respectively. The independence of the energy, seen in fig. 1e), is progressively removed from top to bottom and the populated domain shifts towards lower values of the total charge of the reaction products taken into account. For $`Z5`$, for example, the shapes of the distribution are very alike: gaussians which are slightly broadening with increasing energy, but centred at lower and lower values of the sum $`Z_{\mathrm{bound}}`$ of the fragment charges. For $`Z8`$ and the largest incident energies (45 and 50 AMeV) shoulders are observed at low $`Z_{\mathrm{bound}}`$, which indicates the onset of events with a high degree of fragmentation: those events have only one fragment with a charge in the range 8-14.

The evolution of the exit channel average measured quantities with the available energy is also synthesized in Table 1. These average values concern: the total multiplicity, $`M_{tot}`$, the fragment ($`Z5`$) multiplicity, $`M_f`$, the total charge emitted in fragments, $`Z_{\mathrm{bound}}`$, and the fragment atomic number, $`Z_f`$. The growth of the average total multiplicity, when the incident energy increases, is accounted by the growth of the light nuclei ($`Z4`$) average multiplicities, while the average fragment multiplicity barely changes. The ratio of the bound charge to the total charge of the system goes down with increasing excitation energy. The size of the produced fragments is diminishing too, as a proof that the collisions are becoming more violent. At a given incident energy, however, $`Z_{\mathrm{bound}}`$ increases roughly linearly as a function of fragment (Z$`5`$) multiplicity with a slope varying from 3 to 4.5 when the energy increases from 25 to 50 AMeV. This again indicates that at higher energy the system is highly fragmented.

### 3.2 Charge distributions

The charge distributions, normalized to the number of single source events, can be compared all together in fig. 2. When the available energy increases, the emission of small fragments ($`Z<10`$) increases while that of fragments with $`Z>15`$ decreases producing steeper and steeper distributions; quite equivalent rates of production for fragments with $`10Z15`$ are observed between 32 and 50 AMeV. The domain of the heaviest fragments $`Z50`$, populated at the lowest energy, gradually vanishes at higher energies. This behaviour is quite remarkable, particularly seeing that the number of fragments is the same at all energies (fig. 1c). Roughly, one observes the following trend: the curves are gradually passing from a regime with two slopes, at the lowest incident energy, towards a regime with one slope, at the highest incident energy. The average measured charge of the fragments can be found in Table 1.

## 4 Charge correlation functions and enhanced production of nearly equal-sized fragment events

A few years ago a new method called higher order charge correlations Mor96  was proposed to enlighten any extra production of events with specific fragment partitions. The high sensitivity of the method makes it particularly appropriate to look for small numbers of events as those expected to have kept a memory of spinodal decomposition properties. Thus, such a charge correlation method allows to examine model independent signatures that would indicate a preferred decay into a number of equal-sized fragments in events from experimental data or from simulations. In this section and the following only events with 3 to 6 fragments (Z $``$ 5) at 32, 39, 45 and 50 AMeV will be considered. They represent in all cases 90-92 % of the selection previously considered. At these four energies, the detection efficiency is independent of the fragment charge, in the selected event samples. This is not so well verified at 25 AMeV.

### 4.1 Methods

The classical two fragment charge correlation method considers the coincidence yield $`Y(Z_1,Z_2)`$ of two fragments of atomic numbers $`Z_{1,2}`$, in the events of multiplicity $`M_f`$ of a sample. A background yield $`Y^{}(Z_1,Z_2)`$ is constructed by mixing, at random, fragments from different coincidence events selected by the same cut on $`M_f`$. The two particle correlation function is given by the ratio of these yields. When searching for enhanced production of events which break into equal-sized fragments, the higher order correlation method appears much more sensitive. All fragments of one event with fragment multiplicity $`M_f`$ = $`M=_Zn_Z`$, where $`n_Z`$ is the number of fragments with charge $`Z`$ in the partition, are taken into account. By means of the normalized first order:

$$Z=\frac{1}{M}\underset{Z}{}n_ZZ$$
(1)
and second order:

$$\sigma _Z^2=\frac{1}{M}\underset{Z}{}n_Z(ZZ)^2$$
(2)
moments of the fragment charge distribution in the event, one may define the higher order charge correlation function:

$$1+R(\sigma _Z,Z)=\frac{Y(\sigma _Z,Z)}{Y^{}(\sigma _Z,Z)}|_M$$
(3)
Here, the numerator $`Y(\sigma _Z,Z)`$ is the yield of events with given $`Z`$ and $`\sigma _Z`$ values. Because the measurement of the charge belonging to a given event is not subject to statistical fluctuations, we use here expression (2) rather than the “nonbiased estimator” of the variance, $`\frac{1}{M1}_Zn_Z(ZZ)^2`$, as proposed in Mor96  and used in our previous paper I31-Bor01 . Note that this choice has no qualitative influence on the forthcoming conclusions. The denominator $`Y^{}(\sigma _Z,Z)`$, which represents the uncorrelated yield of pseudo-events, was built in Mor96 , as for classical correlation methods, by taking fragments at random in different events of the selected sample of a certain fragment multiplicity. This Monte-Carlo generation of the denominator $`Y^{}(\sigma _Z,Z)`$ can be replaced by a fast algebraic calculation which is equivalent to the sampling of an infinite number of pseudo-events Des02 . Its contribution to the statistical error of the correlation function is thus eliminated. However, owing to the way the denominator was constructed, only the fragment charge distribution $`\mathrm{d}M/\mathrm{d}Z`$ of the parent sample is reproduced but the constraints imposed by charge conservation are not taken into account. This has, in particular, a strong effect on the charge bound in fragments $`\mathrm{d}M/\mathrm{d}Z_{\mathrm{bound}}`$ distribution. This fact makes the denominator yield distributions as a function of $`Z`$ wider and flatter than those of the numerator T27Tab00 . Consequently, even in the absence of a physical correlation signals, the ratio (3) is not a constant equal to one. The correlations induced by the finite size of the system (charge conservation) distorts the amplitude, or even may cancel other less trivial correlations. Therefore, a new method for the evaluation of the denominator Des02 , based on the “intrinsic probability” of emission of a given charge, was adopted. It minimizes these effects and replicates all features of the partitions of the numerator, except those (of interest) due to other reasons than charge conservation.

The goal of the method is to take into account in a combinatorial way the trivial correlations due to charge conservation. If there is no correlation between the charges, each charge can be fully described by an emission probability (referred to as intrinsic probability and noted $`{}_{}{}^{\mathrm{intr}}P_{Z}^{}`$). This fact was largely demonstrated for the bulk of multifragment emission around the Fermi energy Mor97 . Without charge conservation constraint, the intrinsic probabilities would be equal to the emission probability. Charge conservation makes the emission probability of high charge fragments smaller than their intrinsic probability, whereas it is the contrary for small charges. The probability to observe a given partition ($`𝐧:(n_1,\mathrm{},n_{Z_{\mathrm{tot}}})`$), at a given total multiplicity, $`M_{tot}=m=_Zn_Z`$, is obtained by the multinomial formula. If the total charge is fixed ($`Z_{\mathrm{tot}}=_ZZn_Z`$), the partition probabilities are given by:

$$P(𝐧|m)=\alpha m!\underset{Z}{}\frac{{}_{}{}^{\mathrm{intr}}P_{Z}^{n_Z}}{n_Z!}\delta _{Z_{\mathrm{tot}},_ZZn_Z},$$
(4)
where $`\alpha `$ is the normalization constant (so that $`_nP(𝐧|m`$) = 1) and $`\delta `$ is the Kronecker symbol. For the reason previously mentioned, no Z-dependent efficiency term is needed in this formula. These partition probabilities contain all information relative to the charges and their correlations. For example, the denominator $`Y^{}(\sigma _Z,Z)`$ is obtained by summing the probabilities of the partitions with given mean charge and standard deviation. Of course, the intrinsic probabilities are not direct experimental observables, they have to be evaluated by inversion of Eq. (4). However, this inversion is possible only if the physical correlations (not due to charge conservation) are weak. If the data sample contains only trivial correlations, then the higher order charge correlation function (or any other correlation function) is everywhere equal to 1. Local positive (negative) physical correlations will appear as peaks (holes). In the case of experimental events, it has to be noted that a set of intrinsic probabilities exists only if the selected sample comprises essentially single source events having reached thermal equilibrium (the source has to be described by a unique set of intrinsic probabilities). Otherwise the results would be equal to the convolution of the partition probabilities corresponding to the different source sizes, each of them being described by a different set of intrinsic probabilities. Hence, the convergence of the inversion procedure (Eq. 4) is a strong indication that these two conditions are fulfilled by the sample under study Des02 . The data samples correspond to multifragmenting single sources, supplemented by pre-equilibrium and secondary decay particles and light fragments. Hence, only fragments with charge greater than or equal to 5 have been retained for the analysis of the multifragmenting source. This feature has to be taken into account in the calculation of the denominator. The fragment partition ($`𝐍:(n_5,\mathrm{},n_{Z_{\mathrm{tot}}})`$) probability can be calculated as the sum over the complete partitions $`𝐧`$ which include $`𝐍`$. Noting $`M`$ the fragment multiplicity, it can be straightforwardly shown from Eq. (4) that the fragment partition probability reads:

$$P(𝐍)=f(M,Z_{\mathrm{bound}})P^{}(𝐍),$$
(5)
with

$$\begin{array}{c}\begin{array}{c}\hfill f(M,Z_{\mathrm{bound}})=\alpha \underset{𝐧:(n_1,\mathrm{},n_4)}{}\left(\genfrac{}{}{0pt}{}{m}{M}\right)(mM)!\\ \hfill \times \underset{Z=1}{\overset{4}{}}\frac{{}_{}{}^{\mathrm{intr}}P_{Z}^{n_Z}}{n_Z!}\delta _{Z_{\mathrm{tot}}Z_{\mathrm{bound}},_{Z=1}^4Zn_Z}\end{array}\hfill \\ P^{}(𝐍)=M!_{Z=5}^{Z_{\mathrm{tot}}}\frac{{}_{}{}^{\mathrm{intr}}P_{Z}^{n_Z}}{n_Z!},\hfill \end{array}$$
(6)
where $`\left(\genfrac{}{}{0pt}{}{m}{M}\right)`$ are the binomial coefficients. The value of $`Z_{\mathrm{tot}}`$ was fixed at the total charge of the system, $`Z_{\mathrm{projectile}}`$ \+ $`Z_{\mathrm{target}}`$; it has been noticed that, whereas the total charge conservation has to be explicitely included in the calculation, the results are only modified when $`Z_{\mathrm{tot}}`$ is lower than 85. This new method to build the denominator will be denoted as the Intrinsic Probability Method (IPM) in what follows. However, the explicit calculation of the intrinsic probabilities may not be the only method for building a denominator including only the correlations induced by charge conservation (another procedure is proposed in I41-Cha02 , see Appendix for a comparison with the IPM method). In any case, by definition, all valid methods lead to a partition probability described by Eq. (4).

A comparison of the results obtained in the framework of the higher order correlation method, with the analytical denominator, and with IPM, will be drawn in the following for an event sample resulting from a simulation of 32 AMeV 129Xe + 119Sn central collisions.

### 4.2 Stochastic mean-field and spinodal decomposition

Dynamical stochastic mean-field simulations have been proposed for a long time to describe processes involving instabilities like those leading to spinodal decomposition Ran90 ; Cho91 ; Bur91 . In this approach, spinodal decomposition of hot and dilute finite nuclear systems can be mimicked through the Brownian One-Body (BOB) dynamics Cho94 ; Gua97 ; I29-Fra01 , which consists in employing a Brownian force in the kinetic equations. Simulations have been performed for head-on 129Xe on 119Sn collisions at 32 AMeV. The ingredients of simulations can be found in I29-Fra01  as well as a detailed comparison between filtered simulated events (to account for the experimental device) and experimental data. A good agreement between both is revealed.

To refine the comparison higher-order charge correlations have been calculated for the simulated events, keeping the compact presentation proposed in I31-Bor01 : charge correlation functions are built for all events, whatever their multiplicity, by summing the correlated yields for all $`M`$ and by replacing the variable $`Z`$ by $`Z_{\mathrm{bound}}=M\times Z=_ZZn_Z`$. Uncorrelated events are constructed and weighted in proportion to real events of each multiplicity. This presentation is based on the experimental observation that the peaks observed independently for each fragment multiplicity correspond to the same $`Z_{\mathrm{bound}}`$ region I31-Bor01 . The variance bin was chosen equal to one charge unit. We recall that in the considered domain of excitation energy, around 3 MeV per nucleon I11-Mar98 ; I29-Fra01 , secondary evaporation leads to fragments one charge unit smaller, on average, than the primary $`Z1020`$ ones, with a standard deviation around one T27Tab00 . If a weak enhanced production of exactly equal-sized fragments exists, peaks are expected to appear in the interval $`\sigma _Z=01`$, because of secondary evaporation. This interval in $`\sigma _Z`$ is hence the minimum value which must be chosen to look for nearly equal-sized fragments. Any (unknown) intrinsic spread in the fragment size coming from the break-up process itself may enlarge the $`\sigma _Z`$ interval of interest. This will be discussed at the end of the section, for the moment we will consider only events with $`\sigma _Z<1`$, which corresponds to differences of at most two units between the fragment atomic numbers in one event.

Fig. 3 shows the correlation function calculated using the analytical denominator (a) or the denominator given by the IPM (b). Both functions are drawn versus the variables $`Z_{\mathrm{bound}}=M\times Z`$ and $`\sigma _Z`$. In fig. 3a, the equal-sized fragment correlations in the first bin are superimposed over trivial correlations due to the finite size of the system. For this reason, the ratio (3) is generally different from one and smoothly varies with the variables $`Z_{\mathrm{bound}}`$ and $`\sigma _Z`$. For each bin in $`Z_{\mathrm{bound}}`$ (fixed at 6 atomic number units), an exponential evolution of the correlation function is observed from $`\sigma _Z=78`$ down to $`\sigma _Z=23`$. This exponential evolution is thus taken as a “background” empirically extrapolated down to the first bin $`\sigma _Z=01`$ The amplitude of the correlation function in the domain $`Z_{\mathrm{bound}}=3660`$ is well above the background, with a confidence level higher than 90$`\%`$, proving thus a statistically significant enhancement of equal-sized fragment partitions. Of the 1% of events having $`\sigma _Z<1`$, (0.13$`\pm `$0.02)% (called extra-events from now on) are in excess of the background. In fig. 3b, as expected, all correlations due to the charge conservation are suppressed and the correlation function is equal to 1 (within statistical fluctuations) wherever no additional correlation is present. Again one observes peaks for $`\sigma _Z<1`$. The percentage of extra-events is 0.36$`\pm `$0.03%, higher than the one obtained with the previous method. Moreover, with this method, peaks also appear at the maximum values of $`\sigma _Z`$ for a given $`Z_{\mathrm{bound}}`$. They correspond to events composed of one big (a heavy residue) and several lighter fragments (sequentially emitted from the big one). In that case fusion-multifragmention does not occur and the peaks reveal the small proportion (0.15%) of events which undergo the fusion-evaporation process. The IPM approach may thus reveal other correlations not seen with the previous one. Therefore, we shall use this approach in the following analysis of the experimental data.

To conclude this part we can say that, although all events in the simulation arise from spinodal decomposition, only a very small fraction of the final partitions have nearly equal-sized fragments. Different effects: beating of different modes, coalescence of nascent fragments, secondary decay of the excited fragments and, above all, finite size effects are responsible for this fact Jac96 ; Colo97 . The signature of spinodal decomposition can only reveal itself as a “fossil” signal.

### 4.3 Experimental results

We shall present now higher-order charge correlations for the selected experimental events. This will be done for four incident energies (32, 39, 45 and 50 AMeV) in the framework of the IPM for the denominator. The first step consists in determining the intrinsic probabilities of fragments for each multiplicity and at each incident energy. These probabilities are obtained by a recursive procedure of minimization. The minimization criterion is the normalized $`\chi ^2`$ between the experimental fragment probabilities and the fragment partition probabilities given by (5). The calculated $`\chi ^2`$ were always lower than one: the lower is the incident energy and the larger is the multiplicity, the lower is the $`\chi ^2`$ value.

Charge distributions experimentally observed for the different fragment multiplicities are shown as an example in fig. 4; they correspond to the 39 AMeV incident energy. Dashed lines refer to the intrinsic probabilities calculated with IPM and the corresponding charge distributions are the full lines. The experimental charge distributions are faithfully described. Note that the ratio of the experimental to the calculated curves does not reveal any anomalous enhancement in a preferential domain of fragment atomic number. The charge distributions, summed over fragment multiplicities 3 to 6, are displayed in fig. 5, for the four incident energies. For each incident energy the intrinsic probabilities have been calculated independently for the different fragment multiplicities (see fig. 4) and weighted in proportion to real events of each fragment multiplicity. The slight differences for the intrinsic probabilities corresponding to the different multiplicities reveal small differences in the average excitation energies of the multifragmenting sources (on the other hand, it has been possible to fit all events from the BOB simulation, irrespective of their multiplicity, using the same intrinsic probability distribution due to the fact that the events correspond to the same initial conditions). We note again the excellent agreement between calculations and data at all incident energies. The contributions to the Z distribution of the three heaviest fragments of each partition are as well described as in fig. 4, and the charges bound in fragments (not shown) are also perfectly reproduced.

Fig. 6 illustrates the higher-order correlation functions measured for the different fragment multiplicities. It concerns single sources selected at 39 AMeV incident energy. To make the effects more visible, peaks with confidence level lower than 80% were flattened out. We observe significant peaks in the bin $`\sigma =01`$ for each fragment multiplicity. For M=6, peaks are essentially located in the bin $`\sigma _Z=12`$. As observed in simulations, peaks corresponding to events composed of a heavy residue and light fragments ($`\sigma _Z`$ in the region 5-10 associated with low $`Z`$) are also visible.

Let us now present the results for the different incident energies, summed over multiplicities 3-6. We recall again that IPM denominators of correlation functions are weighted in proportion to real events of each multiplicity. In fig. 7 the measured functions are displayed. If we exclude some peaks with low confidence levels (as done in fig. 8) correlation functions are equal or close to one except at low $`\sigma _Z`$ values. A summary of average-charge domains contributing to the correlation peaks in the first bin in $`\sigma _Z`$ is given in Table 2 as a function of fragment multiplicity. All multiplicities, associated to the largest $`Z`$ ranges (20-6 and 18-6) contribute to the peaks at 39 and 45 AMeV. At 32 AMeV the peaks only come from the low multiplicites ($`M`$ = 3-4), with a smaller $`Z`$ range : 21-11; at 50 AMeV, the situation is completely different and we only observe a contribution from $`M`$ = 6 with $`Z`$ = 7-9. Except in this latter case, the $`Z`$ domains for which correlation peaks are present are similar but slightly shifted towards lower values with increasing energy.

To be quantitative, correlation functions for the first bin in $`\sigma _Z`$ are displayed in fig. 9 with their statistical errors. Full points (open points) correspond to a confidence level higher (lower) than 90%. Two percentages of events relatively to the single source events ($`M_f=36`$) are presented in columns 2 and 3 of Table 3: they refer to the total number of events and of extra-events (taking into account correlations and anti-correlations) in the bin. The present analysis fully confirms at 32 AMeV incident energy the previous one I31-Bor01  and the extra-percentage of events with nearly equal-sized fragments is maximum at 39 AMeV.

A closer examination of fig. 3b and 8 reveals an islet of peaks, with a high confidence level, in the second and the third bins in $`\sigma _Z`$. They are located in the upper region of the $`M\times Z`$ ranges indicated in table 2. They correspond to events with a broader spread of charges. Observing such events in the simulation may indicate that the intrinsic spread (which is unknown) due to spinodal decomposition is larger than that coming from secondary evaporation, and hence that these peaks also sign the original process. The percentages of events with $`\sigma _Z<3`$ are also reported in table 3. The conclusions are the same as above: while more events have small values of $`\sigma _Z`$ when the incident energy increases, the percentage of extra-events shows a maximum at 39 AMeV, and vanishes at 50 AMeV.

## 5 Discussion

The first hint of a bulk effect for producing fragments in central collisions between heavy nuclei at a given moderate excitation energy ($`ϵ^{}7A`$MeV) was the measurement of identical fragment charge distributions for two different system sizes (32 AMeV 129Xe + natSn and 36 AMeV 155Gd + nat U) I12-Riv98 . However this feature could as well be interpreted from a statistical point of view, namely the dominance of phase space Mor97 ; Bon95 . Indeed both dynamical and statistical approaches were able to reproduce the experimental observation. With the same dynamical approach as that used in this paper, the experimental charge and multiplicities distributions were reproduced while the average fragment kinetic energies were accounted for within 20 % I29-Fra01 . These same properties were also well accounted for with the statistical model SMM Bon95 ; T16Sal97 ; LenBol00 ; TabBol00 . From this agreement between data and the two models we learnt that the dynamics involved is sufficiently chaotic to finally explore enough of the phase space and describe fragment production through a statistical approach. To go further more constrained observables were needed. Such was the goal of the present studies.

We have just seen that, at 32 AMeV incident energy, experimental correlation functions are similar to those obtained with events from dynamical simulations, BOB. Above all, in both of them, there is abnormal production of nearly equal-sized fragments pointed out by the peaks in the first $`\sigma _Z`$ bin(s). Supported by this simulation we thus attributed the greatest part of fusion-multifragmentation events to spinodal decomposition I31-Bor01 . These peaks hence appear as fossil fingerprints of the break-up process.

Correlation peaks are also observed for experimental data at higher incident energies: the excitation function is displayed in fig. 10. Information on the associated thermal excitation energies (and extra radial collective energy) involved over the incident energy domain studied can be provided by the SMM model which well describes static and dynamics observables of fragments. Starting from a freeze-out volume fixed at three times the normal volume, the thermal excitation energies of the dilute and homogeneous system, extracted from SMM, vary from 5.0 to 7.0 AMeV and the added radial expansion energy remains low: from 0.5 to 2.2 AMeV T16Sal97 ; T25NLN99 . A rise and fall of the percentage of “fossil partitions” from spinodal decomposition is measured. Fig. 10 reveals some difference between the experimental (full symbols) and simulated events (open symbols): the experimental percentages of extra events are closer to the simulated ones in fig. 10b than in fig. 10a. This means that the charge distributions inside an event are slightly narrower in the simulation than in the experiment either because of the primary intrinsic spread, or because the width due to evaporation is underestimated. For the considered system, incident energies around 35-40 AMeV appear as the most favourable to induce spinodal decomposition; it corresponds to about 5.5-6 AMeV thermal excitation energy associated to a very gentle expansion energy around 0.5-1 AMeV. The qualitative explanation for those numbers can be well understood in terms of a necessary compromise between two times. On one hand the fused systems have to stay in the spinodal region $``$ 100-150 fm/c Colo97 ; Idi94 ; Nor02 , to allow an important amplification of the initial fluctuations and thus permit spinodal decomposition; this requires a not too high incident energy, high enough however for multifragmentation to occur. On the other hand, for a finite system, Coulomb interaction and collective expansion push the “primitive” fragments apart and reduce the time of their mutual interation, which is efficient to keep a memory of “primitive” size properties. Above 45 AMeV incident energy part of trajectories followed by the system in the temperature-density plane may not sufficiently penetrate the spinodal region. Finally, considering secondary decay effects, they should be essentially similar between 39 and 50 AMeV because the primary fragment excitation energy per nucleon was found constant over this incident energy range I39-Hud02 . Particularly the excitation function shown in fig. 10 should not be affected by such effects.

Let us come now to the size of fragments associated to “fossil partitions”. From the theoretical point of view, spinodal instabilities have been mainly studied within semiclassical (as BOB) or hydrodynamical Nor02 ; Nor00 ; Hei88  frameworks. Precise information on the most unstable collective modes, which is needed for discussing fragment size produced by spinodal decomposition, can only be obtained when considering quantal effects. From the limited number of studies investigating such effects Ayi95 ; ChoHirs99 ; Col02  a few trends emerge: i) an increase of the number of unstable collective modes with the size of the system (up to multipolarity $`L`$=5-6 for A=140), ii) a decrease of the $`L_{max}`$ value when the temperature increases and iii) a dominance of the octupole mode ($`L`$=3) which appears as the most unstable. Charge correlations were only studied for events with 3 to 6 fragments, which precludes any information on the quadrupole mode ($`L`$=2). If we except results at the highest incident energy, which strongly differ from others, we do observe the dominance of events with three fragments, associated to the multipolarity $`L`$=3. Indeed the $`Z`$ domains marked in bold in table 2 correspond to the largest proportion of events which populate the first $`\sigma _Z`$ bin, which vary from 70% at 32 AMeV to 40% at 45 AMeV. Then the number of events with higher multiplicities, corresponding to larger $`L`$ values, progressively decreases, as expected.

At 50 AMeV results observed for the fragment sizes ($`\sigma _Z<1`$) are puzzling. Moreover the percentage of extra-events drops to zero when events with $`\sigma _Z<3`$ are considered. Further simulations and theoretical works are needed to progress in the interpretation of these data.

## 6 Conclusions

In conclusion, we have investigated charge correlation functions for compact and heavy fused systems which undergo multifragmentation, as a function of the incident energy, from 32 to 50 AMeV. At the lowest energy, we have confirmed an enhanced production of events with equal-sized fragments. Supported by theoretical simulations we have interpreted this enhancement as a signature of spinodal instabilities as the origin of multifragmentation of those systems in the Fermi energy domain. This fossil signal culminates for incident energy around 35-40 AMeV, which corresponds to the formation of a hot and dilute system at 0.3-0.4 the normal density and temperature around 4-5 MeV TGua96 ; I29-Fra01 . Spinodal decomposition describes the dynamics of a first order phase transition, and the present observations support the existence of such a transition for hot finite nuclear matter MDA00 ; LenBol00 ; MDA02 .

###### Acknowledgements.

The experimental results were obtained by the INDRA collaboration.

## Appendix

Very recently a new method for building charge correlation functions was proposed, of implementation easier than the IPM method: the denominator is built by mixing events through random exchanges of two fragments between two events under the constraint that the sum of the two exchanged fragments is conserved, besides the conservation of $`Z_{\mathrm{bound}}`$ I41-Cha02 . To test the sensitivity ot this method as compared to the IPM method, the present experimental data were also analyzed using the exchange method, by building correlation functions $`1+R(\sigma _Z,Z)`$ with $`\sigma _Z`$ being calculated from eq. 2. The results of the two methods are compared in fig. 11 and fig. 12, displaying for each bin of the plane ($`\sigma _Z`$, M $`\times Z`$), the deviation from 1, $`R(\sigma _Z,Z)`$, of the correlation function, normalized to its statistical error bar, $`\sigma _{1+R(\sigma _Z,Z)}`$, calculated from the numerator only.

The greater sensitivity of the IPM method is clearly seen on fig. 11: the values of the correlation function are closer to one ($`R`$=0) except at low $`\sigma _Z`$ where we observe correlations with a significant confidence level. Conversely, the exchange method, fig. 12, leads to a large dispersion of the values of $`R(\sigma _Z,Z)/\sigma _{1+R(\sigma _Z,Z)}`$, $``$1.6 times broader than with the IPM method at 32, 39 and 45 AMeV. At 50 AMeV both methods lead to similar dispersions.


format: text

text:
# 1 Introduction

## 1 Introduction

In 2-D integrable quantum field theories which can be considered as conformal field theories (CFTs) perturbed by a relevant operator, two-point correlation functions are complicated objects to study. However, using operator product expansion (OPE) in the short-distance limit one can reduce down their expression in terms of vacuum expectation values (VEVs) of local fields. Since four years, important progress has been made in this direction, as exact VEVs either of primary fields  or their first descendent  have been obtained explicitly. However it remains an open important problem to find all higher level VEVs of descendent fields and study their properties. Although a general method is still lacking, a case by case study based on CFT data provides a useful tool in order to determine some of the simplest higher level VEVs. Beyond the technical aspects, the knowledge of any of such quantities improves the analytical prediction for short-distance expansion of two-point functions, which can be better compared with the results obtained from the numerical study of the model (see  for instance).

Recently , we considered the Bullough-Dodd (BD) model and its quantum group restrictions, following the approach of  concerning the sinh-Gordon or sine-Gordon models. In Euclidean space, the action associated with the BD model writes

$`𝒜_{BD}={\displaystyle d^2x\left[\frac{1}{16\pi }(_\nu \phi )^2+\mu e^{b\phi }+\mu ^{}e^{\frac{b}{2}\phi }\right]}.`$ (1)
Here, the parameters $`\mu `$ and $`\mu ^{}`$ are introduced, as the two operators do not renormalize in the same way, on the contrary to any simply-laced affine Toda field theory. The purpose of this letter is to provide an exact expression for the VEV of the third level descendent fields (next to leading order in the UV limit of the two-point function) in the BD model, in order to complete the short-distance expansion of the two-point function calculated in . It should be stressed that differently from sine-Gordon (SG) model, this VEV is nonzero, due to the existence of a local conserved current of spin 3 in the BD model.

Finally, it is well-known that $`c<1`$ minimal CFT with action

$`𝒜`$ $`=`$ $`_{p/p^{}}+\lambda {\displaystyle d^2x\mathrm{\Phi }_{pert}}`$ (2)
perturbed by the operator $`\mathrm{\Phi }_{pert}\{\mathrm{\Phi }_{12},\mathrm{\Phi }_{21},\mathrm{\Phi }_{15}\}`$ can be obtained by a quantum group (QG) restriction of imaginary Bullough-Dodd model  with special values of the coupling. Here we denote respectively $`\mathrm{\Phi }_{12}`$, $`\mathrm{\Phi }_{21}`$ and $`\mathrm{\Phi }_{15}`$ as specific primary operators of the unperturbed minimal model $`_{p/p^{}}`$ and introduce the parameter $`\lambda `$ which characterizes the strength of the perturbation. Using this correspondence and the previous VEVs in the BD model, we will deduce $`0_s|L_3\overline{L}_3\mathrm{\Phi }_{lk}|0_s`$ in the perturbed minimal model (2).

## 2 VEVs of the third level descendent fields

The BD model can be regarded as a relevant perturbation of a Gaussian CFT in which case the field is normalized such that $`\phi _(z,\overline{z})\phi (0,0)_{Gauss}=2\mathrm{log}(z\overline{z})`$. For imaginary coupling $`b=i\beta `$, the perturbation is relevant for $`0<\beta ^2<1`$. Although the model (1) for real coupling is very different from the one with imaginary coupling in its physical content (this latter model contains solitons and breathers), there are good reasons to believe that the expectation values obtained in the real coupling case provide also the expectation values for the imaginary coupling. Then, let us now consider the two-point function in the BD model with imaginary coupling $`𝒢_{\alpha _1\alpha _2}(r)=e^{i\alpha _1\phi }(x)e^{i\alpha _2\phi }(y)_{BD}`$ with $`r=|xy|`$. It can be expanded in the short-distance limit ($`r0`$) which, as mentioned above, contains a term corresponding to the third descendent contribution. The result reads (see  for details)

$`𝒢_{\alpha _1\alpha _2}(r)=𝒢_{\alpha _1+\alpha _2}r^{4\alpha _1\alpha _2}\{1+_{1,2}(\alpha _1\beta ,\alpha _2\beta ,\beta ^2)\mu (\mu ^{})^2r^{63\beta ^2}+{\displaystyle \frac{(\alpha _1\alpha _2)^2}{4}}(\alpha _1+\alpha _2)r^4`$
$`{\displaystyle \frac{\alpha _1^2\alpha _2^2(\alpha _1\alpha _2)^2}{144}}𝒦(\alpha _1+\alpha _2)r^6+O\left(\mu ^2(\mu ^{})^4r^{126\beta ^2}\right)\}`$
$`+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\mu ^nr^{4\alpha _1\alpha _2+4n\beta (\alpha _1+\alpha _2)+2n(1\beta ^2)+2n^2\beta ^2}j_n(\alpha _1\beta ,\alpha _2\beta ,\beta ^2)𝒢_{\alpha _1+\alpha _2+n\beta }\left\{1+O\left(\mu (\mu ^{})^2r^{63\beta ^2}\right)\right\}`$
$`+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\mu _{}^{}{}_{}{}^{n}r^{4\alpha _1\alpha _22n\beta (\alpha _1+\alpha _2)+2n(1\frac{\beta ^2}{4})+\frac{n^2\beta ^2}{2}}j_n({\displaystyle \frac{\alpha _1\beta }{2}},{\displaystyle \frac{\alpha _2\beta }{2}},{\displaystyle \frac{\beta ^2}{4}})𝒢_{\alpha _1+\alpha _2\frac{n\beta }{2}}\left\{1+O\left(\mu (\mu ^{})^2r^{63\beta ^2}\right)\right\}`$
$`+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\mu ^n\mu ^{}r^{4\alpha _1\alpha _2+4(n\frac{1}{2})\beta (\alpha _1+\alpha _2)+2n(12\beta ^2)+2+2n^2\beta ^2}_{n,1}(\alpha _1\beta ,\alpha _2\beta ,\beta ^2)𝒢_{\alpha _1+\alpha _2+(n\frac{1}{2})\beta }`$
$`\times \left\{1+O\left(\mu (\mu ^{})^2r^{63\beta ^2}\right)\right\}`$ (3)
where we defined $`(\alpha )`$ and $`𝒦(\alpha )`$ by the ratios

$`(\alpha )={\displaystyle \frac{(\phi )^2(\overline{}\phi )^2e^{i\alpha \phi }_{BD}}{e^{i\alpha \phi }_{BD}}}\text{and}𝒦(\alpha )={\displaystyle \frac{(\phi )^3(\overline{}\phi )^3e^{i\alpha \phi }_{BD}}{e^{i\alpha \phi }_{BD}}}`$ (4)
and $`𝒢_\alpha =e^{i\alpha \phi }_{BD}`$ is the VEV of the exponential field in the BD model. A closed analytic expression for $`𝒢_\alpha `$ and $`(\alpha )`$ has been proposed in ref.  and ref. , respectively. Their expression involves an integral representation which is well defined if

$`{\displaystyle \frac{1}{2\beta }}<e(\alpha )<{\displaystyle \frac{1}{\beta }}`$ (5)
and obtained by analytic continuation outside this domain. Here we used the notations of  for the Dotsenko-Fateev integrals $`j_n(a,b,\rho )`$ and $`_{n,m}(a,b,\rho )`$. In particular, the integrals $`j_n(a,b,\rho )`$ have been evaluated explicitly in  with the result

$`j_n(a,b,\rho )=\pi ^n{\displaystyle \underset{k=0}{\overset{n1}{}}}{\displaystyle \frac{\gamma ((k+1)\rho )}{\gamma (\rho )}}\gamma (1+2a+k\rho )\gamma (1+2b+k\rho )\gamma (12a2b(n1+k)\rho )`$ (6)
where the notation $`\gamma (x)=\mathrm{\Gamma }(x)/\mathrm{\Gamma }(1x)`$ is used. Also, the integral $`_{1,1}(a,b,\rho )`$ can be obtained from the result of . Instead, the integral $`_{1,2}(a,b,\rho )`$ is a quite complicated object, and its explicit calculation goes beyond the purposes of this letter.

In the (Gaussian) free field theory, the composite fields $`(\phi )^3(\overline{}\phi )^3e^{i\alpha \phi }`$ are spinless with scale dimension

$`D\mathrm{\Delta }+\overline{\mathrm{\Delta }}=2\alpha ^2+6.`$ (7)
For $`0<\beta ^2<1`$ the perturbation is relevant and a finite number of lower scale dimension counterterms are sufficient to cancel the divergences arising in the VEVs of third level descendent fields. However, this procedure is regularization scheme dependent, i.e. one can always add finite counterterms. For generic values of $`\alpha `$ this ambiguity in the definition of the renormalized expression for these fields can be eliminated by fixing their scale dimensions to be (7). It exists however a set of values of $`\alpha `$ for which the ambiguity still remains. In the BD model with imaginary coupling, this situation arises if two fields, say $`𝒪_\alpha `$ and $`𝒪_\alpha ^{}`$, satisfy the resonance condition

$`D_\alpha =D_\alpha ^{}+2n(1\beta ^2)+2n^{}(1{\displaystyle \frac{\beta ^2}{4}})\text{with}(n,n^{})`$ (8)
associated with the ambiguity

$`𝒪_\alpha 𝒪_\alpha +\mu ^n\mu _{}^{}{}_{}{}^{n^{}}𝒪_\alpha ^{}.`$ (9)
In this case one says that the renormalized field $`𝒪_\alpha `$ has an $`(n|n^{})`$-th resonance  with the field $`𝒪_\alpha ^{}`$. Due to the condition (5) we find immediately that a resonance can appear between the third level descendent field $`(\phi )^3(\overline{}\phi )^3e^{i\alpha \phi }`$ and the following primary fields:

$`(i)`$ $`e^{i(\alpha \beta )\phi }\text{i.e.}(n|n^{})=(1|4)\text{for}\alpha ={\displaystyle \frac{1}{\beta }}{\displaystyle \frac{\beta }{2}};`$ (10)
$`(ii)`$ $`e^{i(\alpha +3\beta )\phi }\text{i.e.}(n|n^{})=(3|0)\text{for}\alpha =\beta ;`$
$`(iii)`$ $`e^{i(\alpha \frac{\beta }{2})\phi }\text{i.e.}(n|n^{})=(0|1)\text{for}\alpha ={\displaystyle \frac{2}{\beta }};`$
$`(iv)`$ $`e^{i(\alpha \frac{3\beta }{2})\phi }\text{i.e.}(n|n^{})=(0|3)\text{for}\alpha ={\displaystyle \frac{\beta }{2}}.`$
If we now look at the expression (3), we notice that the contribution brought by the third level descendent field in (4), and that of any of the exponential fields in $`(i)`$, $`(ii)`$, $`(iii)`$ and $`(iv)`$, have the same power behavior in $`r`$ ($`r^{4\alpha _1\alpha _2+6}`$) at short-distance for the corresponding values of $`\alpha `$. The integrals which appear in these contributions are, respectively:

$`(i)_{1,4}(\alpha _1\beta ,\alpha _2\beta ,\beta ^2),(ii)j_3(\alpha _1\beta ,\alpha _2\beta ,\beta ^2),`$
$`(iii)j_1({\displaystyle \frac{\alpha _1\beta }{2}},{\displaystyle \frac{\alpha _2\beta }{2}},{\displaystyle \frac{\beta ^2}{4}}),(iv)j_3({\displaystyle \frac{\alpha _1\beta }{2}},{\displaystyle \frac{\alpha _2\beta }{2}},{\displaystyle \frac{\beta ^2}{4}}).`$
As we will see, $`𝒦(\alpha )`$ (and similarly for the real coupling case) exhibits the same poles in order that the divergent contributions compensate each other. This last requirement leads for instance to a set of relations for $`𝒦(\alpha )`$. The third one reads

$`{\displaystyle \frac{\alpha _1^2\alpha _2^2(\alpha _1\alpha _2)^2}{144}}es_{\alpha =\frac{2}{\beta }}𝒦(\alpha )=\mu ^{}{\displaystyle \frac{𝒢_{\alpha \beta /2}}{𝒢_\alpha }}|_{\alpha =\frac{2}{\beta }}es_{\alpha =\frac{2}{\beta }}j_1({\displaystyle \frac{\alpha _1\beta }{2}},{\displaystyle \frac{\alpha _2\beta }{2}},{\displaystyle \frac{\beta ^2}{4}}),`$ (11)
which is used to fix the $`\alpha `$-independent part (normalization) of $`𝒦(\alpha )`$.

On the other hand, to determine the explicit form of the $`\alpha `$-dependent part of $`𝒦(\alpha )`$, we use the reflection relations method. Indeed, the BD model (1) can be regarded as two different perturbations of the Liouville field theory . First, one can consider the Liouville action where the perturbation is identified with $`e^{\frac{b}{2}\phi }`$. The holomorphic stress-energy tensor

$`T(z)={\displaystyle \frac{1}{4}}(\phi )^2+{\displaystyle \frac{Q}{2}}^2\phi `$ (12)
which ensures the local conformal invariance of the Liouville field theory with coupling $`b`$ can be written in terms of the standard Virasoro generators $`T(z)=_nL_nz^{n2}\text{and}\overline{T}(\overline{z})=_n\overline{L}_n\overline{z}^{n2}`$. Then, using the OPE of the stress-energy tensor of the Liouville part with any primary field, we have the relation

$`L_3\overline{L}_3e^{a\phi }=\left[\left({\displaystyle \frac{a+Q}{2}}\right)^2^3\phi {\displaystyle \frac{1}{2}}^2\phi \phi \right]\left[\left({\displaystyle \frac{a+Q}{2}}\right)^2\overline{}^3\phi {\displaystyle \frac{1}{2}}\overline{}^2\phi \overline{}\phi \right]e^{a\phi }.`$ (13)
Furthermore, taking the expectation value of the combination above and using the (Gaussian) equations of motion $`\overline{}\phi =0`$ we obtain

$`L_3\overline{L}_3e^{a\phi }_{BD}={\displaystyle \frac{a^2}{16}}(a+1/b)^2(a+b)^2(\phi )^3(\overline{}\phi )^3e^{a\phi }_{BD}.`$ (14)
Alternatively, we can consider $`e^{b\phi }`$ as a perturbation. Using both pictures and CPT framework, we deduce reflection relations between operators with the same quantum numbers. We report the reader to  for details about this approach. Consequently, if we denote

$`K(a)={\displaystyle \frac{(\phi )^3(\overline{}\phi )^3e^{a\phi }_{BD}}{e^{a\phi }_{BD}}},`$ (15)
then we obtain the following two functional relations

$`K(a)`$ $`=`$ $`\left[{\displaystyle \frac{(b+1/ba)(b+2/ba)(2b+1/ba)}{a(a+1/b)(a+b)}}\right]^2K(Qa),`$ (16)
$`K(a)`$ $`=`$ $`\left[{\displaystyle \frac{(b/2+2/ba)(b/2+4/ba)(b+2/ba)}{a(a+2/b)(a+b/2)}}\right]^2K(Q^{}+a).`$
Notice that these equations are invariant with respect to the symmetry $`b\frac{2}{b}`$ with $`aa`$ in agreement with the well-known self-duality of the BD-model. Assuming that $`K(a)`$ is a meromorphic function in $`a`$, we find that the “minimal” solution which follows from (11), (16) is:

$`K(a)`$ $`=`$ $`{\displaystyle \frac{1}{a^2}}\left[{\displaystyle \frac{m\mathrm{\Gamma }(\frac{b^2}{h})\mathrm{\Gamma }(\frac{2}{h})}{\mathrm{\Gamma }(\frac{1}{3})\sqrt{3}2(Q+Q^{})^2}}\right]^6\gamma \left({\displaystyle \frac{2ba+b^2+2}{h}}\right)\gamma \left({\displaystyle \frac{2ba2}{h}}\right)\gamma \left({\displaystyle \frac{2bab^2+4}{h}}\right)\gamma \left({\displaystyle \frac{2ba2b^2}{h}}\right)`$
$`\times \gamma \left({\displaystyle \frac{2ba+2b^22}{h}}\right)\gamma \left({\displaystyle \frac{2ba4}{h}}\right)\gamma \left({\displaystyle \frac{2ba+b^2+2}{h}}\right)\gamma \left({\displaystyle \frac{2bab^2}{h}}\right)`$
where $`h=6+3b^2`$ is the “deformed” Coxeter number . Here we have used the exact relation between the parameters $`\mu `$ and $`\mu ^{}`$ in the action (1) and the mass of the fundamental particle $`m`$  :

$`m={\displaystyle \frac{2\sqrt{3}\mathrm{\Gamma }(1/3)}{\mathrm{\Gamma }(1+b^2/h)\mathrm{\Gamma }(2/h)}}\left(\mu \pi \gamma (1+b^2)\right)^{1/h}\left(2\mu ^{}\pi \gamma (1+b^2/4)\right)^{2/h}.`$ (17)
Notice that $`K(a)`$ is invariant under the duality transformation $`b2/b`$ as expected, and contains all the expected poles. Accepting this conjecture and taking $`a=0`$, we obtain for instance:

$`L_3\overline{L}_3𝕀_{BD}={\displaystyle \frac{m^2}{2^{10/3}}}{\displaystyle \frac{\mathrm{\Gamma }^2(1+2/h)\mathrm{\Gamma }^2(1+b^2/h)\mathrm{\Gamma }^2(2/3)}{\gamma (1/2+2/h)\gamma (1/2+b^2/h)\gamma (1/3+6/h)\gamma (1/3+3b^2/h)}}f_{BD}^2`$ (18)
where $`f_{BD}`$ is the bulk free energy of the Bullough-Dodd model, obtained in .

## 3 Application to perturbed conformal field theories

For imaginary value of the coupling $`b=i\beta `$, with the substitutions $`\mu \mu `$ and $`\mu ^{}\mu ^{}`$ in (1) the BD model possesses quantum group symmetry $`U_q(A_2^{(2)})`$ with deformation parameter $`q=e^{i\pi /\beta ^2}`$ . At roots of unity, it is used to describe $`\mathrm{\Phi }_{12}`$, $`\mathrm{\Phi }_{21}`$ or $`\mathrm{\Phi }_{15}`$ perturbed CFTs (2). Let us consider the first case i.e. the $`\mathrm{\Phi }_{pert}\mathrm{\Phi }_{12}`$ perturbation, obtained for $`\beta ^2=p/p^{}`$ with $`1<p<p^{}`$ relative prime integers. In the following, $`\mathrm{\Phi }_{lk}`$ will denote a primary field of the minimal model $`_{p/p^{}}`$. The exact relation between the parameters $`\lambda `$ in (2) and the mass of the fundamental kink $`M`$ can be found in . Here we denote

$`\xi ={\displaystyle \frac{p}{p^{}p}}.`$ (19)
For unitary minimal models $`\xi >1`$ which, for $`m(\lambda )=0`$, corresponds to a massive phase . Using the particle-breather identification  $`m=2M\mathrm{sin}\left(\frac{\pi \xi }{3\xi +6}\right)`$ and parameter $`a=i\left(\frac{l1}{2\beta }\frac{k1}{2}\beta \right)`$ in $`K(a)`$ it is then straightforward to get the VEV:

$`{\displaystyle \frac{0_s|L_3\overline{L}_3\mathrm{\Phi }_{lk}|0_s}{0_s|\mathrm{\Phi }_{lk}|0_s}}`$ $`=`$ $`\left[{\displaystyle \frac{2^{2/3}\pi M\mathrm{\Gamma }(\frac{2+2\xi }{3\xi +6})}{\sqrt{3}\mathrm{\Gamma }(\frac{1}{3})\mathrm{\Gamma }(\frac{\xi }{3\xi +6})(1+\xi )}}\right]^6{\displaystyle \frac{1}{\xi ^2(1+\xi )^2(3\xi +6)^2}}`$ (20)
$`\times {\displaystyle \frac{\gamma (\frac{\eta 4\xi 3}{3\xi +6})\gamma (\frac{\eta 4\xi 3}{3\xi +6})\gamma (\frac{\eta +1+\xi }{3\xi +6})\gamma (\frac{\eta +1+\xi }{3\xi +6})}{\gamma (\frac{\eta +2\xi +3}{3\xi +6})\gamma (\frac{\eta +2\xi +3}{3\xi +6})\gamma (\frac{\eta 2\xi +1}{3\xi +6})\gamma (\frac{\eta 2\xi +1}{3\xi +6})}}.`$
Here $`|0_s`$ is one of the degenerate ground states of the QFT (2) (see  for a detailed discussion of the vacuum structure of the model).

For the second restriction $`\beta ^2=p^{}/p`$ which leads to the action (2) with $`\mathrm{\Phi }_{pert}\mathrm{\Phi }_{21}`$, the exact relation between the parameter $`\lambda `$ and the mass of the fundamental kink $`M`$ has been obtained in . The VEV of the third order descendent field immediately follows from (20) with the replacement $`\xi 1\xi `$.

Another subalgebra of $`U_q(A_2^{(2)})`$ is the subalgebra $`U_{q^4}(sl_2)`$. One can again restrict the phase space of the complex BD with respect to this subalgebra for a special value of the coupling $`\beta ^2=4p/p^{}`$ with $`2p<p^{}`$ relative prime integers in order to describe the third case, i.e. $`\mathrm{\Phi }_{pert}=\mathrm{\Phi }_{15}`$. The calculations are straightforward so we will not report them here.

## 4 Conclusion

In conclusion, we have conjectured an exact expression for the VEV of the third level descendent in the BD model. The highly nontrivial check of the residue conditions corresponding to the poles (10), (11) strongly supports our conjecture.

As explained above, the computation of the (UV behavior) of the two-point function involves an infinite tower of VEVs of descendent fields. It is not clear how to solve this problem in general. Even in the simplest case as the SG theory, a system of functional equations appear for the 4-th level descendent. A solution of this problem is still an open question.

Yet, for practical reasons, the computation of any higher order descendent VEV gives new information. Our results for example can be used in improving the comparison with the numerical computations in some interesting statistical models around their critical point: the critical Ising model in a magnetic field , the tricritical Ising model perturbed by its energy operator  or by its subleading magnetic operator,…

### Aknowledgements

We are grateful to Al. B. Zamolodchikov for discussions. MS’s work is supported by KISTEP exchange program 12-69-002. PB’s work is supported by JSPS fellowship.


text:
# DNA uptake into nuclei: Numerical and analytical results

## 1 Introduction

In eukaryotic cells molecular transport between the nucleus and the cytoplasm takes place through the nuclear pore complexes (NPCs) . Ions and small proteins can diffuse through these large channels, but large proteins and protein-RNA complexes are transported actively via a signal-mediated process. DNA, on the other hand, is rarely transported under normal circumstances. However, this latter process might be relevant during viral infection or genetic therapy. In a recent elegant experiment the uptake of single DNA molecules into the nucleus was studied by the group of M. Elbaum . In the present paper we provide a computational and theoretical study to understand and interpret the results of this experiment.

## 2 Description of the experiment

In the experiment double stranded bacteriophage $`\lambda `$ DNA was transported into cell-free nuclei of Xenopus laevis eggs. The DNA molecules were linked to polystyrene beads of $`1`$ or $`3\mu \text{m}`$. Beads lying near a nucleus were trapped by optical tweezers and pulled away. In some cases dragging the bead pulled the whole nucleus behind it at a distance (Fig. 1A). After the bead was released, it started a biased Brownian motion due to thermal fluctuations and its attachment to the DNA segment. Such stretches were repeated on the same bead at two minutes intervals, and the maximum extension decreased with each measurement (Fig. 1B). This suggests that some mechanism inside the nucleus pulled the DNA inward. The remaining extended length vs. time curve has two main features: (i) it starts linearly, and (ii) the length does not go to zero, instead, it converges to about $`1\mu \text{m}`$. The uptake itself is very slow (compared to the time scale of biochemical processes in living organisms), it takes several hundred seconds for a $`10\mu \text{m}`$-long DNA.

## 3 Theoretical estimates for the entropic forces

To explain the convergence of the remaining length vs. time curve and to give some estimations for the uptake force and the effective friction coefficient of the pore, Salman et. al  proposed a model in which the shrinking remaining segment of the DNA applies an increasing resistance against the uptake. They argue that because the entropic spring constant of a chain of length $`L`$ is proportional to $`L^1`$ and the average separation of the two ends of the remaining DNA segment (caused by volume exclusion by the bead and the nucleus) is proportional to $`L^{1/2}`$, the resulting entropic repulsion diverges as their product, $`L^{1/2}`$.

However, a detailed statistical analysis (described below) shows the opposite: the entropic repulsion imposed by the shrinking DNA segment decreases as $`L^1`$. This result makes it necessary to reconsider the proposed kinetics of the uptake of DNA by the nucleus.

It has long been known that a polymer of length $`L`$, anchored at one of its ends to a wall (with the other end being free), pulls on the wall by the force

$$F_{\text{free ended}}=F_{\text{free ended}}^{\mathrm{}}(1\gamma )\frac{k_\text{B}T}{L}$$
(1)
at the anchoring point, on average . $`T`$ denotes the absolute temperature, $`k_\text{B}`$ is the Boltzmann constant, and the length independent force term $`F_{\text{free ended}}^{\mathrm{}}>0`$ depends on the way of anchoring and on other properties of the polymer. $`\gamma 0.69`$ for self-avoiding and $`\gamma =1/2`$ for non-self-avoiding polymers. Due to force balance the rest of the polymer pushes the wall by the same amount of force.

Below we show that for non-self-avoiding polymers Eq. (1) can be obtained in a way much simpler than the original derivation, using a simple lattice chain model. Then we develop further this approach to take into account the role of the bead.

If the persistence length of the polymer is $`L_\text{p}`$, the lattice constant of the equivalent three-dimensional lattice chain is $`a=2L_\text{p}`$, which is called the Kuhn length . Let us suppose that the wall is located at the $`x=0`$ plane, and the anchoring point is at the origin $`A=(0,0,0)`$, as depicted in Fig 2. The constraint that the chain is not allowed to pierce through the wall is modeled by the condition that none of the $`x`$ coordinates of the chain points can be negative: $`0x_i`$, where $`0iN`$ and $`N=L/a`$ is the number of Kuhn segments. The partition function $`𝒵`$ of this system is simply the sum of all possible chain configurations that satisfy the above condition.

Let us choose an arbitrary point $`B`$ in the $`0x`$ half-space. According to the so-called reflection principle, each configuration that connects $`A`$ and $`B`$ but does not satisfy the constraint, has one and only one “mirror configuration” that connects $`A^{}=(2a,0,0)`$ and $`B`$. Illustrated in Fig. 2a, the mirror configuration is obtained by reflecting every chain point from $`A`$ and to the first unsatisfied point about the $`x=a`$ plane. Denoting the number of unconstrained configurations that connect points $`P`$ and $`Q`$ by $`𝒩(PQ)`$, the partition function can be expressed as

$$𝒵=\underset{x_N0,y_N,z_N}{}\left\{𝒩\left[A(x_N,y_N,z_N)\right]𝒩\left[A^{}(x_N,y_N,z_N)\right]\right\},$$
(2)
where each configuration that is counted mistakenly in the first term is canceled by its mirror configuration in the second term. The partition function can further be written as

$`𝒵`$ $`=`$ $`{\displaystyle \underset{x_N0,y_N,z_N}{}}\left\{𝒩\left[A(x_N,y_N,z_N)\right]𝒩\left[A(x_N+2a,y_N,z_N)\right]\right\}`$
$`=`$ $`{\displaystyle \underset{x_N\{0,a\},y_N,z_N}{}}𝒩\left[A(x_N,y_N,z_N)\right]`$
$``$ $`(2d)^N𝒫(0x_N<2a)(2d)^N\mathrm{\hspace{0.17em}2}a\rho (0),`$
where $`d=3`$ is the dimension of the space, $`𝒫(0x_N<2a)`$ is the probability that the chain ends between the planes $`x=0`$ and $`x=2a`$ in the continuum limit, and $`\rho (x)=(2\pi \sigma ^2)^{1/2}\mathrm{exp}[x^2/(2\sigma ^2)]`$ is the probability density of the $`x`$ coordinate of the end position with $`\sigma ^2=a^2N/d`$. The approximations are valid for large $`N`$.

Thus, the partition function can be expressed as

$$𝒵\pi ^{1/2}(2d)^N\left(\frac{2d}{N}\right)^{1/2},$$
(4)
the free energy is

$$=k_\text{B}T\mathrm{ln}(𝒵),$$
(5)
and the force acting on the anchored end of the chain is

$$F_{\text{free ended}}=\frac{}{L}=\frac{1}{a}\frac{}{N}k_\text{B}T\frac{\mathrm{ln}(2d)}{2L_\text{p}}\frac{1}{2}\frac{k_\text{B}T}{L},$$
(6)
in agreement with Eq. (1).

In the experiments, however, the distal end of the polymer is not free, but is bound to a polystyrene bead. Because the bead is large compared to the remaining DNA coil when the uptake stalls, we can approximate it with a mobile wall. The lattice chain calculation can be repeated even for this system. The partition function, i.e., the number of allowed configurations is now

$`𝒵`$ $`=`$ $`{\displaystyle \underset{x_N0,y_N,z_N}{}}\{𝒩[A(x_N,y_N,z_N)]𝒩[A^{}(x_N,y_N,z_N)].`$ (7)
$`.𝒩[A(x_N+2a,y_N,z_N)]+𝒩[A^{}(x_N+2a,y_N,z_N)]\},`$
where the first two terms are from Eq. (2), the third term – based on the reflection principle again (see Fig 2b) – cancels all configurations that go beyond the end point $`(x_N,y_N,z_N)`$, and the forth term compensates for those configurations that have been canceled twice (in both the second and the third terms). Replacing $`𝒩\left[A^{}(x_N,y_N,z_N)\right]`$ by $`𝒩\left[A(x_N+2a,y_N,z_N)\right]`$, and $`𝒩\left[A^{}(x_N+2a,y_N,z_N)\right]`$ by $`𝒩\left[A(x_N+4a,y_N,z_N)\right]`$, and making further simplifications one arrives at

$`𝒵`$ $`=`$ $`{\displaystyle \underset{x_N\{0,a\},y_N,z_N}{}}𝒩\left[A(x_N,y_N,z_N)\right]`$
$`{\displaystyle \underset{x_N\{2a,3a\},y_N,z_N}{}}𝒩\left[A(x_N,y_N,z_N)\right]`$
$``$ $`(2d)^N\mathrm{\hspace{0.17em}2}a[\rho (0)\rho (2a)]\pi ^{1/2}(2d)^N\left({\displaystyle \frac{2d}{N}}\right)^{3/2},`$
leading to the entropic force

$$F_{\text{mobile wall}}k_\text{B}T\frac{\mathrm{ln}(2d)}{2L_\text{p}}\frac{3}{2}\frac{k_\text{B}T}{L},$$
(9)
which indicates that a bead (or mobile wall) attached to the end of the DNA renders the entropic resistance even smaller. This can be understood intuitively by noticing that near the end of the uptake process the geometric constraints drastically reduce the available configuration space for the remaining segment, so it becomes more and more favorable for the DNA to slide through to the other side of the membrane.

Although the reflection principle can only be applied to non-self-avoiding polymers, we expect a qualitatively similar behavior for self-avoiding ones, too. In this latter case the reduction of the configuration space should also decrease the entropic resistance.

To verify our theoretical estimations we have performed numerical simulations. As Fig. 3 indicates, the numerical results for the free ended filament and for the one with the mobile wall are in good agreement with the relevant terms $`k_\text{B}T/(2L)`$ and $`3k_\text{B}T/(2L)`$ of the theoretical predictions, respectively. The case with the polystyrene bead, just as expected, is in between the above two extremes, closer to the one with the mobile wall.

## 4 Simulation technique

We modeled the double stranded DNA as a three-dimensional “bead-rod” chain , and used molecular dynamics simulations to follow its dynamics (see Fig. 4). Note that the “beads” of the bead-rod chain are artificial units and have nothing to do with the polystyrene bead in the real experiment. Monte Carlo methods were ruled out because they usually fail to reproduce the correct dynamics.

The overdamped equations of motion of the beads are

$$\mathrm{\Gamma }\dot{𝐑}_j=𝐅_j^{\text{constr}}+𝐅_j^{\text{bend}}+𝐅_j^{\text{metr}}+𝐅_j^{\text{ext}}+\xi _j,$$
(10)
where $`\mathrm{\Gamma }`$ is the friction coefficient and $`𝐑_j`$ is the position of the $`j`$th bead ($`1jN+1`$). The bond vectors can be defined as $`𝐫_j=𝐑_{j+1}𝐑_j`$ (for $`1jN`$). The forces on the r.h.s. of the equations are as follows:

* $`𝐅_j^{\text{constr}}`$ is the constraint force. These forces are responsible for keeping the length of the rods (or bonds) fixed: $`|𝐫_j|=b`$. The constraint forces can be decomposed into the bond tension forces: $`𝐅_j^{\text{constr}}=𝐅_j^{\text{tens}}𝐅_{j1}^{\text{tens}}`$ with $`𝐅_j^{\text{tens}}=T_j𝐫_j`$ and $`𝐅_0^{\text{tens}}=𝐅_{N+1}^{\text{tens}}=0`$ . The tensions $`T_j`$, represent $`N`$ new variables, however, we have exactly this number of additional equations expressing that the bond lengths are constant: $`𝐫_j\dot{𝐫}_j=0`$.
* $`𝐅_j^{\text{bend}}`$ is the bending force. The bending resistance is represented by angular springs between the adjacent bonds. The bending energy is $`E^{\text{bend}}=\frac{\kappa (b)}{2b}_{i=1}^{N1}\theta _i^2`$, where $`\kappa (b)`$ is the bond length dependent bending modulus, $`\theta _i=\mathrm{arccos}(𝐫_i𝐫_{i+1}/b^2)`$ is the included angle between bonds $`i`$ and $`i+1`$. The bending force is thus $`𝐅_j^{\text{bend}}=_jE^{\text{bend}}`$. In the literature usually the bending modulus of the continuous polymer, $`\kappa _0`$, is used to calculate the bending energy even in a discretized polymer model. However, we realized that for a discretized bead-rod polymer the bending modulus should depend on the bond length to give correct statistics. The detailed calculations will be presented elsewhere, we show here only the result: $`\widehat{\kappa }(\widehat{b})1\widehat{b}/6+\widehat{b}^2(4/\pi ^213/24)+\widehat{b}^3(1/62/\pi ^2)+\widehat{b}^4/96`$, where $`\widehat{b}=b/L_\text{p}`$ is the dimensionless bond length and $`\widehat{\kappa }(\widehat{b})=\kappa (b)/\kappa _0`$ is the dimensionless bending modulus. The persistence length $`L_\text{p}`$ of the polymer is related to the bending modulus through the expression $`L_\text{p}=\kappa _0/(k_\text{B}T)`$. A feature of this bond length dependent bending modulus is that $`\kappa (2L_\text{p})=0`$, which reflects the known fact that a continuous worm-like chain with persistence length $`L_\text{p}`$ can be well approximated by a freely jointed chain with bond length $`b=2L_\text{p}`$ on the large scale.
* $`𝐅_j^{\text{metr}}`$ is the metric or pseudo force, and is needed because the motion of a bead-rod chain in a thermal bath is different from the limiting case of very stiff bonds . The latter is a better model for real biopolymers, and the metric force turns the rigid rods into very stiff springs. $`𝐅_j^{\text{metr}}=k_\text{B}T_j\mathrm{ln}\sqrt{|g|}`$, where $`g`$ is the metric tensor of the hypersurface of the constrained $`N+1`$ beads in the originally $`3(N+1)`$ dimensional configuration space. The explicit form of this metric tensor for a linear bead-rod chain can be found in Refs. .
* $`𝐅_j^{\text{ext}}`$ is an optional external force acting on bead $`j`$. In the application presented in this paper, the polystyrene bead is attached to the end of the polymer by applying an attracting external force on the last polymer bead, while this bead and the wall of the nucleus exert repelling forces on the polymer. Depending on the situation, the bond length $`b`$ may be too large to represent the filament contour faithfully. For example, in the simulation the DNA molecule goes through a $`10\text{nm}`$ wide pore embedded in a $`100\text{nm}`$ thick membrane wall, while the bod length is $`100\text{nm}`$. The spatial resolution of the bead-rod chain in this case would be very low compared to the size of the pore. For this reason we introduced $`M`$ “virtual” beads along each rod, spaced uniformly. Their sole role is that external forces can act on them, and they appropriately redistribute these forces to the original beads at both ends of the bond they are sitting on. Hence their presence improve the spatial resolution, as far as external forces are concerned. We used $`M=10`$ in the simulations.
* $`\xi _j`$ represents spatially and temporally uncorrelated thermal noise, obeying the fluctuation-dissipation relation $`\xi _j(t)\xi _k(t^{})=2k_\text{B}T\mathrm{\Gamma }\delta _{jk}\delta (tt^{})`$.

In our model we neglected the excluded volume effect for the DNA filament itself, although we accounted for the exclusion of the filament by the polystyrene bead and the nuclear membrane wall. Hydrodynamic interactions were also neglected.

In our simulations the persistence length was $`L_\text{p}=50\text{nm}`$ and the friction coefficient of the beads was $`\mathrm{\Gamma }\pi \eta b`$, where the viscosity of the water is $`\eta =10^3\text{kg}\text{s}^1\text{m}^1`$ and the temperature is $`T=300\text{K}`$. We used the second-order Runge-Kutta (or midpoint) method to solve the system of first-order ordinary equations . The bond length $`b`$ varied between $`50\text{nm}`$ and $`100\text{nm}`$.

One can rely on the results of a dynamical simulation only if it is numerically stable. Therefore, it is important to be aware of the physical processes that determine the largest possible time step. In our case the diffusion of the polymer beads and the relaxation of the angular springs pose an upper limit for the time step. The typical distance a diffusive particle takes during time $`\mathrm{\Delta }t`$ is $`\mathrm{\Delta }x=\sqrt{2D\mathrm{\Delta }t}`$, where the diffusion constant is $`D=k_\text{B}T/\mathrm{\Gamma }`$. To ensure that the typical change of the bond length is smaller than a small fraction $`\epsilon _Db`$ of the bond length, the time step cannot be larger than

$$\mathrm{\Delta }t_D=\frac{\stackrel{~}{\mathrm{\Gamma }}\epsilon _D^2}{2k_\text{B}T}b^3,$$

where $`\stackrel{~}{\mathrm{\Gamma }}=\mathrm{\Gamma }/b`$ is the friction coefficient per unit length. Here the constraint forces were not taken into account. However, we expect that even with the constraint forces in effect the numerical error is still proportional to the typical diffusive displacement of the beads in one time step.

The relaxation time of the angular springs is

$$\tau _{AS}=\frac{\stackrel{~}{\mathrm{\Gamma }}}{6\kappa (b)}b^4.$$

Apart from the prefactor, this can be easily seen from dimensional analysis. Thus, the time step cannot be larger than $`\mathrm{\Delta }t_{AS}=\epsilon _{AS}\tau _{AS}`$ either, where $`\epsilon _{AS}1`$. In our simulations we used $`\epsilon _D=\epsilon _{AS}=0.05`$. The maximum allowable time step is then the minimum of $`\mathrm{\Delta }t_D`$ and $`\mathrm{\Delta }t_{AS}`$.

## 5 Discussion

A DNA filament threading through a nuclear pore is pulled by two opposing entropic forces from the two sides of the membrane. If the media on the two sides are identical, the length independent force terms cancel each other. Because the length dependent contributions (being proportional to $`k_\text{B}T/L`$) disappear quickly for large $`L`$, there must exist some active (energy consuming) mechanism responsible for the pulling of the DNA. Since the integral of the length dependent term (from $`L_\text{p}`$ to the total length of the DNA) is in the order of only a few $`k_\text{B}T`$, this force term is not even strong enough to trap the DNA in a configuration where almost the entire DNA is on one side of the membrane.

### 5.1 Chemical composition difference

The simplest way to implement an active pulling mechanism is to produce a difference between the length independent terms of the entropic forces on the two sides of the membrane. This can be achieved, e.g., via a difference in the ionic strengths or the pH values of the two media. Even the simple lattice chain model can account for this effect if the chemical composition of the medium affect the persistence length $`L_\text{p}`$ of the DNA \[see Eqs. (6) and (9)\].

### 5.2 Ratchet mechanism

Another possibility for the pulling is based on a ratchet mechanism . The ratchet mechanism as a means of active transport was first proposed by the group of G. Oster . The idea is that inside the nucleus some kind of molecule is able to bind to the DNA with a rate constant $`k_{\text{on}}`$ and to detach from it with a rate constant $`k_{\text{off}}`$. When such a molecule is bound to the DNA at some location, the DNA cannot slide back beyond this point.

Assuming that the possible binding sites for the molecule are separated by a distance $`b`$, we can estimate the effective pulling force produced by this ratchet mechanism. The grand canonical partition function per binding site is $`𝒵=1+e^{(E_b\mu )/(k_\text{B}T)}`$, where $`E_b`$ is the binding energy and $`\mu `$ is the chemical potential for the molecule. The probability that the site is unoccupied is $`p_{\text{off}}=1/𝒵`$, and the probability that it is occupied is $`p_{\text{on}}=e^{(E_b\mu )/(k_\text{B}T)}/𝒵`$. Thus, the partition function can be written as $`𝒵=1+p_{\text{on}}/p_{\text{off}}=1+k_{\text{on}}/k_{\text{off}}`$. The free energies of a binding site inside and outside the nucleus are $`_{\text{in}}=k_\text{B}T\mathrm{ln}𝒵`$ and $`_{\text{out}}=k_\text{B}T\mathrm{ln}1=0`$, respectively. Finally, an effective pulling force can be defined as

$$F_{\text{ratchet}}=(_{\text{in}}_{\text{out}})/b=\mathrm{ln}\left(1+k_{\text{on}}/k_{\text{off}}\right)k_\text{B}T/b.$$
(11)
Note that this is the force necessary to stop the pulling. The pulling speed is proportional to this effective force only if the equilibration time of the binding of the molecules $`(k_{\text{on}}+k_{\text{off}})^1`$ is shorter than the diffusion time of the DNA over the distance $`b`$. Otherwise the pulling speed is smaller.

### 5.3 Friction in the pore

The friction of the DNA inside the nuclear pores can significantly slow down the uptake process. Let us just estimate how long it would take for an entire DNA to diffuse through the pore if there was no friction inside. The diffusion can be envisioned as the subsequent diffusion of $`N=L_{\text{DNA}}/a`$ Kuhn segments (where $`a=2L_\text{p}`$ is the Kuhn length). The diffusion time of a Kuhn segment can be approximated as $`\tau _\text{K}=(1/2)a^2(\pi \eta a)/(k_\text{B}T)`$, thus the diffusion time of the entire DNA is about $`\tau =N^2\tau _\text{K}=L_{\text{DNA}}^2L_\text{p}(\pi \eta )/(k_\text{B}T)`$. As expected for a polymer without self-exclusion, the diffusion time is proportional to the square of its length .

The diffusion time $`\tau `$ for a $`15\mu \text{m}`$ long DNA is in the order of 10 sec. Because the uptake occurs much slower, in the order of 1000 sec, the friction in the pore must be at least 100-fold greater than the hydrodynamic friction: $`\mathrm{\Gamma }_{\text{pore}}>100L_\text{p}\pi \eta 10^8`$kg/s. Such a strong friction can be explained by weak binding of the DNA to the nuclear pore .

### 5.4 Uptake kinetics I: weak pulling

The experiments show that the uptake slows down in the end and the remaining length of the DNA converges to $`L\text{conv}1\mu \text{m}`$. A constant pulling force $`F`$ (produced by either a chemical composition difference or a ratchet mechanism) corresponds to a linear potential $`U(L)=FL`$ with a reflecting boundary at $`L=0`$. In such a potential the expectation value of $`L`$ at equilibrium is given by $`k_\text{B}T/F`$. Thus, the measured equilibrium distance $`L\text{conv}`$ is compatible with an $`F=k_\text{B}T/L\text{conv}4`$fN pulling force. From the average pulling speed, $`v2\times 10^8`$m/s, one can also calculate the friction coefficient: $`\mathrm{\Gamma }_{\text{pore}}=F/v2\times 10^7`$kg/s, which satisfies the condition given in the preceding paragraph.

The pulling force and the friction coefficient would be even smaller if the length dependent term of the entropic force was also taken into account in the calculation of the equilibrium distance.

We studied the uptake process by the ratchet mechanism using numerical simulations. We assumed that the molecules could bind only to the beads of the bead-rod chain. The simulations started from a configuration in which one end of the DNA had already got through the pore, and the first bead was not allowed to slip back through the pore. To mimic the experiment, a large bead of radius $`R_\text{b}`$ was attached to the other, outer end of the DNA. The diameter of the pore was $`10\text{nm}`$, the width of the wall was $`100\text{nm}`$, and the friction at the pore was neglected. Fig. 5 shows a sequence of such a simulation.

Identifying the extended length of the remaining DNA by its arc length, we plotted the extended length vs. time curve, averaged over 20 runs. It resembles the experimental curve: it starts linearly and converges to a constant value, about $`1\mu \text{m}`$, which coincides with the experimental one. The time scale of the simulation is shorter than that of the experiments, because we omitted the friction inside the pore, in order to save computational time.

The negative side of this interpretation is that the uptake is stochastic. Individual runs show back-slips at both the linear and the horizontal part of the curve. The stationary state of the system is not an energy minimum but a dynamic equilibrium between the fluctuations of the DNA filament and the weak pulling of the transport mechanism. Since such fluctuations are not apparent in the experiments, we examine another possibility for the kinetics of the uptake.

### 5.5 Uptake kinetics II: strong but slow ratcheting

A monotonous uptake can be achieved either by a strong continuous pulling in conjunction with a strong friction or by a strong but slow ratchet mechanism. In such a ratchet mechanism the speed of the uptake would be set by the small binding rate $`k_{\text{on}}`$, however, the much smaller unbinding rate $`k_{\text{off}}k_{\text{on}}`$ would ensure that unbinding is negligibly rare and, thus, the ratcheting is practically irreversible. Note that the large $`k_{\text{on}}/k_{\text{off}}`$ ratio corresponds to a large effective pulling force.

Due to the negligibly rare unbinding events the ratchet mechanism is also supported by the fact that the uptake can be stopped but never reversed in the experiments  by applying pN forces on the polystyrene bead. A continuous pulling would be incompatible with this observation.

The only phenomenon that cannot be explained by such a ratchet mechanism is the convergence of the remaining length to a non-zero value. There are, however, alternative explanations for this:

* the non-zero distance might be an artifact of the image processing;
* there might be unknown structures on the nuclear envelope that prevent the polystyrene bead from getting in contact with the nucleus;
* the pN forces used to extend the DNA might be sufficiently strong to deform the membrane by hundreds of nm-s at the vicinity of the pore ;

or a combination of these.

## 6 Conclusions

Our numerical and analytical studies of the uptake of DNA filaments into the nucleus suggest two possible mechanisms for the uptake process. None of them can be ruled out on the basis of existing experimental data. One of the mechanisms explains the convergence of the remaining length of the DNA to a finite value, but has an intrinsically stochastic nature (back-slips occur). This mechanism alone cannot account for the stopping of the uptake under retarding forces in the pN range. The other mechanism seems more plausible, it predicts a monotonous uptake, but requires additional elements to explain the convergence. Only further experiments can unambiguously distinguish between these mechanisms, focusing on their differences.

We would like to acknowledge Collegium Budapest for its hospitality. We also thank M. Elbaum for many useful discussions and helpful comments. This research has been supported by the Hungarian National Research Grant Foundation (OTKA T033104).

## References