text:
# Structure optimization in an off-lattice protein model

## Abstract

We study an off-lattice protein toy model with two species of monomers interacting through modified Lennard-Jones interactions. Low energy configurations are optimized using the pruned-enriched-Rosenbluth method (PERM), hitherto employed to native state searches only for off lattice models. For 2 dimensions we found states with lower energy than previously proposed putative ground states, for all chain lengths $`13`$. This indicates that PERM has the potential to produce native states also for more realistic protein models. For $`d=3`$, where no published ground states exist, we present some putative lowest energy states for future comparison with other methods.

Predicting the structure of a protein, given its sequence of amino acids, is one of the central problems in computational biology. Since the problem is too difficult to be approached with fully realistic potentials derived from first principles, many authors have studied it in various degrees of simplifications. This involves in particular neglect of solvent water, simplifying the interactions, lumping together smaller groups of atoms, and putting everything on a discrete lattice. Among the most radically simplified models is the HP model of Dill and coworkers dill  where each amino acid is treated as a point particle on a regular (quadratic or cubic) lattice, and only two types of amino acids – hydrophobic (H) and polar (P) – are considered. Apart from the forces responsible for the connectedness of the chain, the only forces are contact forces between nearest lattice neighbours which are different for HH, HP, and PP pairs.

Even in this highly simplified model it is far from trivial to predict the native state for a given amino acid sequence yd93\_95 ; bd96 ; chikenji ; lw01 ; hmng03 ; backofen . The most efficient algorithms are either deterministic and cannot be generalized to more realistic models at all backofen , or use sequential importance sampling with re-sampling in the form of the pruned-enriched-Rosenbluth method (PERM) hmng03 . Although it was shown that the latter can be applied also to off-lattice homopolymers at higher temperatures g97 , it is not obvious that it will be efficient for off-lattice heteropolymers at the low temperatures needed for protein folding.

While there are a large number of benchmark cases for lattice protein models in the literature, there exist very few simple off-lattice models with known lowest energy states that can be used as benchmarks for efficient algorithms. One such model is the so-called AB model by Stillinger et al shh93 ; sh95  which also uses only two types of monomers, called now “A” (hydrophobic) and “B”(polar). The distances between consecutive monomers along the chain are held fixed to $`b=1`$, while non-consecutive monomers interact through a modified Lennard-Jones potential. In addition, there is an energy contribution from each angle $`\theta _i`$ between successive bonds. More precisely, the total energy for a $`N`$ monomer chain is expressed as

$$E=\underset{i=2}{\overset{N1}{}}E_1(\theta _i)+\underset{i=1}{\overset{N2}{}}\underset{j=i+2}{\overset{N}{}}E_2(r_{ij},\zeta _i,\zeta _j),$$
(1)
where

$$E_1(\theta _i)=\frac{1}{4}(1\mathrm{cos}\theta _i),$$
(2)
$$E_2(r_{ij},\zeta _i,\zeta _j)=4(r_{ij}^{12}C(\zeta _i,\zeta _j)r_{ij}^6).$$
(3)
Here $`r_{ij}`$ is the distance between monomers $`i`$ and $`j`$ (with $`i<j`$). Each $`\zeta _i`$ is either $`A`$ or $`B`$, and $`C(\zeta _i,\zeta _j)`$ is $`+1,+\frac{1}{2},`$ and $`\frac{1}{2}`$ respectively, for $`AA,BB,`$ and $`AB`$ pairs, giving strong attraction between $`AA`$ pairs, weak attraction between $`BB`$ pairs, and weak repulsion between $`A`$ and $`B`$.

This model has been studied in several papers shh93 ; sh95 ; ipp ; ipp3d ; gorse02  For its 2-d version, putative ground states for various $`AB`$ sequences and for various chain lengths are given in shh93 ; sh95 ; gorse02 . Similar models were also studied in ipp ; ipp3d ; tlp01 ; gorse01 , but putative ground states for these generalizations were not given at all or for very short chains only. The methods used to find low energy states of the AB model include neural networks shh93 , conventional Metropolis type Monte Carlo procedures sh95 , simulated tempering ipp3d , multicanonical Monte Carlo ipp , biologically motivated methods gorse01 ; gorse02 , and molecular dynamics tlp01 . In all cases the stochastic minimization can only lead to some state in the neighbourhood of a local (and hopefully also global) minimum. A greedy deterministic method such as conjugate gradient descent is subsequently applied to reach the minimum itself.

It is the purpose of the present paper to see whether PERM can be efficient for energy minimization in the AB model. In particular we shall use the new variant of PERM presented in hmng03 . We shall restrict ourselves to the subclass of “Fibonacci sequences” studied also in sh95 , defined recursively by

$$S_0=A,S_1=B,S_{i+1}=S_{i1}S_i$$
(4)
. Here * is the concatenation operator. The first few sequences are $`S_2=AB,S_3=BAB,S_4=ABBAB`$, etc. They have lengths given by $`N_{i+1}=N_{i1}+N_i`$, i.e given by the Fibonacci numbers. Hydrophobic residues ($`A`$) occur isolated along the chain, while $`B`$s occur either isolated or in pairs. The fraction of $`B`$s tends to the golden mean $`\gamma =0.618033`$ as the length $`N\mathrm{}`$.

Although PERM gives also detailed information about excited states and thermodynamic behaviour at temperatures $`T>0`$, we shall not discuss this here. For studying the dynamics of the folding transition, in contrast, we would have to assume some realistic microscopic dynamics. Just like other advanced sampling methods such as simulated annealing or parallel tempering, PERM sacrifices the realism of the dynamics for efficiency. In addition, as in the studies mentioned above, PERM is only used for coming close to the native state, and conjugate gradient descent is then used to reach the minimum energy state itself.

PERM is a biased chain growth algorithm with “population control”, i.e. a sequential importance sampling method with re-sampling liu , implemented recursively in a depth first fashion g97 . While chains grow, they acquire weights that include both Boltzmann factors and bias correction (“Rosenbluth” rr55 ) factors. During the growth, samples with large weight are cloned, while chains with too small weight are pruned out. Except for the depth-first implementation and for the fact that it gives the correct Gibbs-Boltzmann statistics, PERM resembles therefore genetic algorithms. While the original version of PERM was quite successful for lattice proteins g98 ; bast98  and for a host of other applications perm-review , it worked rather poorly for minimization of off-lattice polymer models fg99 .

In this paper we therefore employ an improved version called nPERMis in hmng03  (for ”new PERM with importance sampling”). Basically, instead of making exact clones of high weight chains and hoping that these clones will evolve differently during the subsequent growth (as in original PERM), we now branch such that the last monomers are different at the point of branching. Thereby we now force the two copies to be distinct, and we avoid the loss of diversity that also plagues genetic algorithms when the evolution pressure is too high.

For lattice polymers hmng03  one has for each partially grown chain a finite number of “candidate directions” for the next step. One first estimates the total weight of all these one-step continuations. Based on this estimate, one decides on the number of clones to be made. If, say, one wants to make $`k`$ clones, one scans all possible $`k`$-tuples of possible different candidate directions, and selects one of these tuples according to its weight. For off-lattice polymers one proceeds exactly in the same way, with one exception: The candidates are now no longer the lattice bonds, but one has to choose $`K`$ candidate locations for the next monomer randomly. The number $`K`$ is an important parameter. While $`K5`$ was optimal for 3-d polymers near the $`\mathrm{\Theta }`$-point g97 , we found that lowest energies were reached in the AB model for $`K50`$. While it was necessary to make clonings with very many siblings in simulating the HP model with the old version of PERM g98 ; bast98 , we now obtained good results by restricting ourselves to $`k`$-tuples with $`k3`$.

Another important parameter is the temperature at which the simulation is run. We typically used temperatures well below the collapse transition, $`kT0.1`$ or even lower. In order to speed up the ground state search, we also modified the Lennard-Jones potential by putting $`E_2(r)=+\mathrm{}`$ for $`r<1`$. This hard core constraint reduces the available phase space, but has no effect on ground state configurations (we did not use it in the conjugate gradient minimization, and we checked that it was satisfied after minimization). For chain deformation algorithms it could slow down the dynamics, since the hard cores could act as barriers, but it can only improve any pure chain growth algorithm. Finally, as a last trick, we used equally spaced azimuthal angles for all candidates (with one overall angle chosen at random, for each group of candidates), in order to make them cover the unit sphere more uniformly. All simulations were done on Linux and UNIX workstations. CPU times were up to 2 days, but their precise values are not very significant. Exact timings would involve frequent comparisons of the minimizer basins of attraction reached by PERM, which we considered as too time consuming.

In Table I we list the lowest energies thus obtained for the two- and three-dimensional $`AB`$ model for all Fibonacci sequences with $`13N55`$. The latter is equal to the length of the longest sequence studied in sh95 . Let us first discuss the case $`d=2`$. For comparison we quote also the putative ground state energies from Table II of sh95 . For $`N<13`$, our energies agree perfectly with those of sh95 . Except for the shortest chain with $`N=13`$, already PERM gave in all cases shown in Table I lower energies than those found in sh95 . In all these cases already PERM by itself showed that the topologies shown in sh95  are not the native ones. While the subsequent gradient descent improved the energies substantially, it changed in no case the overall topology.

The latter is true also for $`d=3`$, although there the subsequent minimization gave even larger energy changes than in $`d=2`$. This shows that in $`d=3`$, too, PERM is able to find states very close to the native ones. Since there exist no published ground state energies for the $`3d`$ AB model, we are unable to compare PERM with other methods.

The configurations corresponding to the energies shown in Table I are shown in Figs. 1 (for $`d=2`$) and 2 (for $`d=3`$). For $`d=2`$ we see that none of the configurations, except the one for $`N=13`$, have single hydrophobic cores. Instead, the hydrophobic (A) monomers form clusters of typically 4 to 5 particles. This is easily explained by the fact that hydrophobic monomers always are flanked by polar monomers along the chain. Thus a clean separation into hydrophobic and polar regions is impossible. This shows that the AB model with Fibonacci sequences would be a very poor model for real proteins in $`d=2`$. From Fig. 3 we see that the same is true to a lesser degree in $`d=3`$. There the chains with $`N=21`$ and $`N=34`$ fold into configurations with single hydrophobic cores (except for a single A monomer which keeps out in both cases), and only the chain with $`N=55`$ forms 2 clearly disjoint main hydrophobic groups.

In conclusion we have extended the PERM algorithm to an off-lattice two-species protein model. We have shown that it performs well, indeed we are able to refute with it previous claims for putative ground states.

The chosen model is not very realistic. Partly this follows from the restriction to two types of monomers, partly from the fact that we did not include, as in ipp3d , more realistic local (bond angle and torsion) forces, and partly from the restriction to Fibonacci sequences. Each of this features could have been easily avoided, and PERM works indeed equally well if we modify any of them. But it was not our aim to present a realistic model. Rather we wanted to treat a model which is suitable for benchmarking, because of it is defined in a simple way and because it was already studied in detail before.

It is less obvious whether PERM would also perform well for all-atom models with realistic potentials, or even with explicit solvents. Typically, its performance decreases quite rapidly with the number of degrees of freedom, but presumably it shares this with other modern methods like multicanonical sampling and parallel tempering. To answer this question, we have started such simulations with the ECEPP force field implemented in SMMP smmp . But it is still too early to draw any conclusions.

Acknowledgments: We thank Walter Nadler for numerous very fruitful discussions and for critically reading the manuscript.


format: text

text:
# Radiative transfer in 3D model stellar atmospheres

## 1. Introduction

Reality is (at least) 3-dimensional. This obvious statement is equally applicable to stellar atmospheres, in particular those of late-type stars where the surface convection zone reaches the region responsible for the emitted stellar radiation. Even a casual glimpse of the surface of our own Sun immediately reveals the existence of photospheric inhomogeneities in the form of granulation, in addition to magnetic phenomena like sunspots and flux tubes which probably at some level also are related to the convective motions. Granulation is the direct observational manifestation of convection: concentrated, rapid downdrafts of cool material (intergranular lanes) in the midst of broad, slow upwellings of warm gas (granules). Furthermore, the solar atmosphere is distinctly time-dependent: the individual granules are continuously evolving with typical life-times of a few minutes and waves excited by the convective motions are regularly propagating into the upper layers of the photosphere and the chromosphere. Similar features are expected on other late-type stars and may in fact often be more severe than in the Sun. In view of the dynamic nature of stellar atmospheres, one has to worry whether the use of time-independent hydrostatic 1D model atmospheres as regularly done in for example stellar abundance analyses could lead to severe systematic errors. Here we will describe efforts currently being made in developing more physically motivated models of stellar atmospheres and the line formation processes, which should hopefully place the findings on firmer footing.

With the advent of supercomputers it is now feasible to perform highly realistic time-dependent 3D hydrodynamical simulations of stellar surface convection (e.g. Nordlund & Dravins 1990; Stein & Nordlund 1998; Asplund et al. 1999, 2000a,b,c; Asplund & García Pérez 2001). Without invoking any of the free parameters traditionally employed in 1D analyses the 3D simulations successfully reproduce a large number of observational constraints such as the granulation topology and detailed line shapes and asymmetries. In this context the energy exchange between the gas and the radiation field is absolutely crucial by largely determining the atmospheric structure and ultimately driving the convection. In the present contribution, special emphasis will be placed on how realistic radiative transfer is achieved both in constructing the 3D convection simulations and in predicting the emergent spectra since different approaches, assumptions and approximations are necessary in the different cases. It is important to emphasize that although the present 3D models are not perfect, the main limitations are of a computational nature, which will be naturally addressed in the future when more powerful computers become available, and due to a lack of manpower. In contrast, 1D models are flawed in a physical sense, since the intrinsic restriction to 1D can not adequately describe all aspects of what is inherently a 3D phenomenon. That is of course not to say that 1D models do not have their merits but it is important to investigate their potential shortcomings, which can only be achieved by comparisons with more sophisticated modelling.

## 2. 3D Hydrodynamical Simulations of Stellar Surface Convection

### 2.1. Ingredients

The 3D hydrodynamical simulations of stellar surface convection described herein have been computed with a 3D, time-dependent, compressible, explicit, radiative-hydrodynamics code (Stein & Nordlund 1998). The hydrodynamical equations for conservation of mass, momentum and energy are solved on a non-staggered Eulerian mesh with gridsizes of $`100^3`$. The physical dimensions of the grids are sufficiently large to cover many ($`>10`$) granules simultaneously in the horizontal direction and many ($`>10`$) pressure scale-heights in the vertical. The simulations are thus only covering a small but representative fraction of the whole stellar surface and convection zone, which gives statistically significant results for long time sequences. In terms of continuum optical depth the simulations extend at least up to log $`\tau _{\mathrm{Ross}}5`$ which for most purposes are sufficient to avoid numerical artifacts of the open upper boundary on spectral line formation which could otherwise skew the results. The lower boundary is located at large depths to ensure that the inflowing gas is isentropic and featureless, while periodic horizontal boundary conditions are employed. The temporal evolution of the simulations cover several convective turn-over time-scales to allow thermal relaxation to be established. The simulations include a hyper-viscosity for stabilization purposes whose parameters are fixed from standard hydrodynamical test cases. No magnetic fields have been included although the code is capable of simulating magneto-convection. The simulations described here thus only apply to magnetically inactive stars such as the Sun.

Since the intention is to make the simulations as realistic as possible, state-of-the-art input physics is used. The adopted equation-of-state comes from Mihalas et al. (1988), which accounts for ionization, excitation and dissociation of the most important atoms and molecules. The Uppsala opacity package (Gustafsson et al. 1975 with subsequent updates) provide the continuous opacity sources and the extensive lists of Kurucz (1998, private communication) the line opacities. The input parameters characterizing the models are the surface gravity (log $`g`$), the overall metallicity (\[Fe/H\]) and the entropy at the lower boundary. The effective temperature ($`T_{\mathrm{eff}}`$) is a consequence of the granulation properties, which depend on the entropy structure. To construct a model with a prespecified $`T_{\mathrm{eff}}`$ is therefore a rather time-consuming task involving making several simulations until exactly the right entropy at the bottom has been found which produces the targeted $`T_{\mathrm{eff}}`$. An alternative approach is to settle for a $`T_{\mathrm{eff}}`$ in the vicinity of the desired one and carry out a differential analysis with 1D models with exactly the same stellar parameters. This procedure is justified since 3D effects will not change drastically with minor adjustments of the parameters and has been adopted in most cases by us. The reader is referred to Stein & Nordlund (1998), Asplund et al. (1999, 2000a,b) and Asplund & García Pérez (2001) for further details of the 3D convection simulations employed here.

### 2.2. Radiative Transfer

An accurate representation of the energy balance between the gas and the radiation field is crucial in order to achieve a realistic atmospheric structure. We therefore compute detailed 3D radiative transfer at each time-step for a selected number of inclined rays (typically around 10) using a Feautrier-like (1964) long characteristic formal solution. To improve the numerical accuracy, the radiative transfer is solved for a finer depth grid and the results subsequently interpolated back to the original hydrodynamical grid. For computational reasons, the assumption of local thermodynamic equilibrium (LTE) without continuous scattering terms ($`S_\lambda =B_\lambda `$) is made throughout. The effects of line-blanketing is accounted for through opacity binning (Nordlund 1982) using normally four opacity bins. The original 2748 wavelength points are sorted into groups representing continuum, weak lines, intermediate strong lines and strong lines. For each bin pseudo-Planck functions are computed and stored in a lookup-table together with the opacities and equation-of-state. Even with the much simplifying assumptions of LTE and only four opacity bins, it is noteworthy that the most time-consuming task in the simulations is the solution of the radiative transfer while the time-integration of the hydrodynamical equations only require a fraction of the time for the highest resolution simulations. We are currently designing methods to switch to a selected opacity sampling technique with $`50100`$ wavelength points instead of opacity binning (Trampedach et al., these proceedings).

### 2.3. Stellar Surface Convection

A recent realisation which has emerged from convection simulations like those described here is the fundamental importance of radiative cooling in a thin surface layer for driving the convection (Stein & Nordlund 1998). As the ascending gas approaches the optical surface radiation starts to escape, lowering the temperature of the gas. In these regions the opacity is very sensitive to the temperature, which means that the opacity decreases rapidly. As a consequence more photons can leak out bringing the temperature down further. This strong feedback causes a sudden drop in the gas temperature, bringing the gas from around $`\mathrm{10\hspace{0.17em}000}`$ K to around $`\mathrm{5\hspace{0.17em}000}`$ K in about a minute for the Sun. The cooled, entropy-deficient material now has a negative buoyancy and is therefore being accelerated downward. It is noteworthy that almost all buoyancy work is done by the downflows rather than the upflows, which instead are mainly being pushed upwards by mass conservation. This picture is clearly conceptually very different from the traditional view of convection as rising bubbles, as in the mixing length theory or relatives thereof.

Granulation in other solar-type stars is qualitatively similar to the solar case with broad upflows and concentrated downdrafts according to numerical simulations. This is also the case at low metallicities, but there are distinct differences in the optically thin layers. The upflowing gas overshoots into the convectively stable layers above the optical surface. The gas temperature is then determined by a competition between radiative heating and adiabatic cooling due to the expansion following the density stratification. At solar metallicity, the large number of spectral lines (continuum absorption is very inefficient in these optically thin layers) retain the temperature close to the radiative equilibrium value. At low metallicity however, the severe lack of lines shifts the balance towards adiabatic cooling, leading to distinct sub-radiative equilibrium temperatures (Asplund et al. 1999). The temperature difference between a 3D model and a corresponding 1D model can easily reach $`\mathrm{1\hspace{0.17em}000}`$ K in the upper layers, which naturally have dramatic consequences for any spectral lines formed in those regions, such as molecular lines.

## 3. 3D LTE Spectral Line Formation

### 3.1. Radiative Transfer

Using the 3D convection simulations as a time-series of 3D model atmospheres, spectral line formation calculations can be performed which subsequently can be spatially and temporally averaged to produce disk-integrated line profiles. In LTE, the solution of the radiative transfer equation is straightforward as the level populations are directly obtained from the Boltzmann and Saha distributions. This formal solution is computed using a Feautrier (1964) scheme with long characteristics without any continuum scattering terms but accounting for the Doppler shifts caused by the convective motions in the simulations. Before the radiative transfer calculations, the vertical resolution is improved but the horizontal resolution of the original simulation is decreased to typically $`N_\mathrm{x}N_\mathrm{y}=5050`$ to ease the computational burden, which however does not encumber the results, as extensive tests show. The calculation of a 3D flux profile require some 100 wavelength points ($`N_\lambda `$), 10-20 angles ($`N_{\mathrm{angles}}`$) and at least 50 snapshots ($`N_\mathrm{t}`$) to obtain statistically significant profiles in terms of line asymmetries; for most abundance purposes significantly fewer snapshots are necessary. Thus a single 3D LTE flux profile require typically $`N_\mathrm{t}N_\mathrm{x}N_\mathrm{y}N_{\mathrm{angles}}N_\lambda >10^8`$ 1D radiative transfer calculations, which however is still achievable on modern workstations. In addition, the calculations are normally performed for three different element abundances to enable interpolation to the requested line strength. For computational reasons, the monochromatic continuous opacities and equation-of-state data are pretabulated and interpolated during the line calculations. For each snapshot the spatially resolved intensity profiles (Fig. 2) for the different angles are stored, from which a disk-integration can be performed to account for rotational broadening (Dravins & Nordlund 1990) in the final spatially and temporally averaged flux profile. It is important to emphasize that no micro- or macroturbulence enter the calculations.

### 3.2. Spectral Line Shapes and Asymmetries

The 3D hydrodynamical models described above may be more sophisticated than existing 1D hydrostatic models but they are not necessarily more realistic. Special efforts have therefore been made in comparing in great detail the predictions with observations in terms of for example line shapes and asymmetries to verify how reliable the simulations are. Considering the amazing assortment of strengths, shapes and shifts of spatially resolved profiles, it is obvious that the final averaged profile will be an excellent probe of the atmospheric structure and a very challenging test of the models. Here, the Sun is a perfect test-bench as comparisons can be made without the interference of rotational broadening and information is available for spatially resolved profiles.

As has been well-known since long, 1D models can neither predict the strengths and shapes of lines without invoking the ad-hoc fudge factors micro- and macroturbulence and of course fail utterly to describe line asymmetries. In stark contrast, the 3D line profiles agree almost perfectly with observations, as shown in Figs. 3 and 4 (see also Asplund et al. 2000b, Nissen et al. 2000; Allende Prieto et al. 2002). The self-consistently calculated Doppler shifts arising from the convective motions and to a minor part the oscillatory motions broaden the lines sufficiently. This is also evident from tests when artificially removing all velocity information in the simulations (Fig. 3), giving profiles which closely resemble 1D profiles without micro- and macroturbulence. Thus, the concepts of micro- and macroturbulence introduced in 1D analyses have nothing to with neither turbulence nor microscopic motions but stems from gradients in the convective velocities which are resolved at the present generation of 3D model atmospheres (Asplund et al. 2000a). This is partly due to the fact that the line formation process is heavily biased towards upflows (large area coverage, high continuum intensity, steep temperature gradients) where the divergent nature of the flow prevents large amount of turbulence to develop in spite of the high Reynolds-numbers of the plasma (Nordlund et al. 1997). The excellent overall agreement between predicted and observed line profiles in general makes it very easy to detect minor blends which otherwise may well have gone unnoticed (Allende Prieto et al. 2001).

Even the detailed asymmetries and shifts of lines are very well reproduced with the 3D simulations, both for the Sun and for other solar-type stars (Asplund et al. 2000b; Allende Prieto et al. 2002), as illustrated in Fig. 4. The accuracy is typically $`0.050.1`$ km s-1 for weak lines, which is similar to the uncertainties in the atomic data. In fact, in a few cases the comparisons between predicted and observed solar line profiles have revealed minor typographical errors in tables of laboratory wavelengths. For a large sample of weak solar Fe i lines, Asplund et al. (2000b) found an average difference of $`0.000\pm 0.053`$ km s-1 for the observed and theoretical line shifts (with the gravitational redshift removed). Strong lines show a systematic discrepancy in the line core which is likely attributable to effects of the artificial upper boundary or a less realistic temperature structure in the uppermost atmospheric layers. Also the observed behaviour of spatially resolved solar line profiles are well matched by the predictions from the simulations (Kiselman & Asplund, in preparation).

## 4. 3D Non-LTE Spectral Line Formation

### 4.1. Radiative Transfer

Line formation does not of course obey LTE in general. To make optimal use of the 3D hydrodynamical model atmospheres described above it is necessary to develop tools that can also handle detailed 3D non-LTE line formation. This challenging task has only recently been addressed with a very promising outlook for the future (Kiselman 1997; Uitenbroek 1998; Botnen & Carlsson 1999; Asplund et al. 2002). The major complication in any non-LTE calculation is of course that the level populations depend on the radiation field which in turn depends on the level populations. The radiative transfer equation must therefore be solved for all relevant wavelengths together with the rate equations of the various levels, normally under the simplifying assumption of time-independence (so-called statistical equilibrium, $`\mathrm{d}n_\mathrm{i}/\mathrm{d}t=0`$). To obtain a consistent result, an iterative procedure is adopted using normally an accelerated lambda-iteration technique allowing convergence typically within 10 iterations ($`N_{\mathrm{iter}}`$). For a normal size model atom there may 100 radiative transitions ($`N_{\mathrm{tran}}`$), each described with 50-100 wavelengths ($`N_\lambda `$). In total, each 3D non-LTE calculations therefore require $`N_\mathrm{t}N_\mathrm{x}N_\mathrm{y}N_{\mathrm{angles}}N_{\mathrm{tran}}N_\lambda N_{\mathrm{iter}}`$ corresponding 1D radiative transfer solutions. Hence it should come as no surprise that the formal solution of the radiative transfer equation with the current estimate of the source function is by far the most dominating task in the calculations. It is not possible to perform such 3D non-LTE calculations for long time-series as can be done in LTE and instead the computations are restricted to a few individual snapshots. Fortunately, 3D non-LTE effects are quite robust entities which differ little from snapshot to snapshot, as verified by various test calculations.

The 3D non-LTE code multi3d (Botnen & Carlsson 1999; Asplund et al. 2002) is essentially a 3D version of the widely used multi-code for 1D statistical equilibrium calculations (Carlsson 1986). The formal solution of the radiative transfer equation is computed using a short-characteristic method utilizing the horizontal boundary conditions of the 3D atmospheres and accounting for Doppler shifts introduced by the convective motions. The equation-of-state and continuous opacities are provided by the Uppsala package (Gustafsson et al. 1975 and subsequent updates) and allowance for line-blanketing in the photo-ionization rates is possible.

### 4.2. Li Abundances in Metal-Poor Stars

Given the possibility to estimate the total baryon density of the Universe using the Li abundances in pristine metal-poor halo stars, a large effort has gone into this endeavour ever since the first detection of the Li i 670.8 nm resonance line in such stars (Spite & Spite 1982; Ryan et al. 1999). In view of the large differences in the temperature structures of 3D and 1D models of halo stars (Asplund et al. 1999, see also Sect. 2.3.), there is clearly a great need to study the Li line formation in detail in 3D model atmospheres. Indeed, preliminary calculations assuming LTE indicated that existing 1D analyses had overestimated the primordial Li abundance by as much as $`0.20.35`$ dex (Asplund et al. 1999). Extreme caution must however be exercised in interpreting these results as Li may be prone to severe over-ionization effects which could compensate the 3D LTE corrections. Indeed our new 3D non-LTE calculations reveal pronounced 3D non-LTE effects which bring the derived Li abundances almost exactly back to those determined with 1D models (Asplund et al. 2002).

With multi3d we have recently investigated possible departures from LTE for Li in the two halo stars HD 140283 and HD 84937 using a 21-level Li model atom with 90 radiative transitions (Carlsson et al. 1994). Colllisional excitation and ionization due to electrons are accounted for but not the corresponding case of hydrogen collisions, which however have limited impact on the final results. Compared with the LTE case, the Li line is much weaker in non-LTE, as shown in Figs. 5 and 6. In LTE the line is generally strong above the warm upflows while this is not the case in non-LTE, which implies that the population density is largely controlled by the radiation field rather than the local temperature. As anticipated, the main non-LTE mechanism is over-ionization, which decreases the line opacity significantly. In terms of abundances, the 3D non-LTE calculations yield Li abundances well within 0.05 dex to the standard 1D LTE case, at least for these two halo stars. It is noteworthy though that the subordinate 610.4 nm line show an abundance correction of about $`+0.15`$ dex compared with 1D, which should be detectable with high enough $`S/N`$ observations. It would be premature however from these limited calculations to conclude that the apparent trend in Li abundances with metallicity derived with 1D models (Ryan et al. 1999) is due to galactic chemical evolution of Li. The exact magnitude of the 3D non-LTE effect will depend on the stellar parameters, in particular $`T_{\mathrm{eff}}`$ and \[Fe/H\]. We are planning to investigate the metallicity dependence of the 3D corrections in the near future.

## 5. Concluding Remarks

With realistic 3D hydrodynamical model atmospheres of late-type stars and corresponding 3D LTE and non-LTE spectral line formation tools now becoming available many potential systematic errors in standard analyses can finally be removed eliminated. In many cases the results may well be dramatically revised (e.g. Asplund & García Pérez 2001; Allende Prieto et al. 2001; Nissen et al. 2002), which can have large impact on the inferred conclusions on stellar, galactic and cosmic evolution. However, the significant efforts in pursuing this ambitious endeavour are justified also for those cases where 3D only confirms previous 1D results as the uncertainties are greatly reduced. The reliability of existing results can clearly only be judged after more sophisticated calculations have been performed.

#### Acknowledgments.

We are greatly indebted to the continuing efforts of a large number of wonderful collaborators in the quest of developing 3D model atmospheres and 3D line formation and their applications to astrophysical problems. In particular we would like to mention Carlos Allende Prieto, Mats Carlsson, Ana García Pérez, Nicolas Grevesse, Dan Kiselman, David Lambert, Poul Erik Nissen, Åke Nordlund, Francesca Primas, Jacques Sauval, Bob Stein, and Regner Trampedach. Finally, we are grateful to the editors for their patience.

## References

Allende Prieto, C., Lambert, D.L., Asplund, M. 2001, ApJ, 556, L65

Allende Prieto, C., Asplund, M., García López, R.J., Lambert, D.L. 2002, ApJ, 567, 544

Asplund, M., García Pérez, A.E. 2001, A&A, 372, 601

Asplund, M., Carlsson, M., Botnen, A.V. 2002, submitted to A&A

Asplund, M., Nordlund, Å., Trampedach, R., Stein, R.F. 1999, A&A, 346, L17

Asplund, M., Ludwig, H.-G., Nordlund, Å., Stein, R.F. 2000a, A&A, 359, 669

Asplund, M., Nordlund, Å., Trampedach, R., et al. 2000b, A&A, 359, 729

Asplund, M., Nordlund, Å., Trampedach, R., Stein, R.F. 2000c, A&A, 359, 743

Botnen, A.V., & Carlsson, M. 1999, in Numerical astrophysics, ed. S.M. Miyama et al., 379

Carlsson, M. 1986, Uppsala Astronomical Report No. 33

Carlsson, M., Rutten, R.J., Bruls, J.H.M.J., & Shchukina, N.G. 1994, A&A, 288, 860

Dravins, D., Nordlund, Å. 1990, A&A, 228, 203

Feautrier, P. 1964, C. R. Acad. Sci. Paris, 258, 3189

Gustafsson, B., Bell, R.A., Eriksson, K., Nordlund, Å. 1975, A&A, 42, 407

Kiselman, D. 1997, ApJ, 489, L107

Mihalas, D., Däppen, W., Hummer, D.G. 1988, ApJ, 331, 815

Nissen, P.E., Asplund, M., Hill, V., D’Odorico, S. 2000, A&A, 357, L49

Nissen, P.E., Primas, F., Asplund, M., Lambert, D.L. 2002, A&A, 390, 235

Nordlund, Å. 1982, A&A, 107, 1

Nordlund, Å., Dravins, D. 1990, A&A, 228, 155

Nordlund, Å., Spruit, H.C., Ludwig, H.-G., Trampedach, R. 1997, A&A, 328, 229

Ryan, S.G., Norris, J.E., & Beers, T.C. 1999, ApJ, 523, 65

Spite, F., & Spite, M. 1982, A&A, 115, 357

Stein, R.F., Nordlund, Å. 1998, ApJ, 499, 914

Uitenbroek, H. 1998, ApJ, 498, 427


text:
# 1 Introduction

## 1 Introduction

There exists a widespread understanding in scientific community that the moment of rupture in heterogeneous materials is somewhat similar to a critical point. The fracturing process is then a kind of a second-order phase transition . This underlies the intuitive assumption which compares global failure in disordered materials at mesoscopic scale with percolation at microscopic scale. So, as percolation, global failure is assumed to be a critical phenomenon. The fracturing process is then classified as a self-organised random irreversible process. Extensive numerical simulations and experimental measurements have confirmed this conclusion. This concerns the global failure of materials as well as earthquakes , since an earthquake can also be considered as a global failure of a large material mass. The classical power law, defining the energy release rate, has been verified in the critical region, close to the time of either a large earthquake or a global materials failure. The related critical exponent is, in general, a complex number, which is due to taking into account long-range elastic interactions transported by stress fields around defects and cracks. The imaginary part of the exponent gives rise to the so-called log-periodic corrections which have been identified quite early in renormalization group solutions for critical phase transitions . Experimental results of acoustic emission measurements and seismograms unambiguously revealed the existence of oscillatory corrections before materials fractures and earthquakes, respectively . If one would know the general law describing this acoustic emission, one would be able to make an early prediction of the global failure in heterogeneous materials. It is a principal, not yet solved problem, to find such a general law, which would be valid not solely in the asymptotic vicinity of the critical point, but in the whole interval of pressures or times, both near the critical point of failure as well as far from it. And this is the main aim of the present paper to suggest a mathematically grounded derivation for such a law and to confront it to available experimental data . For this purpose, we employ the recent progress in the self-similar approximation theory \[12–18\] which has been successfully applied to the theory of critical phenomena and to time series forecasting \[19–21\]. We would like to stress here that this theory is general and can be used for any given series of experimental data, whether this concerns phase transitions, financial markets, or any other series.

## 2 Modelling

Let us be interested in the behaviour of the energy release $`E(p)`$ as a function of pressure $`p`$, considered between $`p0`$, i.e. at the initial stage when the very early damages occur, and up to the point $`p=p_c`$ of global failure. Note that $`p`$ refers to fracturing under applied spatially uniform loading, as in experiments where water pressure inside a spherical tank increases up to the material failure. Experimentally, the applied pressure is usually increased linearly with time. The behaviour of $`E(p)`$ at the vicinity of the critical point of rupture is

$$E(p)E(p_c)+A(p_cp)^\alpha +B(p_cp)^\alpha \mathrm{cos}[\omega \mathrm{ln}(p_cp)+\phi ],$$
(1)
as $`pp_c`$, where $`A,B,\alpha ,\omega `$ and $`\phi `$ are parameters ($`\alpha `$ and $`\omega `$ are the real and imaginary parts, respectively, of the so-called critical complex exponent). The problem is how to find the behaviour of $`E(p)`$ in the whole range $`0pp_c`$? Here we consider the case of one variable, pressure. A more general loading, involving not one variable, like $`p`$, but several variables can also be treated by the theory. However, before being involved in such generalisations, we would like to show on a simpler example how to accomplish the extrapolation. Our main aim in this letter is to demonstrate the general possibility of realising such an accurate extrapolation for the energy release.

Let us transform Eq. (1) by introducing the dimensionless variable $`x(p_cp)/p_c`$, $`(0x1)`$ and the reduced energy release $`f(x)E(p)/E(p_c)`$ which is a dimensionless function. Then, expansion (1) takes the form

$$f(x)1+\stackrel{~}{a}(x)x^\alpha ,\stackrel{~}{a}(x)\lambda [1+\mu \mathrm{cos}(\omega \mathrm{ln}x+\beta )]$$
(2)
whose parameters $`\beta `$, $`\lambda `$, and $`\mu `$ can be easily expressed through $`p_c`$, $`E(p_c)`$, $`A`$, $`B`$, $`\omega `$, and $`\phi `$. Note that function $`\stackrel{~}{a}(x)`$ cannot be expanded in powers of $`x`$ near zero. Therefore, expansion (2) can be treated as a generalised asymptotic series, as defined by Poincare , with $`\stackrel{~}{a}(x)`$ considered as a coefficient. Employing the self-similar approximation theory \[12–21\], the asymptotic series (2) can be extrapolated to

$$f^{}(x)=\mathrm{exp}[c(x)x^\alpha ],c(x)=[1+\mu \mathrm{cos}(\omega \mathrm{ln}x+\beta )]\tau ,$$
(3)
where $`\tau `$ has to be defined from an optimization or boundary condition. As is clear, the boundary condition $`E(0)0`$ is not convenient here, leading to large errors in defining $`\tau `$ because of the physical impossibility to determine the very early precursors with a high accuracy. A more judicious choice is to consider a global physical quantity. In the same way as treating $`E(p)`$ as the integrated energy release rate, we may define $`F(p)`$ as the integral energy release $`_0^pE(q)𝑑q`$. Given $`p_0`$, the parameter $`\tau `$ then can be defined from the sum optimization rule

$$_{y_0=\frac{p_cp_0}{p_c}}^1f^{}(x)𝑑x=\frac{F(p_0)}{p_cE(p_c)}(y_00,p_0p_c0).$$
(4)
In general, the lower limit in the integral (4) can pertain to the whole interval . In reality, the trustful experimental data are available only below $`p_c`$. One has to choose the data at the highest available pressure, so that the lower limit in the sum rule (4) be sufficiently small.

In this way, the self-similar approximant (3) becomes

$$f^{}(x)=\mathrm{exp}\left\{\left[1+\mu \mathrm{cos}(\omega \mathrm{ln}x+\beta )\right]\tau x^\alpha \right\}$$
(5)
with $`\tau `$ given by the sum rule (4). Note that this formula contains the same number of parameters as the initial formula (1). Returning, with the help of relations (1) and (2), to the energy release, we come to the formula

$$E^{}(p)=E(p_c)f^{}(x),$$
(6)
which extrapolates $`E(p)`$ for the whole interval $`0pp_c`$. It is worth stressing that the law (6) is obtained as a self-similar extrapolation of the acoustic emission signals observed at the vicinity of the rupture. Another possibility would be to base this kind of extrapolation, starting from the signals existing at the initial stage of the process, being yet very far from the rupture. This latter way was used by Gluzman et al. , who started with a given polynomial expansion at the early time of acoustic emissions, treating the power law behaviour close to rupture as a boundary condition. To our mind, this opposite approach has two weak points: First of all, the early acoustical precursors, as is known, are very small, being embodied in the acoustical noise, because of which they can be hardly measured with a good accuracy. Second, it looks difficult or even impossible to extract information on log-periodic corrections from the early acoustic signals. Another phenomenological expression for the energy release in the noncritical region has been suggested by Sornette and Andersen  whose arguments were based on the existence of a scaling of the macroscopic elastic modulus and the elastic energy release rate in the thermal fuse model . The following phenomenological form has been proposed

$$E(p)E(p_c)A\left(\mathrm{tanh}\frac{p_cp}{\tau _0}\right)^\alpha +B\left(\mathrm{tanh}\frac{p_cp}{\tau _0}\right)^\alpha \mathrm{cos}\left\{\omega \left(\mathrm{tanh}\frac{p_cp}{\tau _o}\right)+\phi \right\},$$
(7)
which is a pure power law, like Eq. (1), in the critical region close to rupture and exponentially relaxes in the noncritical region far from rupture, where only a few damages occur. However, the suggested form (7) possesses an unphysical behaviour, being negative in a wide range of its variable. At the same time, our self-similar formula (5) yields the energy release (6) that is always positive. In Fig. 1, we compare the behaviour of $`E(p)`$ , in the dimensionless form, for both Eqs. (6) and (7), with the parameters chosen so that to have similar variations when approaching the critical point. In this figure, one can clearly see that Eq. (7) possesses unphysical negative values for the major part of the interval , while Eq. (6) is everywhere positive.

## 3 Comparison with experimental data

In order to fit our self-similar formula to experimental results of acoustic emission measurements, we refer to the most accurate data, available nowadays, that is those of Anifrani’s team  at Aerospatiale-Matra Inc., which are reported by Johansen and Sornette . In these measurements, the acoustic emission has been recorded during the pressurisation of spherical tanks of kevlar or carbon fibres pre-impregnated in a resin matrix wrapped up around a thin metallic liner (steel or titanium) fabricated and instrumented by Aerospatiale-Matra Inc. It has been found that the seven acoustic emission recordings of seven pressure tanks, which were brought to rupture, exhibit a clear acceleration in agreement with a power law divergence as expected by the critical point theory. A strong evidence of oscillatory corrections, forming the intermittent succession of bursts of acoustic emission, when approaching the rupture, has been clearly identified.

In Fig. 2, we present the energy release rate, which is the derivative of our formula (6) for the cumulative energy release. The result is obained by fitting to the experimental curve reported at the right bottom in figure 1 of the paper by Johansen and Sornette (2000). These experimental data explicitly show an oscillatory behaviour of the energy release rate before the rupture. The best fit for the five parameters $`\alpha `$, $`\beta `$, $`\mu `$, $`\tau `$, and $`\omega `$ have been obtained for 0.7, 1.9, -0.045, 7.5 and 15, respectively, under the given pressure of rupture at 673 bars. The parameter $`\tau `$ satisfies Eq. (4). The upper pressure $`p_0665`$ bars, measured before the rupture, corresponds to the value $`y_00.0119`$, which is sufficiently small, in agreement with condition (4).

## 4 Conclusion

The parameters can change for different materials as well as for different experimental setups. However, the three key physical parameters: $`\alpha `$, $`\omega `$, and $`\tau `$ should be universal, not essentially depending on the details of the external stress in the critical region. To prove this assumption, it is necessary to accomplish a set of experiments for appropriate materials. Selecting particular materials, we should keep in mind that: (i) the width of the critical region is linearly proportional to the strength of disorder, (ii) long-range elastic interactions are responsible for log-periodic corrections, (iii) the more brittle is the fracturing process, the stronger are the elastic interactions and stronger the energy release fluctuates, (iv) the fluctuations in the energy release are also stronger when a high value of a material characteristic is alternated with a low value of that characteristic. Thus, elastic-porous materials are expected to be the best choice.

Figure Captions

Fig. 1. General form of the dimensionless energy release corresponding to Eq. (6), with $`\alpha =0.7`$, $`\omega =6.9`$ and $`\tau =6.6`$ (continuous line) and to Eq. (7), with $`\alpha =0.7`$, $`\omega =8`$, and $`\tau _0=0.8`$ (broken line).

Fig. 2. Dimensionless energy release rate from Eq. (5), fitted to the experimental data reported by Johansen and Sornette (2000).


text:
# Using Wavelets to reject background in Dark Matter experiments

## 1 Introduction

Experiments searching for rare event phenomena, like those looking for non-baryonic dark matter particles (axions or WIMPs), supposedly filling a substantial part of the galactic haloes, or those looking for the neutrinoless double beta decay have become of most relevance in Particle Physics and Cosmology. To find the nature of Dark Matter, the essential component of the Universe, which in its non-baryonic form might be (according to recent developments) responsible of $``$25-30 % of the dark energy/matter budget of a flat universe, is one of the big challenges of modern Cosmology. On the other hand, the old subject of nuclear double beta decay which has been the objective of many experimental efforts as an unique tool to explore the nature of the electron neutrino and to verify the lepton number conservation, has become recently of even greater importance after the confirmation that the neutrino has a non-zero mass. In fact, the existence of neutrinoless double beta decay may provide information on the absolute scale mass of the neutrino, its mass pattern and possibly on CP violation in the leptonic sector if the sensitivity of the experiments is able to reach the level of a few meV for the Majorana neutrino effective mass upper bound.

The signals to be expected in the two selected examples of rare event physics, have in common their low rates and so the main challenge to achieve the required sensitivity is to disentangle the signal from the various sources of background, which largely hide the expected signals. Consequently, the first basic requirement of these types of searches is to provide an ultra-low background environment, as well as to design specific techniques of shielding and background rejection.

In the particular case of dark matter searches, the expected signal is concentrated in the lowest part of the spectrum  often just above the threshold of the experiment. When the effective threshold of these experiments is fixed by the electronic baseline, random fluctuations of this baseline –electronic noise– can populate the energy range of interest with noise events. Actually this electronic noise becomes a primary source of background which is of special concern.

In this paper we describe a method which has been developed and successfully applied to the last data set of the IGEX dark matter experiment. In the section 2 we will describe briefly the IGEX experiment, just to understand the needs that led to the method presented in this paper. The details of the experiment can be found in ref.  which contain also the scientific results of IGEX. In sections 3 and 4 the wavelet method itself and its calibration are described. We finish with the conclusions in section 5.

## 2 The Experiment

The IGEX experiment , optimized for detecting 76Ge double-beta decay, has been described in detail elsewhere. One of the IGEX detectors of 2.2 kg, enriched up to 86 % in 76Ge, is being used to look for WIMPs interacting coherently with the germanium nuclei. Its active mass is $`2.0`$ kg, measured with a collimated source of 152Eu. The full-width at half-maximum (FWHM) energy resolution is 2.37 keV at the 1333 keV line of 60Co, and the low energy long-term energy resolution (FWHM) is less than 1 keV at the 46.5 line of the 210Pb.

We refer to papers  and  where the latest results of the experiment (regarding dark matter searches) are presented and where the aspects related to the set-up of the experiment, shieldings, the data acquisition system, etc. are described. The threshold of the experiment is about 4 keV and the raw background registered in the region just above threshold is of a few tenths of counts per kg and day and per keV. The background and threshold achievements provide the best WIMP-nucleon cross-section versus WIMP mass ($`\sigma `$,m) exclusion plot for spin-independent interactions ever obtained by an experiment with no mechanism of nuclear/electron recoils discrimination for background rejection. As stressed in , the background level at low energies is crucial to set the sensitivity of the experiment and this is illustrated by the fact that a moderate reduction of the low energy background would lead to a substantial improvement in the final exclusion. To deal with the worrysome noise component in this energy region trying its minimization, we developed the analysis method described in the present paper. The acquisition system for low energy data is implemented by splitting the preamplifier output pulses and routing them through two Canberra 2020 amplifiers having different shaping times. These amplifier outputs are converted using 200 MHz Wilkinson-type Canberra analog-to-digital converters (ADC), controled by a PC through parallel interfaces. The energy spectrum is built using the ADC output of one of the lines. In addition, the output pulse of one of the amplifiers is digitized and recorded for each event using a 800 MHz LeCroy 9362 digital scope. The digitized pulse encompasses a time window of 500 $`\mu `$s with a sampling resolution of 1 $`\mu `$s. This recorded pulse shape is analyzed off-line with the method described below and is formally referred as $`f(t)`$ in the following.

## 3 The method

Wavelets are mathematical tools that have been proven to be extremely useful in a number of different types of analysis . In particular, they have been applied to the recognition and extraction of noise from certain sets of data. From a mathematical point of view, the so-called wavelet transforms can be viewed as extensions of the Fourier transforms, but with the addition of sensitivity to local features of the studied function. In fact, the set of functions which play the role of the $`e^{i\omega t}`$ in the Fourier transform are now $`\psi _{a,b}(\frac{tb}{a})`$ derived from an original localized null-area wavelet function $`\psi (t)`$ by means of translations and contractions parameterized respectively by $`b`$ and $`a`$. Therefore, for a fixed wavelet function $`\psi (t)`$, the biparametric wavelet transform $`[W_\psi f](a,b)`$ of a given function $`f(t)`$ can be defined as:

$$[W_\psi f](a,b)=\frac{1}{\sqrt{a}}_{\mathrm{}}^{\mathrm{}}f(t)\psi \left(\frac{tb}{a}\right)𝑑t$$
(1)
The above mentioned condition of null-area for $`\psi (t)`$, or *admisibility condition*, is derived from the requirement that no information loss must occur in the transformation (i.e. it must be invertible). Apart from this requirement, other properties are desirable in certain contexts, like continuity, gaussian shape, etc. Without entering into the mathematical details of the wavelet theory , and following a practical approach, the ”mexican hat” wavelet function $`\psi (x)`$ was chosen because it fitted most suitably our requirements:

$$\psi (x)(1x^2)\mathrm{exp}\left(\frac{1}{2}x^2\right)$$
(2)
As stated before, $`f(t)`$ is the pulse shape for each event as recorded by the data acquisition system. Following expression (1), the wavelet transform function $`[W_\psi f](a,b)`$ can be calculated numerically for each event. The resulting function is continuous due to the properties of the wavelet function used. So the relative maxima $`(a_i,b_i)_i`$ (where $`i`$ goes from 1 to $`N`$, being N the total number of maxima) can be calculated by means of an appropriate numerical algorithm. Let’s order the maxima $`(a_i,b_i)`$ according to the magnitude of the wavelet transform in that maxima, $`\omega _i[W_\psi f](a_i,b_i)`$, such that the first maximum $`\omega _1`$ is the highest one. In some sense, these maxima contain the information of the pulse. For instance, the position of the highest one in the $`(a,b)`$ plane is determined by the width and position in the time window of the event pulse itself, if present, and the others $`(a_i,b_i)_{i=2\mathrm{}N}`$ are the result of random fluctuations of the baseline along the time window encompassed by $`f(t)`$. That is illustrated in figure 1, where in the lower part a sample event pulse is shown (to give an idea of the vertical scale, the pulse correspond to an event of about 4.5 keV, i.e. just above our threshold). In the upper part of the figure the 2-dimensional wavelet transform of the same pulse is presented, lighter regions corresponding to higher values of the $`[W_\psi f](a,b)`$. The relative maxima are also marked with darker spots, the first one being clearly identified and corresponding to the true pulse, and the others corresponding to random fluctuations of the baseline.

The key idea of the method is based on the comparison between the magnitude of the first maximum $`\omega _1`$ with the distribution of the others $`\{\omega _i\}_{i=2\mathrm{}N}`$. In fact, the mentioned distribution has been checked to be well modelled by a decaying exponential (this fact being related to the kind of wavelet function we are using). Figure 2 shows such a distribution for one event.

Now, the method consists of constructing the exponential distribution $`D(w)=ke^{\alpha w}`$ which follows the baseline induced maxima in every event pulse shape and evaluate the parameter $`P_1`$ of the first maximum, defined as the normalized area of the distribution $`D(w)`$ which lies above the first maximum $`w_1`$:

$$P_1=\frac{_{w_1}^{\mathrm{}}D(w)𝑑w}{_0^{\mathrm{}}D(w)𝑑w}=e^{\alpha w_1}$$
(3)
This parameter $`P_1`$ is the final outcome of the method for each event, and gives us, to a certain extent, an indication of the probability of the main pulse to have been generated randomly by the fluctuations of the baseline, being smaller for pulses farther from a noise-like pulse. Therefore, if $`P_1`$ is large enough, the event must be considered a noise fluctuation and should be rejected, otherwise, the event is kept. The reference value to decide if $`P_1`$ is large enough or not must be fixed by calibration of the method, and is explained in the following section.

## 4 Calibration

In order to fix the rejection criterion this method was applied to a calibration set of data. The calibration was performed with a 22Na source which was introduced inside the shielding through a teflon probe producing a rate of events much higher than the standard background. The source provides an almost pure sample of ”true” events, i.e. events not coming from noise. The result of the method applied to this set of data is illustrated in figure 3 where the parameter $`P_1`$ (obtained as previously explained) versus energy is presented for each event of the calibration set. As expected, this parameter shows a clear correlation with the event energy. It is interesting to see that the same plot for the parameter $`P_2`$, defined similarly as $`P_1`$ but for the second maximum $`w_2`$, shown in figure 4 lacks completely this correlation, as this second maximum is related to some fluctuation of the baseline other than the true pulse. The same type of plot obtained for a background set of data is shown in figure 5. In this plot two different populations are clearly visible. One of them follows the same energy dependence than that of the calibration data and the other does not. This second type of events are identified as noisy fluctuations of the baseline that have triggered the acquisition system but no relevant pulse above threshold is present. From these plots a criterion of $`P_1<0.01`$ can be defined to distinguish the two populations of noise and data. Due to the limited statistics of the set of data available, it is hard to quantify the energy dependence of the loss of efficiency of the cut. However, from the presented plots one can state that it is negligible for events which lay a little above our threshold, whereas only in the very first keV above threshold it could be only a few percent.

Moreover, this method provides another way of rejection. The calculation of the first maximum of the wavelet transform provides the width of the corresponding pulse: $`a_1`$. This parameter is tightly constrained for true pulses due to the shaping. As can be seen in the figure 6, for our present set-up all of the low energy pulse widths are between 7 and 16 $`\mu `$s. This fact can be used to reject events with pulse widths outside this range. This argument allows to reject a small population of events which were not rejected by the previous method, fact that is illustrated in figure 7. They correspond to large lower-frequency fluctuations of the baseline presumably connected with microphonism. Due to the fact that they are not strictly random fluctuations in the sense previously considered, they could not be rejected before. The loss of efficiency of this cut is completely negligible.

## 5 Conclusions

We have presented a method based on wavelet techniques which is able to efficiently reject electronic noise and microphonism from the low energy –just above threshold– region. The method analyzes the recorded pulse shape of the event to extract information about the pulse height in comparison to the fluctuations of the baseline to set if the pulse can be assumed to be randomly generated by the fluctuations. This method was succesfully applied in the last data of the IGEX dark matter direct search. In figure 8 we show the low energy spectrum obtained with the last published 80 kg$``$day of data of IGEX  showing the effective rejection achieved with the present method. About half of the background between 4 and 7 keV and about 20% between 7 and 10 keV was identified as noise and eventually rejected.

## 6 Acknowledgments

This work has been partially supported by MCYT, Spain, under contract No. FPA2001-2437. We are indebted to our IGEX colleagues for their collaboration in other tasks of the experiment. I.G.I. wish to thank Prof. J. Garay for his illuminating discussions on the wavelet transform.


text:
# The Density Matrix Renormalization Group in Nuclear Physics: A Status Report

## I Introduction

The Density Matrix Renormalization Group or DMRG was originally introduced by Steven White in the early 1990s to treat the properties of quantum lattices. It quickly proved to be enormously successful, producing for the ground state energy of the spin-1 Heisenberg chain results that were accurate to 12 significant figures, well beyond what was achievable with any other approximate many-body method. Subsequently, the method was applied with great success to other 1D lattices, including spin chains and t-J and Hubbard models. The model in its original formalism, which worked in terms of real space lattice sites, has also been applied, though with much less success, to some 2D lattices as well.

Subsequently, the DMRG method was reformulated so as not to work solely in terms of real space lattice sites. In this extended version, the method has proven extremely useful in the description of finite Fermi systems, as arise for example in quantum chemistry and in the physics of small metallic grains. The successes achieved in these domains suggests the possible usefulness of the DMRG method in the description of another finite Fermi system, the atomic nucleus. In this paper, we briefly review the current status of our recent efforts to implement such a program.

The outline of the presentation is as follows. In Section II, we briefly review the key elements of the DMRG method, whether for quantum lattices or finite fermi systems. In Section III, we describe our first efforts to apply a variant of the method to real nuclei. As we will see, the results are not particularly good. In Section IV, we describe a possible way to improve the DMRG methodology, to more appropriately tune it to the the physics of nuclei. Finally, in Section V, we summarize the key points of the presentation.

## II Brief Review of the DMRG formalism

The DMRG is a method for systematically taking into account in an approximate fashion all the degrees of freedom of a problem, without letting the problem get numerically out of hand. The method is rooted in Ken Wilson’s “onion” picture, schematically illustrated in figure 1. Each layer of the onion should be thought of as another degree of freedom. It could be a site on a lattice or it could be a single-particle orbital in a problem involving a finite fermi system.

Assume that we have already taken into account a given number of layers of the onion (illustrated by solid lines in the figure) and that the total number of states we have kept to describe that portion of the system, which we call the block, is $`p`$. Assume furthermore that the next layer of the onion (illustrated by a dashed line) contains another $`s`$ states. If we add the new layer to the block, the total number of states in the enlarged block will be $`p\times s`$. The basic idea of all Renormalization Group Methods, whether Wilson’s numerical Renormalization Group (RG) or White’s DMRG, is to somehow truncate from these $`p\times s`$ states to the most important $`p`$ states, the same number we had before enlarging the block. We would then add the next layer of the onion, truncating afterwards to the most important $`p`$ states again. And we would do this until all layers have been treated. We would carry out the calculation as a function of $`p`$, stopping when the change with increasing $`p`$ is sufficiently small.

To this point, we have not yet said how to implement the truncation. Indeed, the answer is different in the usual Wilson RG method and in White’s DMRG method.

In Wilson’s method, one simply diagonalizes the hamiltonian in the enlarged block, including the extra layer, and then truncates to the lowest $`p`$ eigenstates.

In White’s DMRG method the idea is different, and better. Here, one considers the enlarged block in a medium that approximates the rest of the system. The entire system - enlarged block + medium - is referred to as the superblock. The hamiltonian is diagonalized in the superblock, yielding a ground state wave function

$$|\mathrm{\Psi }>=\underset{i=1,p\times s}{}\underset{j=1,t}{}\mathrm{\Psi }_{ij}|i>_B|j>_M$$
(1)
where $`B`$ denotes the states in the block, $`M`$ the states in the medium and $`t`$ the number of states in the medium.

The reduced density matrix for the enlarged block in the ground state,

$$\rho _{ii^{}}^B=\underset{j=1,t}{}\mathrm{\Psi }_{ij}\mathrm{\Psi }_{i^{}j}^{},$$
(2)
is then constructed and diagonalized,

$$\rho ^B|u^\alpha >_B=\omega _\alpha ^B|u^\alpha >_B.$$
(3)
Those eigenstates $`|u^\alpha >_B`$ with the largest eigenvalues $`\omega _\alpha ^B`$ are the most important states of the enlarged block in the ground state of the superblock, i.e. of the system. One thus truncates to the $`p`$ states of the enlarged block with the largest density matrix eigenvalues.

Details on how this can be efficiently done can be found in the Proceedings of this series of symposia in 2001 .

### A The finite versus infinite algorithm

The method described above, where one passes through the set of onion layers a single time, is referred to as the infinite DMRG algorithm. Depending on the manner in which correlations between layers fall off, this method can sometimes lead to an accurate representation of the ground state of the system and perhaps some excited states as well . Usually, it does not, however, since the early layers know nothing of the physics of those treated subsequently. Thus the blocks are only being optimized with respect to the layers that make them up. This suggests the use of a “sweeping” algorithm, whereby once all layers have been sampled we reverse direction and update the blocks based on the information of the previous “sweep”. Such a sweeping algorithm can be iteratively implemented until acceptable convergence in the results has been achieved. The latter procedure, in which we sweep back and forth through the onion, is called the finite algorithm. This is the method usually needed when dealing with finite fermi systems, like nuclei.

### B The p-h DMRG method

In the description of finite fermi systems, it is natural to use a basis of single-particle states and to implement the DMRG by iteratively adding their effects. An important feature of such systems is the presence of a fermi surface, which divides the set of single-particle levels into those that are primarily occupied and those that are primarily filled. It is possible to incorporate this into the DMRG algorithm in the following way.

Consider the set of single-particle levels shown in figure 2. Note that they are divided into two sets. There are particle levels - those that are primarily empty - and hole levels - those that are primarily filled. In the p-h DMRG method, one starts with the levels nearest to the fermi surface, which make up particle and hole blocks, and then gradually adds to them levels further away. The figure is meant to signify that we have already created particle and hole blocks involving the first two available levels, respectively, and that we then enlarge them by adding the third level(s). The particle block then serves as the medium for the holes and the hole block as the medium for the particles as we implement the DMRG truncation strategy. This process is continued until all particle and all hole levels have been treated.

The p-h DMRG method is particularly useful when the infinite algorithm can be used. When it is necessary to introduce sweeping, it becomes too cumbersome to implement because of the preponderance of blocks that have to be considered. This is especially true for nuclei in which one has two kinds of particles and thus twice as many blocks. In such cases, there are four distinct blocks – proton particle, proton hole, neutron particle and neutron hole – and this makes it very hard to implement sweeping.

Nevertheless, the first calculations carried out for nuclei using the DMRG method were based on the p-h algorithm. The hope was that the correlations would fall off sufficiently rapidly as we progress away from the fermi surface, making the infinite algorithm acceptable. It is these results, the first that we have obtained using the DMRG for realistic nuclei, that will be presented in the next section.

## III Application of the p-h DMRG method to $`{}_{}{}^{24}Mg`$

As a first application of the p-h DMRG method in nuclear structure, we considered the nucleus $`{}_{}{}^{24}Mg`$, with four neutrons and four protons outside doubly-magic $`{}_{}{}^{16}O`$. As is usual in the shell model, we assumed that $`{}_{}{}^{16}O`$ is inert and distributed the remaining 8 nucleons over the orbits of the $`2s1d`$ shell only. This shell-model problem is small – perhaps too small – and can be solved trivially by exact shell-model diagonalization.

The results we show were obtained using a Hartree Fock single-particle basis to define the single-particle levels and their order. These levels are illustrated in figure 3.

There is a great calculational simplification when one uses a single-particle basis in which all levels are essentially the same. And indeed when one uses an axially deformed HF basis, all levels are doubly degenerate and thus very similar in structure.

The calculations were done using the usual USD effective hamiltonian for the $`2s1d`$ shell.

In figure 4, we show results for the energies of the lowest four states of the nucleus. The results are presented as a function of the quantity $`p`$ introduced earlier, namely the number of states kept in a block. For this problem the largest $`p`$ that can be achieved is $`p=64`$. The solid line in each panel gives the exact result.

Two points should be noted. First, the results indeed converge to the exact results for $`p=64`$, as they should. What is more important to note, however, is that the convergence is very slow. We need to have values of $`p`$ above $`40`$ to get reasonably accurate results when compared with the exact values. For the values of $`p`$ needed to get an accurate reflection of the full results we must treat hamiltonian matrices not much smaller than those of the full problem. Clearly, the method – as just described – is not working very well.

### A What’s wrong?

From our perspective, there are at least two things wrong with the calculational framework that was used. On the one hand, since we worked in the m-scheme, we did not preserve angular momentum when we imposed truncation. If we want to work in the m-scheme and not lose angular momentum conservation we would have to (1) work in terms of spherical single-particle states, (2) include all states from a given orbit in a single step, and (3) make sure never to cut the truncation within a set of degenerate density matrix eigenvalues. This is possible, but extremely difficult, especially for systems with both neutrons and protons.

The second limitation of the method we employed is that it used the infinite algorithm and thus did not implement sweeping. As noted earlier, it is not practical to include sweeping in the p-h algorithm because of the preponderance of blocks that need to be coupled together. We have indeed implemented sweeping for $`{}_{}{}^{24}Mg`$, and it did improve the results slightly, but it is clear that this is not the way to proceed for more complex nuclei.

## IV What next?

Based on these considerations, we now believe that the minimum requirements for a useful DMRG strategy in nuclear physics are (1) that it works in a J-scheme or angular momentum conserving basis and (2) that it implements sweeping. In the following subsection, we sketch how this can be done. Further discussion of symmetry-conserving DMRG methods can be found in ref. .

### A The J-DMRG

#### 1 Initialization

We will illustrate the J-DMRG method through a problem involving five neutron orbits and five proton orbits, as illustrated in figure 5. These are spherical shell-model orbits, with definite $`n`$, $`l`$ and $`j`$. However, we only show $`j`$ for simplicity.

$`j_5\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}`$    $`\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}j_{10}`$  $`j_4\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}`$    $`\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}j_9`$  $`j_3\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}`$    $`\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}j_8`$  $`j_2\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}`$    $`\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}j_7`$  $`j_1\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}`$    $`\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}j_6`$  $`\nu `$    $`\pi `$

Input to the calculation includes the number of neutrons and protons, the shell-model hamiltonian, including one- and two-body terms, and a set of single-shell matrix elements for each active orbit.

The single-shell matrix elements are the reduced matrix elements of all sub-operators of the hamiltonian within the orbit. All can be readily obtained from the Coefficients of Fractional Parentage for the orbit. Such single-shell reduced matrix elements are a common feature of all J-scheme shell model codes.

#### 2 The warm-up phase

The next step is the warm-up phase, in which we carry out a first pass through the orbits and store an initial set of reduced matrix elements associated with the various possible blocks (groups of orbits). This includes for neutrons the orbits $`j_1j_2`$, $`j_1j_3`$, $`j_1j_4`$ and $`j_1j_5`$, with corresponding blocks for protons.

$`j_5\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}`$    $`\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}j_{10}`$  $`j_4\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}`$    $`\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}j_9`$  $`j_3\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}`$    $`\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}j_8`$  $`j_2\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}`$    $`\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}j_7`$  $`j_1\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}`$    $`\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}j_6`$  $`\nu `$    $`\pi `$

Consider for example the schematic illustration shown in figure 6. At this point, we have treated the first two orbits for both neutrons and protons. In each of the two blocks, we have a given number of states with each $`n`$ (the total number of particles) and $`J`$ (the total angular momentum). The number of states is analogous to the quantity $`p`$ introduced earlier. Furthermore, in the two blocks the reduced matrix elements of all sub-operators of the hamiltonian have already been stored.

We then add the next neutron level, as illustrated in figure 7.

$`j_5\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}`$    $`\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}j_{10}`$  $`j_4\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}`$    $`\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}j_9`$  $`j_3\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}`$    $`\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}j_8`$  $`j_2\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}`$    $`\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}j_7`$  $`j_1\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}`$    $`\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}j_6`$  $`\nu `$    $`\pi `$

To do so, we first construct states of good angular momentum for the enlarged neutron block. We then calculate the matrix elements of all neutron operators in the enlarged neutron block, using standard formulae given for example in ref. . The calculation requires reduced matrix elements from the block made up of orbits $`j_1`$ and $`j_2`$, which were stored in the previous iteration, and those of the orbit $`j_3`$ , which were stored in the input stage.

We then couple neutrons and protons together to states of good angular momentum, calculate the hamiltonian matrix in this basis and diagonalize. We then construct the ground state density matrix for neutrons and use it to truncate to the same number of states as we had before level $`j_3`$ was added. Then we transform all reduced matrix elements to the truncated block and store them.

And then we add the next proton level, and then the next neutron level, etc, until all levels have been included. At that point we have a first guess as to the optimum states associated with the successively larger blocks of levels.

#### 3 The sweeping phase

We now turn to the sweeping phase, schematically illustrated in figure 8. It is during this phase that we systematically improve our description of the physics of each of the blocks, by taking into account the influence of the other orbits in the problem.

The idea of the picture is as follows. We have just treated the block of proton orbits $`j_9`$ and $`j_{10}`$. We now wish to add to it the proton orbit $`j_8`$ to create an enlarged proton block. We will carry out a truncation in this enlarged block, following the density matrix strategy. To do so, we consider the block $`j_6j_7`$ as the proton medium and the block $`j_1j_5`$ as the neutron medium. We then couple the states of the enlarged proton block to the two parts of the medium to obtain the superblock, an approximation to the entire system. We calculate the hamiltonian in the associated superblock, using only stored information We then diagonalize the superblock hamiltonian and determine the density matrix for the enlarged proton block, orbits $`j_8j_{10}`$. We then use this to truncate the enlarged proton block to the same number of states as we had before enlarging it. Then, we calculate all reduced matrix elements in the truncated proton block and store them. And then we add the next orbit, $`j_7`$. We do this for all proton blocks and then for all neutron blocks. And after the sweep is finished, we simply turn around and sweep upwards, continuing the process until the results from one sweep and those from the previous sweep are acceptably close.

$`j_5\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}`$    $`\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}j_{10}`$  $`j_4\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}`$    $`\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}j_9`$  $`j_3\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}`$    $`\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}j_8`$  $`j_2\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}`$    $`\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}j_7`$  $`j_1\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}`$    $`\mathrm{\_}\mathrm{\_}\mathrm{\_}\mathrm{\_}j_6`$  $`\nu `$    $`\pi `$

The formalism, as just described, is in the process of being implemented by two of the authors (SP and JD) in collaboration with Nicu Sandulescu and Larisa Pacearescu.

## V Closing Remarks

In this presentation, we have provided a status report on recent efforts to develop the Density Matrix Renormalization Group method for use in large-scale nuclear shell-model calculations. Following a brief review of the general ideas behind the DMRG, we described the first application of the p-h variant of the method to realistic nuclei. Regrettably, the results were not especially promising. We discussed some of the reasons for the failure and then discussed a possible strategy that might overcome those shortcomings. The basic idea is to implement the DMRG algorithm in an angular-momentum-conserving basis and to include sweeping. While we are guardedly optimistic that this will indeed provide a practical and efficient methodology for large-scale nuclear structure calculations of heavy nuclei, only time will tell.

Acknowledgments This work was supported in part by the US National Science Foundation under grant #s PHY-9970749 and PHY-0140036, by the Spanish DGI under grant BFM2000-1320-C02-02, by NATO under grant PST.CLG.977000, and by the Bulgarian National Foundation for Scientific Research under Contract # $`\mathrm{\Phi }`$-809. One of the authors (SP) wishes to acknowledge fruitful discussions on the J-DMRG method with Nicu Sandulescu.