text:
# What does the potential energy landscape tell us about the dynamics of supercooled liquids and glasses?

## Abstract

For a model glass-former we demonstrate via computer simulations how macroscopic dynamic quantities can be inferred from a PEL analysis. The essential step is to consider whole superstructures of many PEL minima, called metabasins, rather than single minima. We show that two types of metabasins exist: some allowing for quasi-free motion on the PEL (liquid-like), the others acting as traps (solid-like). The activated, multi-step escapes from the latter metabasins are found to dictate the slowing down of dynamics upon cooling over a much broader temperature range than is currently assumed.

The theoretical understanding of supercooled liquids and glasses is one of the important unsolved problems of statistical physics. Since the pioneering work of Goldstein Goldstein (1969), it has been realized that the potential energy landscape (PEL) viewpoint is useful for characterizing supercooled liquids and glasses. To this end, one considers the high-dimensional vector of all particle coordinates as a point moving on the surface of the total potential energy Stillinger (1995); Debenedetti and Stillinger (2001); Sastry et al. (1998). At sufficiently low temperatures, the system resides near the local minima of the high-dimensional PEL. It has turned out that the PEL description of thermodynamics starts to work when cooling below approximately $`T=2T_\text{c}`$, where $`T_\text{c}`$ is the critical temperature of mode-coupling theory Götze and Sjogren (1992). In this temperature regime the statistical properties of PEL minima fully determine all thermodynamic properties Sastry (2001); Sciortino et al. (2000, 1999); Büchner and Heuer (2000); La Nave et al. (2002a).

What is still lacking in the literature is a full quantitative understanding of dynamics, i.e. of the slowing down of molecular motion upon cooling. In this connection, the single-particle diffusion coefficient, $`D(T)`$, is a typical macroscopic quantity of interest. In general, knowledge of thermodynamic properties is not sufficient to predict $`D(T)`$. However, the validity of the Adam-Gibbs relation Adam and Gibbs (1965), connecting $`D(T)`$ to the configurational entropy, suggests the presence of a strong link between thermodynamics and dynamics Saika-Voivod et al. (2001); Scala et al. (2000); Starr et al. (2001). For a model glass-former we have demonstrated via computer simulations that $`D(T)`$ can indeed be inferred from a PEL analysis in a quantitative way. The essential step, however, is to consider whole superstructures of many PEL minima, called metabasins Stillinger (1995); Büchner and Heuer (2000); Doliwa and Heuer (2003a); Denny et al. (2003); Doliwa and Heuer (2003b); Middleton and Wales (2001), rather than single minima. Metabasins are reminiscent of protein-folding funnels Bryngelson et al. (1995) or related structures in small clusters Ball et al. (1996).

The goal of this paper is to derive physical consequences of this mapping between PEL and dynamics and thus to obtain a coherent picture of the glass transition. In particular we will dwell on the interpretation of the activation energies for fragile systems in the Arrhenius diagram, the concept of liquid-like and solid-like behavior Fischer et al. (2002), and the question of a possible crossover at $`T_c`$, as implied by the mode-coupling temperature and the analogy to $`p`$-spin models.

In what follows we analyze a binary mixture of Lennard-Jones particles (BMLJ). All details of the BMLJ have been reported in Doliwa and Heuer (2003a). This system is designed to model a Ni80P20 mixture. Langevin molecular dynamics have been employed with a simulation box of $`N=65`$ particles at density $`\rho =1.2\sigma _{\text{ref}}^3`$. Throughout the paper, all quantities are given in reduced units, which, in the case of Ni80P20, correspond to $`\sigma _{\text{ref}}=2.218`$ Å, $`ϵ_{\text{ref}}=7.762`$ kJ/mol, $`t_{\text{ref}}=1.323`$ fs, $`T_{\text{ref}}=934`$ K, and $`D_{\text{ref}}=\sigma _{\text{ref}}^2/t_{\text{ref}}=0.372\text{cm}^2/s`$. For the temperature regime of our analysis, a system size of $`N=65`$ is provably sufficient to reproduce the dynamic properties of a macroscopic system Büchner and Heuer (1999); Doliwa and Heuer (2003c): We have checked, e.g., that $`D(T)`$ is the same for $`N=130`$ within 10% for $`T>T_\text{c}=0.45T_{\text{ref}}`$. The minima are obtained by a steepest descent minimization. As usual we perform regular quenches of the system to monitor the energies of the corresponding energy minima. In addition, we use temporal interval bisectioning to resolve the elementary transitions between minima Doliwa and Heuer (2003b).

The definition of metabasins (MBs) is motivated by the observation that during the molecular dynamics in configuration space the system performs several back-and-forth jumps between adjacent minima until finally this region is left. The description of the transport may be simplified if these minima are regarded as a single superstructure, i.e. a MB Doliwa and Heuer (2003b). Thus, the MBs correspond to an appropriate tiling of configuration space. Although irrelevant for the understanding of thermodynamic properties, they are of outstanding importance for dynamics. As shown in Saksaengwijit et al. (2003) it is possible to find a strict definition of MBs which can be also used in practice. Similarly to the previous consideration of energy minima the time evolution of the system may be regarded as a continuous sequence of MB visits with individual residence (waiting) times $`\tau `$.

We briefly summarize the main results of our analysis which will be important throughout the paper. They all hold for $`T<2T_c`$. (i) Analyzing the escape characteristics from MBs by repeated simulations from the same MBs at different temperatures $`T`$ Doliwa and Heuer (2003a), it turns out that the mean residence time in MBs of energy $`ϵ`$ is given to a good approximation by

$$\tau (ϵ,T)\tau _0(ϵ)\mathrm{exp}(E(ϵ)/k_\text{B}T).$$
(1)
The Arrhenius form of $`\tau (ϵ,T)`$ suggests that the escapes from MBs are thermally activated, with a barrier height of $`E(ϵ)`$, introduced here as a fitting parameter. The prefactor $`\tau _0(ϵ)`$ turns out to be basically independent of $`ϵ`$. (We define $`ϵ`$ as the energy of the lowest minimum within a MB.) (ii) $`E(ϵ)`$ is directly related to the PEL barriers in the high-dimensional surrounding of the MBs. By a detailed discussion of this issue, incorporating the funnel-like nature of MBs, this was verified Doliwa and Heuer (2003a). Thus, by analyzing the local topology of MBs, $`E(ϵ)`$ can be predicted. We may interpret $`E(ϵ)`$ as the depth of a MB of energy $`ϵ`$. The deeper the MB in the PEL, i.e. the lower $`ϵ`$, the higher the activation energy. In this analysis one has to take into account that the escape from a MB is a multi-step process, i.e. generally comprises hops between several minima. (iii) A crucial quantity is the average residence time $`\tau (T)`$. It is defined as the average over all MBs encountered at a specific temperature, i.e. $`\tau (T)=\text{d}\tau \tau \phi (\tau ,T)`$. By $`\phi (\tau ,T)`$ and $`\phi (ϵ,T)`$, respectively, we denote the distribution of waiting times and energies of visited MBs. (iv) The population of MBs with energy $`ϵ`$ is given by $`p(ϵ,T)=\tau (ϵ,T)\phi (ϵ,T)/\tau (T)`$. Numerically, it is indistinguishable from the population of minima, normally studied in this field and is purely gaussian in the accessible energy range Doliwa (2002). Using Eq.1 the inverse average waiting time can be written as

$`\tau _0/\tau (T){\displaystyle \text{d}ϵp(ϵ,T)\mathrm{exp}(E(ϵ)/k_\text{B}T)}`$
$`\mathrm{exp}(E_{\text{eff}}(T)/k_\text{B}T),`$ (2)
where the energy dependence of $`\tau _0(ϵ)`$ has been neglected. $`E_{\text{eff}}(T)`$ can be interpreted as the typical barrier height encountered at temperature $`T`$. (v) A simple relation exists between $`\tau (T)`$ and $`D(T)`$ via

$$D(T)\frac{a^2}{6N\tau (T)},$$
(3)
with a temperature-independent effective jump length $`a1.0\sigma _{\text{ref}}`$. Thus, the temperature dependence of $`D(T)`$ is exclusively determined by the average waiting time, whereas all spatial aspects of diffusion are temperature independent. Such a simple relation does not hold on the level of single minima. Actually, as a side effect it turns out that the dynamics can be basically described as a random-walk between MBs, thus suggesting a non-topographic view. With the definition $`D_0a^2/6N\tau _0`$ we can thus write

$$D(T)D_0\mathrm{exp}(E_{\text{eff}}(T)/k_\text{B}T).$$
(4)
Note that the ingredients of the constant $`D_0`$, i.e. $`a`$ and $`\tau _0`$, have been obtained from the simulations mentioned above. In the units of our simulation we have $`D_0a^2/6N\tau _01.310^5D_{\text{ref}}`$. If the BMLJ is mapped on a Ni80P20 alloy this corresponds to $`D_04.8\times 10^6\text{cm}^2/\text{s}`$.

Here we discuss the relevant physical implications of Eqs.1-4. In Fig.1(a) we show the barrier height $`E(ϵ)`$ as a function of MB energy, and the population of MB energies, $`p(ϵ,T)`$, for different temperatures. From these quantities we obtain an estimate for the diffusion coefficient, denoted $`D_{\text{est}}(T)`$. A comparison of the estimated with the true diffusion coefficient $`D(T)`$ (directly from molecular dynamics) is shown in (b). The very good agreement at $`T_\text{c}T2T_\text{c}`$ reflects the consistency of our approach and shows that we can indeed understand the macroscopic dynamics from knowledge of the thermodynamics and the local barriers.

From the simulations, we know $`p(ϵ,T)`$ and $`E(ϵ)`$ for $`ϵ>302ϵ_{\text{ref}}`$. Both quantities do not display finite-size effects in this energy range Doliwa and Heuer (2003c). Due to the Boltzmann weighting in Eq.What does the potential energy landscape tell us about the dynamics of supercooled liquids and glasses? one may hope that the available information is sufficient to reliably estimate $`D_{\text{est}}(T)`$ also for $`T<T_c`$. To check this hypothesis, we have calculated $`D_{\text{est}}(T)`$ thereby using extreme candidates for possible $`ϵ`$-extrapolations ($`ϵ<302ϵ_{\text{ref}}`$); see Fig.1. It turns out that $`D_{\text{est}}(T)`$ is very insensitive to $`E(ϵ)`$. For $`p(ϵ;T)`$ we have extended the gaussian description but have set $`p(ϵ;T)=0`$ below different cutoff-energies $`ϵ_{\text{min}}`$. (A lower bound for $`ϵ_{\text{min}}`$ is given by the condition that the number density of minima is of the order of $`1/ϵ_{\text{ref}}`$ Sastry (2001). This leads to $`306ϵ_{\text{ref}}ϵ_{\text{min}}302ϵ_{\text{ref}}`$ Doliwa (2002).) As shown in Fig.1 (b), the resulting uncertainty in the estimation of $`D_{\text{est}}(T)`$ is small up to $`T>0.35T_{\text{ref}}`$ and becomes significant for lower temperatures. Since a cutoff at $`302ϵ_{\text{ref}}`$ is probably far too drastic (see Doliwa (2002) for strong indications), however, the slowing down can be thus predicted semi-quantitatively down to $`TT_g`$, a region which very likely will never be accessible by molecular dynamics simulations.

Moreover, Eq.4 implies that the activation energy which is obtained from connecting $`(1/T=0/T_{\text{ref}},D=D_0)`$ and $`(1/T,D(T))`$ in an Arrhenius plot by a straight line has a simple interpretation: it is the typical barrier height the system experiences at a given temperature Dyre (1995). This non-trivial interpretation has very recently been confirmed by hyperquench experiments Angell et al. (2003), where the typical barrier height was measured by probing the specific heat during the reheating process.

To obtain a more detailed picture of MBs we have determined the distribution of waiting times, $`\phi (\tau ,T)`$, as shown in Fig.2(a). By visual inspection, one can distinguish two different time regimes (separated by the vertical line at $`\tau ^{}5000t_{\text{ref}}`$). In terms of energetics one expects that the long $`\tau `$’s arise from deep traps where the system is caught for a long time. In the opposite limit one may imagine that there are quite shallow MBs which do not strongly confine the system so that it will mainly stay close to the high-dimensional boundary of these MBs (see the sketches in (a)). Thus, short waiting times correspond to just scratching the MBs. This expectation can be verified by analyzing the trajectory of the system during the MB visits; an example of such a computation is given in Fig.1(b). There we plot $`\mathrm{cos}\alpha (\tau )`$ over MBs of lifetime $`\tau `$, where $`\alpha `$ is the angle between the entry point, the lowest minimum of the MB and the exit point of the trajectory (see the sketches in (b)). Again, two time regimes with a temperature-independent crossover time $`\tau ^{}`$ can be identified. Short visits to MBs lead to small values of $`\alpha `$, meaning that the system indeed merely scratches these MBs. For $`\tau >\tau ^{}`$ the value of $`\mathrm{cos}\alpha (\tau )`$ reaches a limiting value of ca. 0.2 , indicating that entry and exit points are largely uncorrelated. This should be the case after a long equilibration inside a MB with many possible exits. Due to the difference in stability, we call the MBs with $`\tau >\tau ^{}`$ solid-like, the other MBs liquid-like. This notation has been borrowed from two-state models where these two types of configurations have been postulated Fischer et al. (2002).

Clearly, the molecular slowing down upon cooling is caused by the enhancement of solid-like configurations. This can be quantified in two different ways. Firstly, one may determine the fraction of solid-like configurations the system encounters, i.e. $`\phi _{\text{sol}}(T)_\tau ^{}^{\mathrm{}}𝑑\tau \phi (\tau ,T)`$, where $`\phi _{\text{sol}}<0.5`$ implies that more liquid-like than solid-like configurations are visited. Secondly, we can specify the fraction of time spent in solid-like configurations, which can be expressed as $`p_{\text{sol}}(T)_\tau ^{}^{\mathrm{}}\text{d}\tau p(\tau ,T)`$ with $`p(\tau ,T)=\tau \phi (\tau ,T)/\tau (T)`$. $`p_{\text{sol}}>0.5`$ implies that the system is mostly residing in solid-like configurations. In Fig.3 we show the temperature dependence of $`\phi _{\text{sol}}`$ and $`p_{\text{sol}}`$. Three different temperature regimes can be distinguished (see the sketch in Fig.3): For $`T>2T_\text{c}`$ both quantities are smaller than 0.5. Thus, the system behaves liquid-like. Interestingly, this temperature regime (defined by dynamics) is exactly the temperature regime for which the minima no longer influence the thermodynamic properties of the system Sastry et al. (1998). Below $`2T_\text{c}`$, $`p_{\text{sol}}`$ is larger than 0.5. Thus, the system mainly resides in solid-like MBs. Finally, below a temperature near $`T_\text{c}`$, also $`\phi _{\text{sol}}`$ exceeds 0.5, i.e. we have a trap-to-trap motion. This crossover, however, is very gradual.

With Eq.What does the potential energy landscape tell us about the dynamics of supercooled liquids and glasses? we can also analyze the question in which temperature regime the dynamics are dominated by activated processes. We call a process activated if the activation energy $`E(ϵ)`$ is larger than $`5k_\text{B}T`$. To this end, we consider the fraction of $`\tau (T)`$ which is made up by activated hops, i.e.

$$\tau (T)^1_{E(ϵ)>5k_\text{B}T}\text{d}ϵ\tau (ϵ,T)\phi (ϵ,T)p_{\text{act}}(T),$$

as shown in Fig.3. The crossover temperature for which $`p_{\text{act}}=0.5`$ is close to $`1.5T_\text{c}`$. Thus, already significantly above $`T_\text{c}`$, the dynamics are dominated by activated processes. Please note that the MBs contributing to the above integral are solid-like since they fulfill $`\tau (ϵ,T)\tau ^{}`$ (see Doliwa and Heuer (2003a)). This is why $`p_{\text{act}}(T)p_{\text{sol}}(T)`$.

The present results seem to be at variance with the current understanding that activated transitions between PEL minima set in at $`T_\text{c}`$, whereas above $`T_\text{c}`$ they are irrelevant. The latter scenario is motivated by the properties of $`p`$-spin models Crisanti and Ritort (2000) and has been backed by the observations that the number of free directions (obtained via instantaneous-normal-mode analyses La Nave et al. (2002b); Chowdhary and Keyes (2002)), is directly related to the diffusion coefficient and that above $`T_\text{c}`$ the system is close to transition states rather than close to minima Angelani et al. (2000); Broderix et al. (2000); Grigera et al. (2002). In our earlier publication Doliwa and Heuer (2003a), we have already pointed out the technical shortcomings of the latter kind of investigations. More importantly, though, the above-cited works do not incorporate the MB structure of configuration space. In doing so, one has focused on the vast majority of intra-MB transitions, which, however, are irrelevant for relaxation.

Finally we note that the results of our work with respect to the crossover to activated behavior as well as the non-topographic nature of inter-MB dynamics is compatible with the implications obtained from a recent analysis of spin-facilitated models of the glass transition Berthier and Garrahan (2003).

We thank C.A. Angell, L. Berthier, J.P. Garrahan, and H.W. Spiess for helpful discussions.


format: text

text:
# Topological thermal instability and the length of proteins

## I INTRODUCTION

The role of geometry has recently been considered as a factor of primary importance for the comprehension of several physical properties of proteins and other biological macromolecules. In particular, since the topology of folded states is known to influence the folding properties Go ; Plaxco ; Makarov ; Baker ; Riddle ; Clementi ; Thirumalai ; Cecconi  , a great deal of work has been devoted to the study of topological aspects describing the networks of links between aminoacids in folded proteins. Kaba ; Park ; Vendru . Futhermore, relevant features of protein conformations seem to follow the geometrical principles of the optimal packing problem Maritan ; Banavar .

Starting from the primary, linear, structure (the sequence of aminoacids), a protein evolves during the folding process until it reaches a final state (native state) whose geometrical shape is crucial to the function of the protein itself. However, the problem of the topological arrangement of proteins in their native states cannot be regarded as purely static issue. Indeed, a massive accumulation of experimental data collected from X-ray, NMR and neutron spectroscopy, has reveiled that protein native states are rather dynamic structures where aminoacids constantly move around equilibrium positions. This dynamics, crucially involved in protein functions Frauen ; Frauen2 , is usually examined and investigated through normal modes analysis (NMA) NMA  and its variants Essential .

The study of collective motions of large proteins is generally of difficult accessibility to realistic all-atoms NMA Levitt  and simplified or approximate approaches are usually welcome. Tirion Tirion  first proposed the possibility to replace, in protein normal mode computations, complicated empirical potentials by Hookian pairwise interactions depending on a single parameter. The approach stems from the observation that low-frequency dynamics, which is mainly associated with protein-domain motion, is generally insensitive to the finer details of atomic interactions. A vast further literature Hinsen ; Bahar ; Atilgan ; Haliloglu  confirmed the success of simple harmonic models in the study of slow vibrational dynamics of large biological macromolecules, until they become a viable alternative to heavy and time-consuming all-atoms NMA. This success relies on the striking agreement with experiments, on the presence of few adjustable parameters and finally on their easy numerical implementation on computers and fast result production. For such reasons harmonic models are suitable also for the systematic analysis of large datasets of proteins.

The topological stability of macromolecules is far from being a pure mechanical problem as it closely involves thermodynamics: indeed, the relevant thermodynamic potential to minimize in order to find the stable configuration is not an energy, but a free energy. This is due to the interaction with the environment (schematized as a thermal bath) which is generally not negligible, specifically for biological macromolecules having a stable phase in a solvent. In particular, water is a very efficient medium for the transfer of thermal energy at microscopic scales (i.e. oscillations and molecular rotations).

Following these considerations, in this work we apply the approach of NMA to investigate the influence of the global native state topology on the thermal stability of proteins.

Vibrational thermal instability is a well-known issue in solid state physics. Since the first classical analysis of Peierls Peierls , it has been recognized that the equilibrium with a thermal bath can dramatically influence the allowed topological arrangement of large geometrical structures. The most striking consequence of Peierls instability has concerned up to now low-dimensional crystals: for one and two-dimensional lattices the mean square displacement of a single atom at finite temperature diverges in the thermodynamic limit, i.e. with increasing number of atoms. When the displacement exceeds the order of magnitude of the lattice spacing, the topological arrangement of the lattice is unstable and the crystal becomes a liquid. For real structures, formed by a finite number of units and far from the thermodynamic limit, the divergence sets a maximal stability size, which is negligible for one-dimensional lattices and typically mesoscopic for two-dimensional lattices.

However, thermal instability is present not only in crystals but also in the case of structurally inhomogeneous systems, such as, glasses, fractals, polymers and non crystalline structures. There, the problem is much more complex. Generalizing the Peierls approach to mesoscopic disordered structures, we are able to apply this kind of arguments to the thermal stability of macromolecules. In the following we describe how this can be done in the case of proteins; we will predict the existence of a critical stability size depending on a global topological parameter (the spectral dimension) and compare our predictions with experimental data.

## II THEORY

In a recent paper BCFV  generalizing the Peierls’ result, we have shown that a thermodynamic instability appears also in inhomogeneous structures and it is determined by the spectral dimension $`\overline{d}`$. The parameter $`\overline{d}`$ Orbach  is defined according to the asymptotic behaviour of the density of harmonic oscillations at low frequencies. More precisely, denoting by $`\rho (\omega )`$ the density of modes with frequency $`\omega `$, then

$$\rho (\omega )\omega ^{\overline{d}1}$$
(1)
for $`\omega 0`$. The spectral dimension is the most natural extension of the usual Euclidean dimension $`d`$ to disordered structures as far as dynamical processes are concerned. It coincides with $`d`$ in the case of lattices, but in general, it can assume non-integer values between $`1`$ and $`3`$. The spectral dimension represents a useful measure of the effective connectedness of geometrical structures at large scales, because large values of $`\overline{d}`$ correspond to high topological connectedness. Moreover, it characterizes not only harmonic oscillations, but it is also related to diffusion, phase transitions and electrical conductivity, allowing a variety of both experimental and numerical methods for its determination Burioni ; Saviot . The relevance of $`\overline{d}`$ in connection with the anomalous density of vibrational modes in proteins has also been considered in refs.Avraham ; Karplus .

In the case of thermal instability, we demostrated that, for $`\overline{d}2`$, the mean square displacement $`r^2`$ of a structural unit (being an atom, a molecule or a supramolecular structure according to the studied case) of a system composed of $`N`$ elements, diverges in the limit $`N\mathrm{}`$. Denoting with $`T`$ the temperature of the heat bath and with $`K_B`$ the Boltzmann constant, the divergence is given by the asymptotic law:

$$r^2K_BTN^{2/\overline{d}1}$$
(2)
for $`\overline{d}<2`$. For $`\overline{d}=2`$, the mean square displacement diverges logarithmically, $`r^2K_BT\mathrm{ln}N`$, as in the case of the Peierls’ result for a two dimensional crystal. Notice that the divergence in $`r^2`$ is only determined by $`\overline{d}`$. Now, at any given temperature $`T`$, there will exist a threshold value $`N(T)`$ beyond which $`r^2^{1/2}`$ exceeds the typical spacing between nearest neighbors, making the solid structure unstable. Therefore at large enough $`N`$, the solid will experience a structural reorganization which can lead either to a homogeneous liquid phase at sufficiently high temperatures or to a disordered 3-dimensional solid, which is homogeneous at large scales and inhomogeneous at small scales. In general, the threshold values of $`N`$ are very small with respect to the typical order of magnitude of macroscopic systems, being rather comparable to the size of large complex macromolecules such as biopolymers.

This poses an intriguing question concerning proteins. Indeed, to exploit their biological function proteins must keep a specified geometrical and topological arrangement and cannot afford any, even partial, geometrical large scale fluctuations as it happens, instead, to swallen polymeric chains in a good solvent Degennes . This makes thermodynamical stability problems of crucial importance and suggests a possible correlation between the spectral dimension and the length of protein chains.

The vibrational stability problem in proteins has been analysed through the gaussian network model (GNM), proposed in Tirion , which yields results in agreement with principal X-ray spectroscopy experiments. This approach generally considers proteins as elastic networks, whose nodes are the positions of the alpha-carbons (Cα) in the native structure and the interactions between nodes are assimilated to harmonic springs. The only information required to implement the method is the knowledge of the native structure, and only two free parameters are introduced, the spring constant and interaction cutoff. The GNM is defined by the quadratic Hamiltonian

$$H=\underset{i}{\overset{N}{}}\frac{𝐩_i^2}{2m}+\frac{\gamma }{2}\underset{ij}{}\mathrm{\Delta }_{ij}(\delta 𝐫_i\delta 𝐫_j)^2$$
(3)
where the first term is the kinetic energy of the system, $`\gamma `$ being the strength of the springs that are assumed homogeneous, $`𝐑_i`$ and $`\delta 𝐫_i`$ indicating the equilibrium position and the displacement with respect to $`𝐑_i`$ of the $`i`$-th Cα atoms. The model is eventually defined by the contact matrix $`\mathrm{\Delta }`$ with entries: $`\mathrm{\Delta }_{ij}=1`$ if the distance $`|𝐑_i𝐑_j|`$ between two Cα’s, in the native conformation, is below the cutoff $`R_0=0.6`$ Å, while is $`0`$ otherwise.

The harmonic spectrum for each structure is given by the set of eigenvalues $`\{\omega _1,\mathrm{},\omega _N\}`$ of the Kirchhoff matrix (or valency-adjaciency matrix) $`\mathrm{\Gamma }_{ij}=\mathrm{\Delta }_{ij}+\delta _{ij}_{li}\mathrm{\Delta }_{il}`$,

Notice that the first eigenvalue $`\omega _1`$ vanishes and corresponds to the constant eigenvector, describing the trivial uniform translation.

The comparison between experimental data and GNM results is obtained via the X-ray crystallographyc B-factors, measuring the mean square fluctuation of Cα atoms around native positions

$$B_i(T)=\frac{8\pi ^2}{3}\delta 𝐫_i\delta 𝐫_i$$

In the GNM approximation, these can be expressed in terms of the diagonal part of the inverse of $`\mathrm{\Gamma }`$ Bfact :

$$\delta 𝐫_i\delta 𝐫_j=\frac{3K_BT}{\gamma }[\mathrm{\Gamma }^1]_{ij}$$

The comparison with crystallographyc data is crucial for setting the correct values of the free parameters $`R_0`$ end $`\gamma `$ of the GNM (see. Methods and Results).

## III METHODS AND RESULTS

We present now the harmonic analysis of the GNM performed over the dataset of protein native structures with different sizes listed in Table 1. and downloaded from the Brookheaven Protein Data Bank. We selected our representative statistical sample according to the following criteria: first, we only considered proteins with a stable large scale geometry. This excludes multiple domains proteins, where the domains can move independently giving rise to large geometrical fluctuations. Moreover, we selected proteins not binded to fragments of DNA, RNA or other substrates because such structures cannot be described with sufficient accuracy in terms of simple harmonic model with only two effective parameters.

Finally, we choose proteins covering uniformly a wide length interval ranging from $`100`$ to $`3600`$ to test our prediction.

The diagonalization of the Kirchhoff matrix $`\mathrm{\Gamma }`$ to obtain its eigenvalues $`\{\omega _1^2,\mathrm{},\omega _n^2\}`$ and eigenvectors has been performed with the standard numerical packages Recipes . The knowledge of the eigenvectors (eigenmodes) allows to compute the GNM B-factors:

$$B_i(T)=\frac{8\pi ^2K_BT}{\gamma }\underset{k}{}\frac{|u_i(k)|^2}{\omega _k^2}$$

where $`i`$ is the residue index, the sum runs over all non-zero frequencies and $`𝐮(𝐤)_i`$ indicates the $`i`$-th component of the $`k`$-th eigenmode.

When generating the contact matrix $`\mathrm{\Delta }`$, the value of the interaction cutoff $`R_0`$ has been chosen in order to have agreement between the experimental and the GNM B-factors for each proteins, as showed in Fig. 1. We found that the average value $`R_0=6\AA `$ was a consistent choice for the whole set of proteins.

The purpose of our analysis is basically:

a) investigating whether there exists a correlation between the spectral dimension of native structures and the length of natural occurring proteins;

b) in the case of positive result for point a), verifying whether the correlation is connected with equation (2).

This expression establishes a rather strong relation between the spectral dimension and the maximum size $`N_{max}`$ of a protein. Since the stability is supposed to fail when the fluctuation $`r^2^{1/2}`$ becomes of the same order of magnitude of the mean distance between non consecutive aminoacids ($`6`$ Å), one has

$$\frac{1}{\mathrm{ln}(N_{max})}=a(2/\overline{d}1)$$
(4)
The proportionality constant $`a`$ in (4) depends on the mean aminoacid spacing, on the spring elastic constant g and temperature T. However, this dependece is very weak (i.e.only logarithmic). This allows for a comparison of different proteins without the computation of the specific parameters.

Actually, equation (4), being based uniquely on thermodynamics stability, can be regarded as an upper bound prediction only.

In Figure 2, we plot the harmonic spectrum, obtained within the GNM, for three proteins with sizes, small, medium and large, respectively. The low frequency regions clearly exhibit a power law behavior whose exponent $`\beta `$ is directly related to the spectral dimension via formula $`\overline{d}=2/\beta `$.

Figure 3 summarises the results of our analysis. We plot the $`1/\mathrm{ln}(N)`$ versus the computed slope from the spectra. The choice of these variables is suggested by relation (4): indeed, if eqn.(3) holds, we should obtain a straight line crossing the origin.

As matter of fact, our data are very well fitted by a straight line, with linear regression coefficient $`0.76`$, but described by the equation

$$\frac{1}{\mathrm{ln}(N)}=a(2/\overline{d}1)+b$$
(5)
with parameters $`a=0.133`$ and $`b=1.134`$.

## IV Discussion and Conclusions

The neat result expressed by (5) deserves several comments. First, this equation is in agreement with the upper bound we obtained Eq. 4, supporting the relevance of topological thermal instability in constraining the protein geometry. More important, not only the upper bound is satisfied, but the experimental points lie on a straight line parallel to the upper bound one Eq. (4). This suggests a more fundamental role of topological stability: in some sense, the protein tends to arrange topologically in such a way to reach the minimum value compatible with stability constraints. In other words, for any fixed length, it tends to the most swollen state which remains stable with respect to thermal fluctuations.

A very interesting point is the meaning of the parameter $`b`$ which is 0 in Eq. (4). Its positive value could have different motivations, but its universal nature (it is a “protein-independent” global shift) must be due to a very general mechanism. A rather obvious reason is the contribution of anharmonic interactions at finite temperatures; a more intriguing one could be an effective longer range interaction due to the presence of bound water molecules around the external aminoacids, which could change the effective form of the interaction matrix $`\mathrm{\Lambda }`$. This hypothesis is also suggested by the physical interpretation of $`b`$ as an anomalous dimension exponent, typically related to a renormalized interaction Golden

However, the most intriguing evidence relies on the regression coefficient Independently of the physical origin of $`b`$, its high value strongly supports the existence of a thermodynamic stability threshold, dependent on the topology of the folded state, for the size of proteins.


text:
# Evaporation residue cross-sections as a probe for nuclear dissipation in the fission channel of a hot rotating nucleus

## Abstract

Evaporation residue cross-sections are calculated in a dynamical description of nuclear fission in the framework of the Langevin equation coupled with statistical evaporation of light particles. A theoretical model of one-body nuclear friction which was developed earlier, namely the chaos-weighted wall formula, is used in this calculation for the 224Th nucleus. The evaporation residue cross-section is found to be very sensitive to the choice of nuclear friction. The present results indicate that the chaotic nature of the single-particle dynamics within the nuclear volume may provide an explanation for the strong shape-dependence of nuclear friction which is usually required to fit experimental data.

1 Introduction

Fission of highly excited compound nuclei produced in heavy ion induced fusion reactions has evoked considerable interest in the recent years mainly due to the fact that it cannot be accounted for by the statistical model of Bohr and Wheeler for nuclear fission . In particular, it is now established that the multiplicities of neutrons, light charged particles and photons emitted by a hot compound nucleus are much higher than those predicted by the statistical model . This implies that the fission life time is considerably underestimated when using the Bohr-Wheeler description of nuclear fission. Consequently, dissipative dynamical models for fission of excited nuclei were developed following the original work of Kramers who considered fission dynamics to be analogous to that of diffusion of Brownian particles over a barrier . These dynamical models were found to be successful in reproducing a large body of experimentally measured multiplicity data.

Dynamical effects are also found essential to calculate the evaporation residue cross-section of highly excited compound nuclei . Specifically, the decay width for an overdamped system was subsequently found to be more appropriate for nuclear fission and this fission width is generally used in the description of the statistical decay of an excited compound nucleus . It turns out that the evaporation residue cross-section depends strongly on the strength of the nuclear dissipation whenever it is a very small fraction of the total fusion cross-section. In fact, the evaporation residue cross-section in such cases can serve the purpose of a sensitive probe for nuclear dissipation.

The dissipative force in the dynamics of fission arises due to the interaction between the large number of intrinsic nuclear degrees of freedom and the few collective or fission degrees of freedom. The strength of the dissipative force is usually treated as an adjustable parameter in order to fit the experimental data. Fröbrich et al.  obtained a phenomenological shape-dependent nuclear dissipation which was able to reproduce the fission probability and prescission neutron multiplicity excitation functions for a number of compound nuclei. The magnitude of this phenomenological dissipation was found to be very small for compact shapes of the nucleus, but a strong increase in its value was needed for large deformations. In a recent work, Diószegi et al.  have analyzed the $`\gamma `$ as well as neutron multiplicities and evaporation residue cross-section of 224Th and have concluded that the experimental data can be fitted equally well with either a temperature or a deformation-dependent nuclear dissipation. Interestingly, the deformation-dependence of the above dissipation also corresponds to a lower value of the strength of the dissipation inside the saddle and a higher value outside the saddle, similar to the phenomenological dissipation of ref.  mentioned in the above.

A theoretical model for nuclear dissipation, namely the wall formula, was developed long ago by Blocki et al.  in a simple classical picture of one-body dissipation. However, the strength of the wall dissipation was found to be considerably higher than that required to fit experimental data. It is only recently that a modification has been incorporated into the wall formula which resulted in a shape-dependent reduction factor in the strength of the friction. The modification essentially arose out of a closer examination of one crucial assumption of the wall formula concerning the randomization of particle motion within the nuclear volume. It was assumed in the original wall formula that the particle motion is fully randomized. This full randomization assumption is relaxed in the modified wall formula in order to make it applicable to systems with partly randomized or chaotic single-particle motion. In what follows, we shall use the term “chaos-weighted wall formula” (CWWF) for this modified dissipation in order to distinguish it from the original wall formula (WF) dissipation. As was shown in ref., the CWWF dissipation coefficient $`\eta _{cwwf}`$ will be given as

$`\eta _{cwwf}=\mu \eta _{wf},`$ (1)
where $`\eta _{wf}`$ is the dissipation or friction coefficient as was given by the original wall formula  and $`\mu `$ is a measure of chaos (chaoticity) in the single-particle motion and depends on the instantaneous shape of the nucleus. The value of chaoticity $`\mu `$ changes from 0 to 1 as the nucleus evolves from a spherical shape to a highly deformed one. The CWWF dissipation is thus much smaller than the WF dissipation for compact nuclear shapes while they become closer at large deformations. It is worthwhile to note here that the above shape-dependent CWWF has features similar to the empirical dissipations discussed in the last paragraph. In a recent work, we have shown that the prescission neutron multiplicity and fission probability calculated from Langevin dynamics using the chaos-weighted wall dissipation agree fairly well with the experimental data for a number of heavy compound nuclei ($`A200`$) over a wide range of excitation energies . This strongly suggests that the shape-dependent CWWF can be considered as a suitable theoretical model for one-body dissipation in nuclear fission. Further, a lack of full randomization or chaos in the single-particle motion can provide a physical explanation for the empirical requirement for reduction in strength of friction for compact nuclear shapes in order to fit experimental data.

In the present work, we shall employ the chaos-weighted wall formula dissipation to calculate the evaporation residue excitation function for the 224Th nucleus. Our main motivation here will be to put the chaos-weighted wall formula to a further test and verify to what extent it can account for the experimental evaporation residue data which is a very sensitive probe for nuclear dissipation. In our calculation, we shall first assume that a compound nucleus is formed when a projectile nucleus completely fuses with a target nucleus in a heavy ion collision. Processes such as fast fission or quasifission thus cannot be described in the model considered here. We shall describe the fission dynamics of the compound nucleus by the Langevin equation while the light particles and photons will undergo statistical emission. The Langevin equation will be solved by coupling it with particle and $`\gamma `$ evaporation at each step of its time evolution. The prescission particle multiplicity and fission probability will be obtained by sampling over a large number of Langevin trajectories. The chaos- weighted wall friction coefficient is obtained following a specific procedure  which explicitly considers particle dynamics in phase space in order to calculate the chaoticity factor $`\mu `$ of eq.(1). There is no free parameter in this calculation of friction. The other input parameters for the dynamical calculation are obtained from standard nuclear models. Calculation will be performed at a number of excitation energies for 224Th formed in the 16O+208Pb system. We have chosen this system essentially because of the availability of experimental data on both evaporation residue and prescission neutron multiplicity covering the same range of excitation energies and the fact that earlier analyses of the evaporation residue excitation function have already indicated the need for a dynamical model for fission of this nucleus .

In the following section, we shall briefly describe the dynamical model along with the necessary input as used in the present work. The details of the calculation will also be given here. The calculated evaporation residue excitation function and prescission neutron multiplicities will be compared with the experimental values in sect.3. A summary of the results along with the conclusions can be found in the last section.

2 Langevin description of fission

An appropriate set of collective coordinates to describe the fission degree of freedom consists of the shape parameters $`c,h`$ and $`\alpha `$ as was suggested by Brack et al.  and we shall employ them in the present calculation. We further simplify our calculation by considering only symmetric fission ($`\alpha =0`$) since the compound nucleus 224Th is much heavier than the Businero-Gallone transition point. The potential landscape in $`(c,h)`$ coordinates is generated from the finite-range liquid drop model  where we calculate the generalized nuclear energy by double folding the uniform density with a Yukawa-plus-exponential potential. The Coulomb energy is obtained by double folding another Yukawa function with the density distribution. We shall further assume in the present work that fission would proceed along the valley of the potential energy landscape. Consequently we shall use an effective one-dimensional potential in the Langevin equation which will be defined as $`V(c)=V(c,h)`$ at valley. This will reduce the problem to one dimension in order to simplify the computation. The Langevin equations in one dimension will thus be given  as

$`{\displaystyle \frac{dp}{dt}}`$ $`=`$ $`{\displaystyle \frac{p^2}{2}}{\displaystyle \frac{}{c}}\left({\displaystyle \frac{1}{m}}\right){\displaystyle \frac{F}{c}}\eta \dot{c}+R(t),`$ (2)
$`{\displaystyle \frac{dc}{dt}}`$ $`=`$ $`{\displaystyle \frac{p}{m}}.`$ (3)
The shape-dependent collective inertia and the dissipation coefficients in the above equations are denoted by $`m`$ and $`\eta `$ respectively. $`F`$ is the free energy of the system while $`R(t)`$ represents the random part of the interaction between the fission degree of freedom and the rest of the nuclear degrees of freedom considered collectively as a thermal bath in the present picture. The collective inertia, $`m`$, will be obtained by making the Werner-Wheeler approximation  assuming an incompressible irrotational flow. The driving force in a thermodynamic system should be derived from its free energy which we will calculate considering the nucleus as a noninteracting Fermi gas . The instantaneous random force $`R(t)`$ is assumed to have a stochastic nature with a Gaussian distribution whose average is zero . The strength of the random force will be determined by the dissipation coefficient through the fluctuation-dissipation theorem. The various input quantities are described in some detail in a recent publication .

2.1 Nuclear dissipation

We shall use the chaos-weighted one-body wall-and-window dissipation in the present calculation. We shall beiefly describe here the essential features of this dissipation, the details of which may be found elsewhere . In the wall-and-window model of one-body dissipation, the window friction is expected to be effective after a neck is formed in the nuclear system . Further, the radius of the neck connecting the two future fragments should be sufficiently narrow in order to make the energy transfer irreversible. It therefore appears that the window friction should be very nominal when neck formation just begins. Its strength should increase as the neck becomes narrower reaching its classical value when the neck radius becomes much smaller than the typical radii of the fragments. We shall approximately describe the above scenario by defining a transition point $`c_{win}`$ in the elongation coordinate at which the window friction will be turned on. We shall also assume that the compound nucleus evolves into a binary system beyond $`c_{win}`$ and accordingly correction terms for the motions of the centers of mass of the two halves will be applied to the wall formula for $`c>c_{win}`$ . However, it may be noted that the window dissipation and the center of mass motion correction tend to cancel each other to some extent. Consequently, the resulting wall-and-window friction is not very sensitive to the choice of the transition point. We shall choose a value for $`c_{win}`$ at which the nucleus has a binary shape and the neck radius is half of the radius of either of the would-be fragments.

The wall-and-window dissipation and its chaos-weighted version will thus be given as

$`\eta _{wf}(c<c_{win})=\eta _{wall}(c<c_{win}),`$ (4)
and

$`\eta _{wf}(cc_{win})=\eta _{wall}(cc_{win})+\eta _{win}(cc_{win}),`$ (5)
while

$`\eta _{cwwf}(c<c_{win})=\mu (c)\eta _{wall}(c<c_{win}),`$ (6)
and

$`\eta _{cwwf}(cc_{win})=\mu (c)\eta _{wall}(cc_{win})+\eta _{win}(cc_{win}).`$ (7)
The detailed expressions for the wall-and-window frictions can be found in ref. .

The chaoticity $`\mu (c)`$, introduced in the earlier section, depends on the instantaneous shape of the nucleus . In a classical picture, this will be given as the average fraction of the nucleon trajectories within the nucleus which are chaotic when the sampling is done uniformly over the nuclear surface. The value of the chaoticity for a given nuclear shape is evaluated by sampling over a large number of classical nucleon trajectories while each trajectory is identified either as a regular or as a chaotic one by considering the magnitude of its Lyapunov exponent and the nature of its variation with time . The shape-dependence of the chaoticity, thus obtained, is shown in fig.1. Further, defining a quantity $`\beta (c)=\eta (c)/m(c)`$ as the reduced dissipation coefficient, its dependence on the elongation coordinate is also shown in fig.1 for the 224Th nucleus. It can be immediately noticed that the CWWF is strongly suppressed compared to the WF dissipation for near-spherical shapes ($`c1`$) and this can be qualitatively understood as follows. A particle moving in a spherical mean field represents a typical integrable system and its dynamics is completely regular. When the boundary of the mean field is set into motion (as in fission), the energy gained by the particle at one instant as a result of a collision with the moving boundary is eventually fed back to the boundary motion in the course of later collisions. An integrable system thus becomes completely nondissipative in this picture resulting in a vanishing dissipation coefficient. It may be noted that the suppression of friction strength in CWWF is qualitatively similar to the shape-dependent frictions found empirically  to fit experimental data. The considerations of chaos (or rather lack of it) in particle motion can thus provide a physical explanation for the reduction in friction strength required for compact shapes of the compound nucleus.

The strong shape-dependence of the CWWF can have some interesting consequences. In a dynamical description of fission, a compound nucleus spends most of its time in undergoing shape oscillations in the vicinity of its ground state shape before it eventually crosses the saddle and proceeds towards the scission point. Since the spin of a compound nucleus formed at a small excitation is also small, its ground state shape is nearly spherical and in this region the CWWF friction is also small. Conversely, higher spin values are mostly populated in a highly excited compound nucleus making its ground state shape highly deformed and thus it experiences a strong CWWF friction. Therefore, if one uses a shape-independent friction in a dynamical model of fission, its strength has to increase with increasing temperature in order to give an equivalent description to that provided by the temperature-independent but shape-dependent CWWF friction. In fact, it was observed in ref.  that a shape-dependent friction fits the experimental data equally well to that achieved by a strong temperature-dependent friction. Since there is a physical justification for shape-dependence in nuclear friction from chaos considerations, it is quite likely that the above strong temperature-dependence, at least a substantial part of it, is of dynamical origin as explained in the above and thus is an artifact arising out of using a shape-independent friction.

It must be pointed out, however, that one would expect a temperature-dependence of nuclear friction from general considerations such as larger phase space becoming accessable for particle-hole excitations at higher temperatures. In a microscopic model of nuclear friction using nuclear response function, Hofmann et al.  have obtained a nuclear friction which depends upon temperature as $`0.6T^2`$ (leading term). This may be compared with the empirical temperature-dependent term of $`3T^2`$ which was found in ref.. It therefore appears that only a small fraction of the empirical temperature-dependence can be accounted for by the inherent temperature-dependence of nuclear friction while the rest of it has a dynamical origin as we have discussed in the above.

In the present work, we shall not consider any empirical temperature-dependence of the CWWF or WF frictions in order to study solely the effects of shape-dependence. In what follows, we shall use both the WF and CWWF dissipations in a dynamical model of fission and shall investigate the effect of the reduction in the CWWF strength on the evaporation residue cross-section.

2.2 Dynamical model calculation

In the dynamical model calculation, the initial spin and excitation energy of the compound nucleus is determined from the entrance channel specifications in the following manner. The fusion cross-section of the target and projectile in the entrance channel usually obeys the following spin distribution

$`{\displaystyle \frac{d\sigma (l)}{dl}}={\displaystyle \frac{\pi }{k^2}}{\displaystyle \frac{(2l+1)}{1+\mathrm{exp}\frac{(ll_c)}{\delta l}}}`$ (8)
where we shall obtain the parameters $`l_c`$ and $`\delta l`$ by fitting the experimental fusion cross-sections. The initial spin of the compound nucleus in our calculation will be obtained by sampling the above spin distribution function. The total excitation energy ($`E^{}`$) of the compound nucleus can be obtained from the beam energy of the projectile, and energy conservation in the form

$$E^{}=E_{int}+V(c)+p^2/2m$$
(9)
gives the intrinsic excitation energy $`E_{int}`$ and the corresponding nuclear temperature $`T=(E_{int}/a)^{1/2}`$ where $`a`$ is the nuclear level density parameter. The centrifugal potential is included in $`V(c)`$ in the above equation.

We shall use the following level density parameter due to Ignatyuk et. al.  which incorporates the nuclear shell structure at low excitation energy and goes smoothly to the liquid drop behavior at high excitation energy,

$`a(E_{int})=\overline{a}(1+{\displaystyle \frac{f(E_{int})}{E_{int}}}\delta M),`$ (10)
with

$`f(E_{int})=1exp(E_{int}/E_D)`$ (11)
where $`\overline{a}`$ is the liquid drop level density parameter, $`E_D`$ determines the rate at which the shell effects disappear at high excitations, and $`\delta M`$ is the shell correction given by the difference between the experimental and liquid drop masses, $`(\delta M=M_{exp}M_{LDM})`$. We shall further use the shape-dependent liquid drop level density parameter given as ,

$`\overline{a}(c)=a_vA+a_sA^{\frac{2}{3}}B_s(c)`$ (12)
where we choose the values for the parameters $`a_v`$, $`a_s`$ and the dimensionless surface area $`B_s`$ following ref..

In order to solve the Langevin equations, the initial values of the coordinates and momenta $`(c,p)`$ of the fission degree of freedom are obtained from sampling random numbers following the Maxwell-Boltzmann distribution. The Langevin equations (eq.(2)) are subsequently numerically integrated following the procedure outlined in ref.. At each time step of integration of the Langevin equations, particle (neutron, proton and alpha) and giant dipole $`\gamma `$ evaporation will be considered following a Monte Carlo sampling technique . For this purpose, the particle and $`\gamma `$ decay widths are calculated using the inverse cross-section formula as given in ref.. After each particle emission, the potential energy landscape of the parent nucleus is replaced by that of the daughter nucleus. The intrinsic excitation energy, mass and spin of the compound nucleus are also recalculated after each emission. The spin of the compound nucleus is reduced only in an approximate way by assuming that each neutron, proton or a $`\gamma `$ carries away $`1\mathrm{}`$ angular momentum while that carried by an alpha particle is $`2\mathrm{}`$. However, only neutron emission is found to be relevant for the yield of the residue cross-section.

A Langevin trajectory will be considered as undergone fission if it reaches the scission point ($`c_{sci}`$) in course of its time evolution. Alternately it will be counted as an evaporation residue event if the intrinsic excitation energy becomes smaller than either the fission barrier or the binding energy of a particle. The calculation proceeds until the compound nucleus undergoes fission or becomes an evaporation residue. This calculation is repeated for a large number of Langevin trajectories and the evaporation residue formation probability is obtained as the fraction of the trajectories which have ended up as evaporation residues. The evaporation residue cross-section is subsequently obtained by multiplying the experimental value for fusion cross-section in the entrance channel with the above formation probability of the evaporation residue. Similarly, the average number of particles (neutrons, protons or alphas) emitted in the fission events will give the required prescission particle multiplicities.

Following the fission dynamics through the Langevin equation during the entire life time of a compound nucleus can however take an extremely long computer time. As an alternative procedure, we shall first follow the time evolution of a compound nucleus according to the Langevin equations for a sufficiently long period during which a steady flow across the fission barrier is established. Beyond this period, a statistical model for compound nucleus decay is expected to be a equally valid and more economical in terms of computation. We shall therefore switch over to a statistical model description after the fission process reaches the stationary regime. This combined dynamical and statistical model, first proposed by Mavlitov et al.  , however, requires the fission width along with the particle and $`\gamma `$ widths in the statistical branch of the calculation. This fission width should be the stationary limit of the fission rate as determined by the Langevin equation. However, it is not possible to obtain this fission rate analytically for the strongly shape-dependent CWWF and WF disipations. We shall therefore use a suitable parametric form of the numerically obtained stationary fission widths using the CWWF (and also WF) dissipations in order to use them in the statistical branch of our calculation. The details of this procedure is given in ref. which we shall follow to calculate all the required fission widths for the present work.

3 Results

We have calculated the prescission neutron multiplicity and the evaporation residue (ER) cross-section for the compound nucleus 224Th when it is formed in the fusion of an incident 16O nucleus with a 208Pb target nucleus. The calculation is done at a number of incident energies in the range of 80 MeV to 140 MeV using both the WF and the CWWF dissipations. Figure 2 shows the calculated prescission neutron multiplicity along with the experimental data . Both the WF and CWWF predictions for multiplicity are quite close to the experimental values and this shows that neutron multiplicity is not very sensitive to the dissipation in fission in the energy range under consideration. It must be pointed out, however, that the CWWF predictions for neutron multiplicity are closer to experimental data compared to those from WF at much higher excitations of the compound nucleus.

It may be mentioned at this point that though we calculate the number of prescission protons, alphas and GDR $`\gamma `$’s, we do not compare them with experimental data because these numbers are rather small with large statistical uncertainties in the present work. In order to obtain the energy spectrum of the $`\gamma `$ multiplicity with a reasonable statistical accuracy, in particular, it is necessary to perform computation using a much larger ensemble of trajectories than the one (comprising of 20000 trajectories) used here. This puts a severe demand on computer time making such computations impractical at present. However, an alternative approach would be to make use of the time-dependent fission widths in a full statistical model calulation of compound nucleus decay. This calculation would be much faster than the present Langevin dynamical model calculation though the time-dependent fission widths would be required as input to this statistical model calculation. We plan to perform such calculations in future.

We shall next consider the results of the evaporation residue calculation. Figuer 3 shows the evaporation residue excitation functions calculated using both the WF and CWWF dissipations. The experimental values of evaporation residue cross-section are also shown in this figure . We first note that the calculated evaporation residue cross-section is very sensitive to the dissipation in the fission degree of freedom. The WF predictions are a few times (typically 2-5) larger than those obtained with the CWWF dissipation. Next we make the important observation that the CWWF predicted excitation function is much closer to the experimental values than that obtained with the WF dissipation. This observation clearly shows that the chaos-weighted factor in CWWF changes its strength in the right direction. We must take note of the fact, however, that the CWWF still considerably overestimates the ER cross-section. Since the present dynamical calculation considers only one (elongation) fission degree of freedom, it is expected that inclusion of the neck degree of freedom will increase the fission probability  further and hence reduce the ER cross-section. We plan to extend our work in this direction in future. We further observe that while a peak appears in the experimental excitation function at about 85 MeV, the same is shifted by 10 MeV in the calculated results. We do not have any explanation for this discrepancy except pointing out that there is no free parameter in our calculation and thus no parameter tuning has been attempted in order to fit experimental data. A similar shift has also been observed in an earlier work .

The structure of the evaporation residue excitation function can also reveal certain interesting features. Since the calculated values of the evaporation residue cross-section are obtained as the product of the fusion cross-section and the probability of evaporation residue formation, the initial rise of the ER cross-section with beam energy essentially reflects the steep rise of fusion cross-section in this energy region . At still higher beam energies, the ER cross-section becomes approximately stable which results from a delicate balance between the increasing trend of the fusion cross-section and the decreasing trend of the probabilty of ER formation. Had the ER formation probability decreased at a rate higher than those obtained in the present calculation, the resulting ER cross-section would have decresed at higher compound nuclear excitations. In fact, such an observation was made in ref. where the ER cross-section obtained from standard statistical model calculation was found to decrease very steeply beyond 100 MeV of beam energy. In order to explore this point a little further, we have calculated the excitation function of the partial width for fission. Since fission can take place at any stage during neutron (or any other light particle) evaporation, the partial widths are calculated for 222Th and 223Th as well at excitation energies reduced by the neutron separation energy after each neutron emission. The compound nuclear spin was taken as $`l_c`$ from eq.(7) while only CWWF dissipation was considered for this calculation.

The calculated excitation functions of the fission partial widths are shown in fig. 4. The calculated values of the ER formation probability ($`P_{ER}`$) are also displayed in this figure. Each partial width excitation function is found to have a minimum around 90 MeV of beam energy after which it starts increasing till this trend is arrested and reversed at higher excitations. Recalling the fact that the above results on partial widths are only indicative while $`P_{ER}`$ is obtained from a full dynamical calculation, it is of interest to note that a bump in the excitation function of $`P_{ER}`$ also appears in the above ($``$ 90 MeV) energy range. Subsequently, the value of $`P_{ER}`$ drops rather sharply before it settles to a steady value at higher excitations. This feature is also complementary to that of the excitation functions of the partial widths of fission. We thus demonstrate in a schematic manner how the structure in the excitation function of the ER cross-section is related to the competition between fission and other decay channels at different stages of fission.

4 Summary and conclusions

We have applied a theoretical model of one-body nuclear friction, namely the chaos-weighted wall formula, to a dynamical description of compound nuclear decay where fission is governed by the Langevin equation coupled with the statistical evaporation of light particles. We have used both the standard wall formula and its modified form with the chaos-weighted factor in order to calculate the prescission neutron multiplicity and evaporation residue excitation functions for the 224Th nucleus. Though the number of the prescission neutrons calculated with either WF or CWWF friction are found to be very close to each other in the energy range considered, the evaporation residue cross-section is found to depend very strongly on the choice of nuclear friction. The evaporation residue cross-section calculated with the CWWF friction gives a much better agreement with the experimental data compared to the WF predictions. This result demonstrates that the consequences of chaos in particle motion give rise to a strong suppression of the strength of the wall friction for compact shapes of the compound nucleus which ,in turn, brings theoretically calculated evaporation residue cross-sections considerably closer to the experimental values. Thus the chaos considerations may provide a plausible explanation for the shape-dependence of the strength of nuclear friction which was found  to be necessary in order to fit experimental data.


text:
# Data compression using correlations and stochastic processes in the ALICE Time Projection chamber

## I Introduction

ALICE (A Large Ion Collider Experiment) is an experiment that will start in 2007 at the LHC (Large Hadron Collider) at CERN TPCTDR ; alicewww . The experiment will study collisions between heavy ions with energies around 5.5 TeV per nucleon. The collisions will take place at the center of a set of several detectors, which are designed to track and identify the produced particles.

One of the main detectors of the ALICE experiment is the Time Projection Chamber (TPC). Its task is track finding, momentum measurement and particle identification by d$`E`$/d$`x`$. Good two-track resolution, required for correlation studies, is one of the main design goals.

The TPC is a large horizontal cylinder, filled with gas, where a suitable axial electric field is present. When particles pass through, they ionize the gas atoms, and the resulting electrons drift in the electric field. By measuring the arrival of electrons at the end of the chamber, the TPC can reconstruct the path of the original charged particles. The electrons are collected by more than 570 000 sensitive pads where they create signals. These signals are amplified by a preamplifier–shaper and digitalized by a 10-bit A/D converter at a sampling frequency of 5.66 MHz. The digitalized signal is processed and formatted by an Application Specific Integrated Circuit (ASIC) called ALTRO (ALICE TPC Read-Out) musa00altro . At this stage, the overall throughput of the 570 000 channels is around 8.4 GByte/s.

The total amount of the TPC data is expected to be about 1 PBy per year. In order to keep the complexity and cost of the data storage equipment as low as possible, we have to reduce the volume of data using suitable data compression methods. The cost reduction of the data storage system is roughly proportional to the data compression factor. Furthermore, it is better to implement the compression system in the front-end electronics at the output of the ALTRO circuit, so that the cost for the optical links, which carry data out of the chamber to the following stages of the acquisition chain, could be also reduced.More sophisticated methods for TPC data compression based on online tracking, which will be used further in data acquisition chain are developed in Bergen and Heidelberg bergen .

The use of a lossy source model, justified by the fact that generally it can provide significantly higher compression ratios compared to lossless models, has the drawback that some deterioration in the reconstruction of data must be accepted. Lossy source models have become very popular in the last decade in the field of audio and video compression for their remarkable performance. Lossy models have been carefully designed so that reconstruction distortions are not perceived using psychovisual or psychoacoustic models or they remain comparable with the intrinsic signal noise.

Obviously, for physical data, psychovisual or psychoacoustic tests are meaningless or even not applicable since the TPC signal is not to be observed by the human eye or ear. In nicolaucig02 , the compression noise introduced on the sample values by the described lossy or quasi lossless techniques has been evaluated in terms of RMS of the introduced Error (RMSE).

However, in this case, the RMSE, despite being a simple and well known distortion measure, is not very useful. The fundamental information that has to be extracted from TPC data are not sample values themselves but the physical quantities that enable the reconstruction of particle trajectories. Therefore, the correct way to evaluate the importance of the distortion introduced by the compression–decompression process has to be related to the high level information that is carried by the data. In particular, TPC data are collected with the objective of measuring particle energy and trajectory.

Therefore, the most effective way to estimate the consequences of the compression distortion error, is to observe how the extraction of energy and trajectories are affected by the compression–decompression process. A simple way to obtain these estimates is to apply the cluster finding and tracking algorithms on both simulated data and their compressed–decompressed version and compare the results.

This article is arranged in the following way. In section II, all stochastic processes relevant for particle detection in ALICE TPC are briefly described. In section III, TPC data format is specified. In section IV, different lossless compression techniques are described and their efficiencies are compared. In section V, the fast one dimensional lossy compression technique is shown and the impact of compression–decompression to the distortion of most important physical quantities is demonstrated.

## II Stochastic processes in TPC

### II.1 Ionization in gas

A charged particle that traverses the gas of the chamber leaves a track of ionization along its trajectory. The collisions with the gas atoms are purely random. They are characterized by a mean free path $`\lambda `$ between ionizing encounters, which is given by the ionization cross-section per electron $`\sigma _{\mathrm{ion}}`$ and the density N of electrons:

$$\lambda =1/(N\sigma _{\mathrm{ion}}).$$

Therefore, the number of encounters along the length L has the mean of $`L/\lambda `$, and the frequency distribution is given by Poisson distribution

$$P(L/\lambda ,k)=\frac{(L/\lambda )^k}{k!}\mathrm{exp}(L/\lambda ).$$

The mean free path $`\lambda `$ is given by the properties of the gas and by charged particle characteristics:

$$\lambda =\frac{1}{N_{\mathrm{prim}}\times f(\beta \gamma )},$$

where $`N_{\mathrm{prim}}`$ is the number of primary electrons per cm produced by a Minimum Ionizing Particle (MIP), and $`f(\beta \gamma )`$ is Bethe–Bloch curve.

### II.2 Generation of secondary electrons

The energy loss $`E_{\mathrm{tot}}`$ released in primary ionization to atomic electrons is a random variable. It can be described by Photo-Absorbtion Ionization model (PAI). In most cases, if one neglects the atomic shell structure, at sufficiently high $`E_{\mathrm{tot}}`$ (the energy where the atomic shell structure is not more important) it obeys $`1/E_{\mathrm{tot}}^2`$ rule.

If the electron produced by the charged particle has sufficient kinetic energy $`E_{\mathrm{tot}}`$, it will produce secondary electrons creating thus electron cluster. The mean total number of electrons in such cluster is given by:

$$N_{\mathrm{tot}}=\frac{E_{\mathrm{tot}}I_{\mathrm{pot}}}{W_{\mathrm{ion}}}+1,$$

where $`E_{\mathrm{tot}}`$ is the energy loss in a primary collision, $`W_{\mathrm{ion}}`$ is the effective energy required to produce an electron–ion pair and $`I_{\mathrm{pot}}`$ is the first ionization potential. The random character of the secondary ionization process smears out structures in $`E_{\mathrm{tot}}`$ spectra, atomic shell structure behavior is suppressed. For example in the gas mixture 90% Ne, 10 % CO2 the $`E_{\mathrm{tot}}^{2.2}`$ effective parametrization at lower $`E_{\mathrm{tot}}`$ can be used.

### II.3 Diffusion of electrons

Produced electrons drift through the gas with an effective constant drift velocity in the direction given by the electric field $`𝐄`$ and magnetic field $`𝐁`$ (which we assume are parallel to z-direction). Drifting electrons are scattered on the gas molecules so that their direction of motion is randomized in each collision. The position of the electron, after drifting over a distance $`L_{\mathrm{drift}}`$, can be described by 3-D Gaussian distribution:

$`P(x,y,z)={\displaystyle \frac{1}{\sqrt{2\pi }\sigma _\mathrm{T}}}\mathrm{exp}\left[{\displaystyle \frac{(xx_0)^2}{2\sigma _\mathrm{T}^2}}\right]`$
$`{\displaystyle \frac{1}{\sqrt{2\pi }\sigma _\mathrm{T}}}\mathrm{exp}\left[{\displaystyle \frac{(yy_0)^2}{2\sigma _\mathrm{T}^2}}\right]`$
$`{\displaystyle \frac{1}{\sqrt{2\pi }\sigma _\mathrm{L}}}\mathrm{exp}\left[{\displaystyle \frac{(zL_{\mathrm{drift}})^2}{2\sigma _\mathrm{L}^2}}\right],`$ (1)
where $`\{x_0,y_0,z_0\}`$ is the electron creation point and transversal diffusion $`\sigma _\mathrm{T}`$ respectively longitudinal diffusion $`\sigma _\mathrm{L}`$ are given by drift length $`L_{\mathrm{drift}}`$ and gas coefficient $`D_\mathrm{T}`$ and $`D_\mathrm{L}`$

$$\sigma _\mathrm{T}=D_\mathrm{T}\sqrt{L_{\mathrm{drift}}},$$

$$\sigma _\mathrm{L}=D_\mathrm{L}\sqrt{L_{\mathrm{drift}}}.$$

### II.4 $`𝐄\times 𝐁`$ and unisochronity effect near the anode wires

It has been assumed that the electric and magnetic fields in the drift volume are uniform and parallel. This, however, is not true close to the anode wires, where the electric field becomes radial. Thus the electrons experience a shift along the wire direction (due to the Lorentz force). If an electron enters the readout chamber at the point $`(x_\mathrm{e},y_\mathrm{e})`$, it is displaced in the x-direction (assuming that the wires are placed along y-axis). The new y-position of the electron is then given by

$$y=y_\mathrm{e}+\omega \tau (xx_\mathrm{e}),$$

where $`x`$ is the coordinate of the wire on which an electron is collected, and $`\omega \tau `$ is the tangent of Lorentz angle ($`𝐄\times 𝐁`$ effect). The drift length which determines z coordinate will be also affected, because of change in the path to the anode wire (unisochronity effect).

### II.5 Signal generation

Inside the readout chamber, as an electron drifts towards the anode wire, it travels in an increasing electric field. Once the electric field is strong enough that between collisions with the gas molecules the electron can pick up sufficient energy for ionization, another electron is created and the avalanche starts. As the number of electrons multiplies in successive generations, the avalanche continues to grow until all the electrons are collected on the wire. The resulting number of electrons created in the avalanche, can be described by an exponential probability distribution

$$P(q)=\frac{1}{\overline{q}}\mathrm{exp}\frac{q}{\overline{q}},$$

where $`\overline{q}`$ is the average avalanche amplitude.

An electron avalanche collected on the anode wire induces a charge on the pad plane. This charge is integrated over the pad area. The time signal is obtained by folding the pad response to the avalanche with the shaping function of the preampamplifier–shaper. This signal is then sampled with a constant frequency. On the top of sampled signal a random electronic noise is superimposed.

As a result a charged particle interacting with gas generates a cluster of amplitudes. This cluster is used for later estimation of local track position and of local energy deposition. The shape of the cluster is used as additional information for the estimation of position uncertainties and for the estimation of the overlap factor between two tracks.

### II.6 Accuracy of local coordinate measurement

The accuracy of the coordinate measurement is limited by a track angle which spreads ionization and by diffusion which amplifies this spread.

The track direction with respect to pad plane is given by two angles $`\alpha `$ and $`\beta `$ (see fig. 1). For the measurement along the pad-row, the angle $`\alpha `$ between the track projected onto the pad plane and pad-row is relevant. For the measurement of the the drift coordinate (z–direction) it is the angle $`\beta `$ between the track and z axis.

The ionization electrons are randomly distributed along the particle trajectory. Fixing the reference x position of a electron at the middle of pad-row, the y (resp. z) position of the electron is random variable characterized by uniform distribution with the width $`L_\mathrm{a}`$, where $`L_\mathrm{a}`$ is given by the pad length $`L_{\mathrm{pad}}`$ and the angle $`\alpha `$ (resp. $`\beta `$):

$$L_\mathrm{a}=L_{\mathrm{pad}}\mathrm{tan}\alpha $$

The diffusion smears out the position of the electron with gaussian probability distribution with $`\sigma _\mathrm{D}`$. Contribution of the $`𝐄\times 𝐁`$ and unisochronity effect is in the case of Alice TPC negligible.

The accuracy of the position measurement can be expressed as:

$`\sigma _\mathrm{z}`$ of cluster center in z (time) direction:

$`\sigma _{\mathrm{z}_{\mathrm{COG}}}^2={\displaystyle \frac{D_\mathrm{L}^2L_{\mathrm{drift}}}{N_{\mathrm{ch}}}}G_\mathrm{g}+`$ (2)
$`{\displaystyle \frac{\mathrm{tan}^2\alpha L_{\mathrm{pad}}^2G_{\mathrm{Lfactor}}(N_{\mathrm{prim}})}{12N_{\mathrm{chprim}}}}+\sigma _{\mathrm{noise}}^2`$
and $`\sigma _\mathrm{y}`$ of cluster center in y(pad) direction:

$`\sigma _{\mathrm{y}_{\mathrm{COG}}}^2={\displaystyle \frac{D_\mathrm{T}^2L_{\mathrm{drift}}}{N_{\mathrm{ch}}}}G_\mathrm{g}+`$ (3)
$`{\displaystyle \frac{\mathrm{tan}^2\beta L_{\mathrm{pad}}^2G_{\mathrm{Lfactor}}(N_{\mathrm{prim}})}{12N_{\mathrm{chprim}}}}+\sigma _{\mathrm{noise}}^2`$
where $`N_{\mathrm{ch}}`$ is the total number of electrons in cluster, $`N_{\mathrm{chprim}}`$ is the number of primary electrons in cluster, $`G_\mathrm{g}`$ is the gas gain fluctuation factor parametrization, $`G_{\mathrm{Lfactor}}`$ is the secondary ionization fluctuation factor and $`\sigma _{\mathrm{noise}}`$ describe the contribution of the electronic noise and ADC quantization to the resulting sigma of the COG.

The typical resolution in the case of ALICE TPC is on the level of $`\sigma _\mathrm{y}`$ 0.8mm and $`\sigma _\mathrm{z}`$ 1.0mm integrating over all clusters in the TPC.

### II.7 Accuracy of the total amplitude measurement

The total charge deposited in the clusters can be used for particle identification. The important value, which is specific for different particle types and different particle momenta, is the number of primary collision per unit length, $`N_{\mathrm{chprim}}`$. $`N_{\mathrm{chprim}}`$ is a random variable described by Poisson distribution. Due to the secondary ionization and gas gain fluctuations the total charge is described by very broad Landau distribution.

## III The ALICE TPC read-out data format

Before describing the compression algorithm, it is necessary to spend a few words on the format of data at the output of ALTRO circuit, in order to understand how the compression algorithms are applied. Such data are indeed the input of the compression system TPCTDR ; musa00altro .

In the ALTRO data format only the samples over a given threshold are considered, while the others are discarded. This means that, if we call bunch a group of adjacent over-threshold samples coming from one pad, the signal can be represented “bunch by bunch”. More precisely, a bunch is described by three fields: temporal information (temporal position of the last sample in the bunch), one 10-bit word, bunch length (i.e. the number of samples in the bunch, one 10-bit word), and sample amplitude values (few 10-bit words).

## IV Lossless compression of TPC signals

The lossless techniques of the data compression are based on the fact that TPC sample values (ADC and temporal) are not equally probable. A theoretical lower limit on the average word size using Huffman codding, or arithmetic coding lossless technique is given by entropy of the data source:

$`E(p)={\displaystyle p(A)\mathrm{log}_2p(A)}`$ (4)
The lossless techniques described in this paper are based mainly on an appropriate probability model for each data field of the ALTRO data format. Specific probability models for each sample in a bunch were developed. These models intend to capture both temporal correlation among samples and the characteristic shape of TPC electrical pulses.

### IV.1 Time information

As already mentioned, in the ALTRO data format time information is represented as the 10-bit number of the time-bin of the last sample of the bunch. The probability distribution of this variable is roughly uniform. In order to achieve better compression ratio this variable is substituted by the distance between two consecutive bunches. The probability of this variable is described by exponential distribution with much lower entropy factor. The entropy of temporal information is given by mean distance between two bunches. It depends on the event multiplicity, noise level and local occupancy, which is known function of the pad-row radius. In order to optimize entropy coding, it will be necessary to investigate probability distribution as a function of track multiplicity. This information will be known from other faster ALICE detectors.

The mean number of bits used for the coding of time information is roughly 4.9 bits for the full event with maximal track density. Using different codes in different places inside TPC, an additional 6% reduction in time information can be achieved.

### IV.2 Bunch length

In the ALTRO data format, the bunch length is represented as a 10-bit code number of samples in the bunch. The bunch length depends on the diffusion, the angular effect and the total deposited energy. There is no apparent correlation with data coded before. Small diffusion for short drift length is compensated by big angular effect. The total deposited energy is known only after coding of the bunch length. Since no apparent correlation with other data (e.g. length of adjacent bunches) exists and no better model (i.e. a model of events with lower entropy) could be found, this information is coded directly.

### IV.3 Sample values coding

Sample values are the main contribution to the resulting data volume. This subsection describes, first a basic model, and then introduces a more sophisticated one, that can provide higher performances in terms of compression efficiency.

Data compression can be obtained by directly applying entropy coding to the sample values without any modelling of the information source. This method will be referred bellow (in table 2) as Entropy Coding (EC).

### IV.4 Coding model based on the sample position

Improvements in compression performance can be obtained by appropriate modelling. A first improvement has been achieved by the fact that the statistics of the signal sample values depend on the position of the sample itself in the bunch.

Due to the pseudo Gaussian shape of most of the bunches, the first and the last sample of each bunch are likely to have a smaller value with respect to those in central positions. Similarly, small values are also expected for isolated samples, i.e. belonging to one-sample bunches (see table 1).

Therefore, a classification of the samples into three classes was chosen: one class for isolated samples, one for samples at the beginning and at the end of multiple sample bunches, and the last for samples in the central positions of a bunch. Using three different probability distributions for entropy coding the sample values can be coded more efficiently than using only one probability distribution. This coding scheme will be referred in table 2 as coding using Sample Position (SP).

### IV.5 Source models exploiting temporal correlation

Improvement on compression performances can be expected by exploiting temporal correlation, i.e. the correlation between consecutive samples; this can be done by implementing a suitable prediction scheme.

This approach is explained on the example, where a three-sample bunch is considered. Let us assume that the first two samples have already been coded and that the third one has to be coded. The code to be used for sample No. 3 may be chosen among eight possible codes according to the value of sample No. 2. In particular, this is done by subdividing the range of sample No. 2 (i.e. 0. . . 1023) into different intervals, and associating a different code (for the third sample) to each of these intervals.

This conditioned probability model can be extended to all the samples that are not in the first position in the bunch and for any bunch length. However, if the real-time implementation constraints are taken into account, and, in particular, the need to reduce the memory size of the model, it is not good to have an exceedingly large number of codes. Consequently, samples are partitioned into four classes only, to keep the complexity of the model low. This limitation does reduce the efficiency of the model but the reduction is only of the order of 0.6%. This coding scheme will be referred in table 2 as coding using Temporal Correlation (TC).

### IV.6 Comparison of different lossless technique

The results of different lossless approaches on simulated TPC data are shown in table 2. It may be noticed that the latter TC technique provides a compression of data down to 49.2% of the original size. Even this best technique provides reduction factor only by 3% better then direct EC technique.

Additional attempt tried to use predicted mean cluster shape information. Knowing the position of the bunch, the diffusion given by drift length ($`L_{\mathrm{drift}}`$) and inclination angle for primary particles are known. However, due to the fluctuation of cluster shape and due to the large amount of secondary particles with unknown angles, this prediction is not very good, and the entropy of the samples is reduced only by additional factor 2%.

### IV.7 Space correlation

In the trial to exploit space correlation, three lossless models have been considered. The first is based on spatially conditioned probability, the second on a predictive model, third on 2-dimensional cluster finder, with residual saving.

The first one is the equivalent, in the spatial domain, to what has been done for time correlation. Different codes are available to code the samples; for each sample, the appropriate code is selected according to the value of the samples in the same time-bin but in adjacent pads. This method provides poorer performance when compared with the one which exploits time correlation (the comparison being done using the same model complexity, i.e. number of probability distributions available in memory). Moreover, these two techniques cannot be easily combined, i.e. it is difficult to exploit both temporal and spatial correlations at the same time, because this would require a very large number of probability distributions (i.e. code tables).

The second method that has been investigated uses the prediction of the sample values from the samples in adjacent pads and coding the error of this prediction. Unfortunately, also for this model, the performance is not very good.

Pulses in one pad-row often resemble temporally shifted versions of those in the adjacent pad-row. The two methods described above have been modified by adding the first stage which shifts pulses so as to increase spatial correlation with adjacent. Although the performance has slightly improved, the increase of the compression efficiency was lower than expected. The correlations are relatively small. The main problems here are in the big amount of secondary particles crossing TPC with unknown $`\beta `$ angle (not pointing to the primary vertex), big spread of the particle momenta (unknown $`\alpha `$ angle) and the Landau fluctuation of deposited energy on different pad rows, which is almost uncorrelated. Moreover, the position of the original track relative to the pad, affects the correlation by a large factor. The signal amplitude in adjacent pads and adjacent pad-rows are very weakly correlated, unless the position and direction of the track is known.

In order to get better knowledge of the track position, two-dimensional cluster finding can be done before. The entropy of the stored residuals is by 30% lower than entropy of the original samples but there are problems with the track overlaps and with description of the cluster topology (i.e. where to store residuals).

Based on these results we conclude that it is not simple to exploit spatial correlation (i.e. correlations between adjacent channels). There might be more sophisticated and complex lossless models able to exploit it, but relatively simple models seem to fail.

## V Lossy compression of TPC signals

### V.1 Fluctuation and accuracy of the amplitude measurement

The number of primary ionization electrons produced by the charged particle in the gas is the random variable described by Poisson distribution with the mean value 14.35 $`\mathrm{cm}^1`$ for minimum ionizing particle in the gas of Alice TPC. The secondary electron production (described by $`E^{2.2}`$ probability distribution) increases the number of produced electrons. Maximum probable value is 25 $`\mathrm{cm}^1`$ of total electrons. This effect also smears probability distribution to the relatively broader Landau distribution.

Due to the angular effect and diffusion, electrons are distributed among several time-bins and pads. The number of electrons which contributes to the given pad and time-bin is described roughly by Poisson distribution. Each of the registered electrons is subject of gas multiplication which is described by exponential probability distribution. Over this, additional electronic noise is superimposed to each signal.

If we fix the track position and the number of primary electrons, the remaining sample uncertainties can be in the first approximation estimated as:

$`\sigma _\mathrm{S}=\sqrt{\sigma _{\mathrm{noise}}^2+G\times A}`$ (5)
where $`\sigma _{\mathrm{noise}}`$ is given by electronic noise and sampling imprecision, and G is the gain conversion factor.

The situation is more complicated, data samples are correlated through the time response function and the pad response function. The relative correlation between the samples depends on the ratio of the width of the response functions to the width given by stochastic processes.

### V.2 Dynamic precision of the digitization

In the following study, dynamic precision of sample quantization was investigated. The quantization was chosen to correspond to the sample deviation, modifying formula (5) to:

$`\delta _\mathrm{d}=\sqrt{K_{\mathrm{off}}^2+K_{\mathrm{cor}}^2\times A}`$ (6)
where $`K_{\mathrm{off}}`$ and $`K_{\mathrm{cor}}`$ factors were chosen as free parameters. $`K_{\mathrm{off}}`$ is proportional to the electronic noise and $`K_{\mathrm{cor}}`$ is given by statistics of the stochastic processes and by correlations. Different combinations of these factors were investigated.

In table 3 the influence of different quantization on the precision of the cluster characteristic determination is shown.

The gain factor $`G=A_\mathrm{t}/N_{\mathrm{el}}`$ ($`A_\mathrm{t}`$ is the total charge in cluster, $`N_{\mathrm{el}}`$ is the number of electrons contributing to the cluster) measures the precision of the local deposited energy determination. This factor is important for d$`E`$d$`x`$ measurement and consequently for particle identification (PID). The influence of the compression on the cluster position determination varies between 0 to 4%, depending on the compression factor, as can be expected. The shape of the cluster ($`\sigma _{\mathrm{PRF}}`$ and $`\sigma _{\mathrm{TRF}}`$), important for cluster quality determination, varies between 1 to 6%.

In table 4 the influence of the compression on the tracking is shown. Reported distortions in $`p_\mathrm{t}`$ and angular resolution are slightly smaller than in the case of the cluster position(0% up to 3.5%) . This is due to the other stochastic processes which contribute to the track parameters e.g. the multiple scattering. For high-momentum particles, where the influence of multiple scattering is not so important, the expected distortion will be determined by the cluster position distortions.

Reducing the number of the possible sample values, vector quantization of bunches were also investigated. Additional reduction factor of $``$6% was achieved on top of the results reported in table 3.

## VI Conlusions

Several methods of lossless TPC data compression was investigated (sec. IV.6). The best one dimensional methods provide compression factor down to 49.2%.

A lossy compression approach for the data generated by the TPC chamber in the ALICE experiment has been also investigated. The main idea was to preserve, on the level intrinsic noise, the three more important local quantities: the cluster position, the deposited energy and the shape of the cluster. Keeping the distortions of the local quantities at a reasonable level, the impact on the most interesting physical quantities, particle momenta and d$`E`$d$`x`$ is minimal. This approach achieves compression rates in the range from 35% down to 30%, depending on the desired precision. In this study we have focused on methods which are easy to implement in the frontend TPC electronics.