text:
# CMBEASY:: an Object Oriented Code for the Cosmic Microwave Background

## I Introduction

The cmbfast computer program for calculating cosmic microwave background (CMB) temperature and polarization anisotropy spectra, implementing the fast line-of-sight integration method, has been publicly available since 1996 Seljak:1996is . It has been widely used to calculate spectra for open, closed, and flat cosmological models containing baryons and photons, cold dark matter, massless and massive neutrinos, and a cosmological constant. It is a very well-tested program that has enabled many cosmologists to compare models of the Universe against the experimental CMB data.

However, from the point of view of code design, there are few programs that could not be improved. This is also true for cmbfast: it is a rather monolithic code that is quite difficult to oversee and modify. For example, variations in the cosmic expansion history due to dark energy necessitate similar changes to the code in numerous separate locations.

In order to address these shortcomings and simplify modifications of the code – in our case the implementation of quintessence models and gauge invariant variables – we have ported the cmbfast package to the C++ programming language. The C++ language is *object* oriented and it turns out that to think in objects (more of this soon), is very advantageous in cosmology.

The program has not been rewritten from scratch, but redesigned step by step. Some people may argue that it is hence not *independent*, i.e. some unknown errors and limitations in cmbfast could be present in the new code. The object oriented modular design, however, ensures that each part of the code is independently testable. If, for instance, one does not trust the integrator, one can use another one to check it, without changing anything else in the package. In addition most of the lines in the code *have* been rewritten, to benefit from the redesign.

cmbeasy calculates the temperature and polarization anisotropies for spatially flat universes, containing dark energy, dark matter, baryons and neutrinos. The dark energy can be a cosmological constant or quintessence. There are three scalar field field models implemented, or alternatively one may specify an equation-of-state history. Furthermore, the modular code is easily adapted to new problems in cosmology.

In the next section II.1 we explain object oriented programming and its use in cosmology. The design of cmbeasy is presented in section II.2, while the graphical user interface is discussed in section II.3. Computing precision is discussed in section II.4, while the documentation for cmbeasy is introduced in II.5. For an introduction to perturbation theory and CMB physics using the conventions of cmbeasy, as well as a summary of the evolution equations, please see the documentation on www.cmbeasy.org.

## II The code

We briefly introduce the concept of object oriented programming before turning to the design of cmbeasy. We do this using examples from the code.

### II.1 Objects

In C++, a class is a user-defined data type stroustrup . As with any data type, there can be many variables with the same type, e.g. floats, arrays, complex numbers. In the case of user-defined types, the specific variables are called objects.

#### II.1.1 Encapsulation

Quite often, some data and functions acting on the data are so tightly connected, that it is sensible to think of them as one object. As an example, let us discuss splines. Given a discrete set of $`n`$ points $`x_i`$ with $`x_i<x_{i+1}`$ and corresponding $`f(x_i)=y_i`$, a spline can smoothly interpolate, i.e. give $`f(x)`$ for any $`x[x_0,x_n]`$. Provided that the sampling is dense enough, arbitrary functions may be described by a spline for all practical purposes. This is widely used in cmbfast. For instance, the $`C_l`$’s are calculated only every $`50`$ l-values for $`l>200`$. As the spectrum is very smooth, this still gives a precise result.

Now, the function describing the photon visibility (see also Figure 1), calculated in the part of the code that tracks the thermal history of the Universe, can be used to define a spline. Without object orientation, one would need to keep track of various variables, most notably arrays for the $`x,y`$ data and derivatives needed for spline interpolation. Also, in order to assure quick access within the spline data table, one either needs to know the precise layout of the data arrays (cmbfast does this), or even more variables (storing for instance the last interpolation $`x`$ value) would be necessary. In total, this sums up to a lot of bookkeeping for a conceptually simple entity like a spline.

Alternatively, one may define a *class* holding all the necessary variables a spline needs, together with *definitions* of an interface with which other parts of the program can access and manipulate the spline data. An object behaves as described by the corresponding class. There can be an arbitrary number of objects of a certain class (just like there is one floating point type float, but many variables of *type* float in a program).111 We usually denote here (and in the code) classes with capital first letter. In some cases where there is only one object of a class used in the code, we denote the object with the same name as the class, but with lower case initial letter. Hence, the line
Cosmos cosmos;
creates an object ‘cosmos’ of the class ‘Cosmos’. The class (in our case) called Spline, can hence be viewed as yet another data type, with no more bookkeeping needed than say for a floating point number. To illustrate this, let us discuss the visibility function $`g\dot{\kappa }\mathrm{exp}(\kappa (\tau )\kappa (\tau _0))`$, where $`\kappa (\tau )`$ is the optical depth. Its typical shape is depicted in Figure 1, and its peak defines the epoch of last scattering. As soon as the Spline called visibility has been given the data, its maximum can be determined by a single line of code:
tau\_ls = visibility.maximum();
Here, the function maximum() acts on the visibility. The important point to notice is that *all* functions defined in the Spline class are immediately available to everyone who uses Splines. So, whenever one needs to find the maximum, integrate a spline, calculate the convolution of two Splines etc, this can be done using very few lines of code: the functionality is fully encapsulated in the implementation of the Spline class. Any increase in performance or sophistication of the Spline class immediately translates over to all Splines used in the program.

#### II.1.2 Inheritance

Tightly connected to the fact that data and methods are combined within one object, is the concept of *inheritance*, which proves very powerful in cosmology. A class can inherit from another class (in this context called base class). All variables and the full functionality that the base class implements are instantly available to the inheriting class,222This is as if a child was born with the whole knowledge of its parents. No training and learning would be necessary. It could instantly go and increase its capabilities starting from the level of its parents. called sub-class. The sub-class can then re-implement functions of the base class to provide a different functionality, or add new functions and variables. The important point to note is that *all* classes deriving from the same base class necessarily need to provide *all* the functions the base class provides. Hence, whenever the code uses the base class, one can substitute an inheriting class for the base class by simply changing one line of code.

As an illustration, let us look at the Perturbation class of cmbeasy which is designed to evolve the perturbation equations for one $`k`$-mode through conformal time and calculate the anisotropy sources. The Perturbation class defines functions to do these calculations which other parts of the program can be sure to find implemented in all sub-classes. In practice, there are four classes that inherit from it, for perturbations in gauge-invariant variables and in synchronous gauge both with and without quintessence (see also Figure 2). From the point of view of the rest of the program, all of them are equally well suited.333Except for the fact that if one wants quintessence, the perturbation class should of course support it. Hence, when given a Perturbation object, one may e.g. ask it to initialize the perturbation variables:
perturbation-\>initialScalarPerturbations();
Yet, depending on the specific sub-class the object belongs to 444Representing the gauge choice and quintessence support, the implementation of initialScalarPerturbations() and hence the code executed will be different. Technically this is called polymorphism.

### II.2 Design

There are roughly three main steps needed to calculate the CMB anisotropy spectrum within the line of sight strategy: (1) solve the expansion and thermal background evolution, (2) propagate the perturbation equations in Fourier space, (3) map the temperature anisotropy onto the sky today. In each step, several classes work together. To illustrate the hierarchy of the most important classes, an overview is given in Figure 2. The core part of the package is the CmbCalc class. It is the scheduler calling other objects to complete all three steps. So let us turn to step (1), the expansion history of the universe.

#### II.2.1 The background evolution

The Cosmos class is the central instance providing background quantities like the energy density $`\overline{\rho }(\tau )`$ of all species etc. This centralization of the background evolution within the Cosmos class facilitates the modification of the code significantly. A different background cosmology (such as quintessence) can be implemented by just inheriting from Cosmos and re-implementing the expansion history part of the code.555All in all 800 lines of a total of 2500 lines of Cosmos. For most other objects in the code, it makes no difference *exactly what* caused the universe to expand, as long as they can rely on Cosmos or any sub-class to provide access to the expansion history. Both the expansion and the thermal history are computed at the beginning of the calculation and stored in splines.

#### II.2.2 Propagation in k-space

All fluctuation equations are handled by Perturbation and its sub-classes. Its most prominent function is propagateScalar($`\tau _1,\tau _2`$)  which evolves the perturbations from conformal time $`\tau _1`$ to $`\tau _2`$. For numerical stability, one usually treats the era of tight coupling between baryons and photons differently. However, this results in a non-analytic jump in the r.h.s of the perturbation equations at $`\tau _{tc}`$. Some integrators for ordinary differential equations get confused by this. Therefore, propagateScalar() will determine if the jump occurs within \[$`\tau _1,\tau _2`$\] and if so, integrate up to $`\tau _{tc}`$, make a tiny666More specific, it adds $`\mathrm{\Delta }\tau `$, such that for double precision variables $`\tau _{tc}+\frac{1}{2}\mathrm{\Delta }\tau =\tau _{tc}`$. jump and integrate on till $`\tau _2`$. At the end of each propagation, the sources of cmb anisotropies are calculated by the scalarSources() function.

#### II.2.3 The Cl’s

The convolution of the sources with the Bessel functions, as well as the final $`k`$-integration are performed in the ScalarIntegrator sub-class. In contrast to cmbfast, we formulate the convolution with the $`j_l^{}s`$ as an ordinary differential equation (ODE). In fact, the Spline class provides the convolution functionality. If higher precision in this step is wanted, one may simply increase the desired precision for Spline, which will pass this through to the integrator of ODE’s. As far as the $`k`$-integration is concerned, we stick to the interpolation points of cmbfast. We have written two other implementations of ScalarIntegrator, although none could compete with cmbfast’s speed. Furthermore, no gain in accuracy justifying the speed tradeoff was found. As far as the $`k`$-integration is concerned, a fast and independent realization would be desirable. With the high level features of C++, algorithms that seem hard to implement in other languages may be feasible.

#### II.2.4 Little helpers

The most valuable class assisting the computation may be Spline. As splines are used all over the code, and quite often in time-critical environments, the main design goal has been speed. In order to perform a spline interpolation at $`x`$, one needs to find the interval $`x[x_i,x_{i+1}]`$. While cmbfast hard-wires the size and layout of the $`x_i`$-data, this has not been an option for a versatile spline implementation. Yet, in all time critical situations within our calculation, two facts come to the rescue: Firstly, interpolations occur in order and close to the last interpolation. Hence, the next value of $`x`$ is likely to lie within the same interval (or close to it). Secondly, many splines with exactly the same $`x`$-data are evaluated at the same $`x`$. For instance, the r.h.s of the perturbation equations need background quantities all at the same $`\tau `$. As all background quantities get their data at the same $`\tau _i`$ in Cosmos, they necessarily share the same $`x`$-data. Therefore, some parts of the interpolation formulae have to be computed only once777With “once”, we mean until interpolation is requested for a different $`\tau `$. for each $`\tau `$. In addition, finding the position of the interpolation interval is simple: it will be the same for the same $`\tau `$.

To accommodate for these groups of splines and caching, splines come in two flavors: with and without its own $`x`$-data. While a spline with its own $`x`$-data can live on its own, a spline without $`x`$-data can only be created if a “mother” spline providing the $`x`$-data is given upon the creation of the “child” spline. Within these “families” of shared $`x`$-data, interpolation at the same $`x`$ is particularly fast.

Not only must a spline be fast, it must also be robust. Hence, a spline will allocate more memory if needed888and free unused memory, check that no access to data occurs outside its bound, check that all data is ordered $`x_i<x_{i+1}`$ and make sure that the interpolation tables have been built before any interpolation takes place. Splines can easily be visualized. To write the data999and interpolated values, if desired to a file, simply call
myspline.dump("filename")
The anisotropy sources, are evaluated on a $`\{k,\tau \}`$ grid. For the final $`\tau `$ and $`k`$ integration, one needs to interpolate at intermediate points. This is conveniently done using the class SplineWeb. It is essentially an array of Splines, one for each $`k`$ and one for each $`\tau `$ value. If one needs to interpolate at a given $`k`$ value, one can ask SplineWeb to return a Spline with x-data in the $`\tau `$ direction:
spline\_at\_k = web.createAlongY(k);
This is used to get the splines which are convoluted with the Bessel functions. It is also very convenient for making slices through the cold dark matter power spectrum (at either fixed $`k`$ or $`\tau `$). Given for instance the cold dark matter SplineWeb power, one can dump it to a file:
power.dumpAlongX($`\tau `$,"filename");
The above command would save the spectrum at arbitrary conformal time $`\tau `$ to disk.

There are some more of these “helpers”, most of them small classes that sometimes make life much easier. The last one we will introduce here is Anchor which keeps track of objects that are dynamically created (as are many splines). The objects it is asked to keep track of are deleted by it on demand, or when the Anchor object itself is deleted. Therefore, it greatly simplifies memory handling and prevents memory leaks.

#### II.2.5 Quintessence Implementation

The different background evolution of quintessence Wetterich:1988fm ; Ratra:1988rm ; Caldwell:1998ii  scenarios is implemented using the QuintCosmos and the Quintessence class. Each subclass of Quintessence corresponds to a certain model, such as the exponential potential Wetterich:1988fm ; Peebles:1988ek ; Ratra:1988rm , inverse power law Ratra:1988rm , leaping kinetic term Hebecker:2001zb  etc. Certainly, a more monolithic design with the quintessence models implemented in the QuintCosmos class would have been possible. However, we believe that the *details* of a model are best kept to a class of its own. For instance tuning model parameters in order to get the right amount of $`\mathrm{\Omega }_0^\phi `$ is different for each model and a monolithic design would have to call *differently named* functions for *different* models. Using sub-classing, QuintCosmos (and Perturbation) always calls functions with the same argument and name for all models. Yet, as the object implementing the function differs for each model, the code executed by calling the function can be totally different. Thus, a new quintessence model can be implemented by simply sub-classing Quintessence and implementing a minimal set of functions, such as one for the scalar field potential.

### II.3 Graphical User Interface

For educational purposes and also to simplify the parameter input and subsequent visualization of results, a graphical user interface (GUI) is of great value. Luckily, there is the very sophisticated and publicly available ‘Qt’ library troll  with which the creation of a GUI is facilitated. Its object oriented C++ design is a perfect match for the cmbeasy package. There is therefore an executable program called ‘cmbeasy’ giving interactive access to almost the full capabilities of the package, including quintessence. As is seen from Figure 3, the spectra are visualized in separate plots arranged in a so called Tab-Widget.101010A widget is a part of a user interface that can interact. Examples are buttons, sliders, etc. One can for instance zoom in, select and save curves or print the plot.

### II.4 Computing precision

Comparing the results from cmbfast to those obtained from cmbeasy, one generally notes accordance to better than $`0.5\%`$ (see also Figure 4).

In principle, the line of sight algorithm is equivalent to evolving a large hierarchy of multipoles. In practice, as one truncates the hierarchy in the line of sight approach at some low $`l`$, there is a reflection of power after the higher multipoles gradually build up. However, as by far the largest contribution towards the anisotropy occurs around recombination (when the higher multipoles are still small), this does not spoil the result dramatically. Yet, to check this, one may increase the number of multipoles involved. This is true both for cmbfast and cmbeasy.

Within cmbeasy, there are five instances in which the precision may be increased. In three of the five cases, changing one variable suffices. The first case is the background evolution. Here one may want to increase (or decrease) the precision requested from the ODE integrator. During perturbation evolution in $`k`$-space, the precision requirement given to Perturbation can be changed in a similar manner. The same is true for the convolution of the sources with the Bessel functions, taking place in ScalarIntegrator. Finally, and a bit more difficult, one may change the number and position of $`k`$-values contributing towards the $`k`$-integration. This can be done by modifying Integrator, the base class of ScalarIntegrator. Alternatively, one could devise a different implementation for the $`k`$-integration111111i.e. a new sub-class of Integrator. Finally, one may increase the number of $`l`$’s at which the $`C_l`$’s are calculated.

### II.5 Documentation

Using the doxygen program, the documentation is automatically generated from the source code of the cmbeasy package. In its html version, it is interactively navigable and includes the full source code. Due to the automatic generation, the documentation and the code are naturally synchronized. Instructions for setting up and running the program are available in this html documentation. In addition, a postscript version of the documentation is generated.

## III Conclusions

We have devised an object oriented code for calculating cosmic microwave background anisotropies. The main benefits of this code derive from its modularity.

The modular design allows one to test and refine each part of cmbeasy individually. Different modules, for different numerical methods and theoretical approaches, can easily be switched. This feature should prove extremely useful in the development and refinement of the code. As cosmological theory reaches for a new level of precision to match the new generation of experiments, map MAP , planck PLANCK  and other experiments will soon put our predictions and hypotheses to the test, so future work will seek to improve the accuracy of the code. (At the present time, cmbeasy and the standard cmbfast are generally in agreement to better than $`0.5\%`$.)

The modular design allows one to more easily develop and build upon cmbeasy, to implement new cosmological models and include more physics. As our theoretical cosmological models are a work-in-progress — as CMB data comes in — so are our tools for modeling the Universe. And so it makes sense to use a program design which makes this development easy. Hence, the functional blocks used in cmbeasy considerably reduce the time to implement new cosmological models. For most purposes, a new model can be achieved by sub-classing from existing solutions, thus eliminating the need to rewrite large parts of the code. Furthermore, new modules can be easily introduced to carry out calculations that rely on the linear perturbation spectrum. Examples include lensing, secondary sources in the CMB, reionization, and cross-correlation with other tracers of large scale structure.

An additional benefit of the cmbeasy package is the graphical front-end. As we have described, the GUI allows a user to visually dial through a range of parameters. We expect that this functionality will be useful for developing understanding and intuition for CMB physics and cosmological models.

As stated at the outset, there is, perhaps, no program design that could not be improved. However, we hope that cmbeasy will serve the cosmology community well.

###### Acknowledgements.

I would like to thank Luca Amendola and Ruth Durrer for their help, Uros̃ Seljak and Matias Zaldarriaga for cmbfast and their permission to distribute cmbeasy under its new name. I would like to thank R. R. Caldwell, Martin Kunz, Christian Müller, Gregor Schäfer and C. Wetterich for useful discussions. Finally I would like to thank the free software community for their excellent compilers, libraries and desktop environments. This work was supported by NSF grant PHY-0099543 at Dartmouth and PHY-9907949 at the KITP and the Graduiertenkolleg des Instituts für Theoretische Physik der Universität Heidelberg.


format: text

text:
# Renormalization of the isovector 𝜋𝑁 amplitude in pionic atoms

## I Introduction

In the past few years there has been renewed interest  in the extraction of the isovector $`s`$-wave $`\pi N`$ scattering amplitude near threshold from pionic atom observables. Of particular interest is the relationship between this in-medium amplitude and the corresponding amplitude for the free $`\pi N`$ interaction at threshold. This renewed interest in pionic atoms in general, and in the $`s`$-wave part of the potential in particular, stems from three recent developments. The first one is the experimental observation of ‘deeply bound’ pionic atom states in the (d,3He) reaction , the existence of which was predicted a decade earlier . The second one is the very accurate measurements of the shift and width of the 1$`s`$ level in pionic hydrogen  and in pionic deuterium  which leads to precise values for the $`s`$-wave $`\pi N`$ scattering lengths (see also Ref. ). The third development is the attempt to explain the ‘anomalous’ $`s`$-wave repulsion  in terms of a density dependence of the pion decay constant , or very recently by constructing the $`\pi N`$ amplitude near threshold within a systematic chiral perturbation expansion  and in particular imposing on it gauge invariance .

The so-called anomalous repulsion of the $`s`$-wave pionic atom potential is the empirical finding, from fits of optical potentials to pionic atom data, that the strength of the repulsive $`s`$-wave potential inside nuclei is typically twice as large as is expected on the basis of the free $`\pi N`$ interaction. It is interesting to note that early indications for a very repulsive isovector $`s`$-wave amplitude go as far back as 1978  when high precision data became available. Subsequent analyses  also indicated such a repulsive amplitude. A key point in this problem are the uncertainties which are associated with the empirical determination of the relevant potential parameters. This problem has been aggravated by recent works in which the ‘deeply bound’ pionic atom states were granted special status and unfounded implicit assumptions were made on the pion-nucleus optical potential (e.g. Refs. ). The present work addresses all these issues with the aim of obtaining reliable parameter values for the $`s`$-wave part of the pion-nucleus optical potential at threshold, together with realistic estimates of uncertainties. The focus of this work is on the extraction of $`b_1`$, the in-medium isovector $`s`$-wave $`\pi N`$ amplitude, from pionic atom data, using the conventional (density-independent) $`\pi `$-nucleus optical potential. It is shown that no definitive conclusions can be reached on whether or not $`b_1`$ is renormalized with respect to its free-space value unless the whole bulk of pionic-atom data across the periodic table is used in the analysis. In particular, the deeply bound 1$`s`$ states observed recently in the (d,3He) reaction  do not offer conclusive evidence for a renormalization of $`b_1`$. That applies also to analyses using the deeply bound states together with ‘normal’ 1$`s`$ states .

Section II presents briefly the conventional pion-nucleus potential  where the various parameters are defined. Section III is devoted to a discussion of the proton and neutron density distributions which are an essential ingredient of the pion-nucleus optical potential. The proton densities obviously are obtained from the quite well known charge distributions of nuclei but the neutron distributions are generally not known to a sufficient accuracy. For this reason simple but physical parameterizations of the neutron density distributions are introduced, in comparison with available data and with nuclear models. The dependence of the isovector $`s`$-wave $`\pi N`$ amplitude on the neutron distributions was only marginally touched upon in recent works  and, therefore, we study in Section IV in some detail the dependence of this amplitude on both the root mean square (rms) radius of the distribution and, to some extent, on the shape of the distribution. This dependence is translated then into an additional uncertainty of the extracted parameter values. Section IV also presents the results of fits to the data, starting with ‘global’ fits to 100 and 120 data points and gradually reducing the extent of the data base in order to assess the accuracy of the results, in comparison with other works using ‘partial’ fits to very limited sets of data  that place special emphasis on the deeply bound 1$`s`$ states in 115,119,123Sn and in 205Pb. The last subsection of Section IV discusses the Seki-Masutani linearization scheme  demonstrating the limits of its applicability to pionic-atom studies. In Section V we explain why the deeply bound pionic atom states generally have failed to provide new information on the pion-nucleus interaction. A brief summary and conclusions are presented in Sect. VI.

## II The pion-nucleus potential

The pion-nucleus potential has been extensively discussed (e.g. Ref. ), but for completeness a short summary follows.

The interaction of pions at threshold with the nucleus is described by the Klein-Gordon equation of the form:

$$\left[^22\mu (B+V_{\mathrm{opt}}+V_c)+(V_c+B)^2\right]\psi =0(\mathrm{}=c=1)$$
(1)
where $`\mu `$ is the pion-nucleus reduced mass, $`B`$ is the complex binding energy and $`V_c`$ is the finite-size Coulomb interaction of the pion with the nucleus, including vacuum-polarization terms. The optical potential $`V_{\mathrm{opt}}`$ used in the present work is of the canonical form due to Ericson and Ericson ,

$$2\mu V_{\mathrm{opt}}(r)=q(r)+\stackrel{}{}\alpha (r)\stackrel{}{}$$
(2)
with

$`q(r)`$ $`=`$ $`4\pi (1+{\displaystyle \frac{\mu }{M}})\{\overline{b}_0(r)[\rho _n(r)+\rho _p(r)]+b_1[\rho _n(r)\rho _p(r)]\}`$ (4)
$`4\pi (1+{\displaystyle \frac{\mu }{2M}})4B_0\rho _n(r)\rho _p(r),`$
$$\alpha (r)=\frac{\alpha _1(r)}{1+\frac{1}{3}\xi \alpha _1(r)}+\alpha _2(r),$$
(5)
where

$$\alpha _1(r)=4\pi (1+\frac{\mu }{M})^1\{c_0[\rho _n(r)+\rho _p(r)]+c_1[\rho _n(r)\rho _p(r)]\},$$
(6)
$$\alpha _2(r)=4\pi (1+\frac{\mu }{2M})^14C_0\rho _n(r)\rho _p(r).$$
(7)
In these expressions $`\rho _n`$ and $`\rho _p`$ are the neutron and proton density distributions normalized to the number of neutrons $`N`$ and number of protons $`Z`$, respectively, and $`M`$ is the mass of the nucleon, $`q(r)`$ is referred to as the $`s`$-wave potential term and $`\alpha (r)`$ is referred to as the $`p`$-wave potential term. The function $`\overline{b}_0(r)`$ in Eq. (4) is given in terms of the local Fermi momentum $`k_\mathrm{F}(r)`$ corresponding to the isoscalar nucleon density distribution:

$$\overline{b}_0(r)=b_0\frac{3}{2\pi }(b_0^2+2b_1^2)k_\mathrm{F}(r),$$
(8)
where $`b_0`$ and $`b_1`$ are minus the pion-nucleon isoscalar and isovector effective scattering lengths, respectively. The quadratic terms in $`b_0`$ and $`b_1`$ represent double-scattering modifications of $`b_0`$. In particular, the $`b_1^2`$ term represents a sizable correction to the nearly vanishing linear $`b_0`$ term. Similar modifications of $`b_1`$  were found by us to be negligibly small. The coefficients $`c_0`$ and $`c_1`$ in Eq. (6) are the pion-nucleon isoscalar and isovector $`p`$-wave scattering volumes, respectively. The parameters $`B_0`$ and $`C_0`$ in Eqs. (4) and (7) represent $`s`$-wave and $`p`$-wave absorption, respectively, on pairs of nucleons and as such they have imaginary parts. Dispersive real parts are found to play an important role in pionic atom potentials. The terms with $`4\rho _n\rho _p`$ were originally written as $`(\rho _n+\rho _p)^2`$, but the results hardly depend on which form is used. The parameter $`\xi `$ in Eq. (5) is the usual Ericson-Ericson Lorentz-Lorenz (EELL) coefficient . An additional relatively small term, known as the ‘angle-transformation’ term (see Eq.(24) of ), is also included.

The pion-nucleus optical potential contains nine parameters but only six parameters were varied in the present fits, namely, all the parameters of the $`s`$-wave part of the potential and the empirical complex $`p`$-wave absorption parameter $`C_0`$. The linear $`p`$-wave parameters $`c_0`$ and $`c_1`$ were kept fixed at their respective free $`\pi N`$ values, as discussed extensively in Ref. . The EELL coefficient was held fixed at the value $`\xi `$ = 1. It was demonstrated in Ref.  that the extracted $`s`$-wave parameters are insensitive, with negligibly small variations, to sizable variations of the $`p`$-wave parameters, provided a fit to the data was made. Note that the extracted parameters $`b_0`$ and $`b_1`$ are some averages over nuclear densities of in-medium quantities, and if they turn out to be different from the corresponding free $`\pi N`$ values then one could refer to a renormalization of the free space values in the nuclear medium.

## III Nuclear density distributions

It is evident from the previous section that the nuclear density distributions $`\rho _p`$ and $`\rho _n`$ are an essential part of the pion-nucleus potential. Nevertheless, their effect in the present context of studying the $`s`$-wave part of the potential and in particular the dependence of the parameter $`b_1`$ on the choice of nuclear densities has been hardly addressed, and only briefly so in Refs. . This point is studied in more detail in the present work.

Nuclear charge densities are known quite well  from studies of elastic scattering of electrons and from finite-size effects in muonic atoms. Density distributions for protons can then be obtained by unfolding the finite size of the charge of the proton. This leads to reliable proton densities, particularly in the surface region, which is the relevant region for producing strong-interaction effects in pionic atoms. The neutron distributions are, however, generally not known to sufficient accuracy. A host of different methods have been applied to the extraction of rms radii of neutron distributions in nuclei , including some studies using pionic atoms , but the results are sometimes conflicting. Several nuclei have been studied more than others, e.g. isotopes of Ca, Sn and Pb  but for others there is no experimental information on neutron densities and one must then rely on calculations based on various models. To complicate things further we note that there is a long history of conflict between values of neutron rms radii derived from experiments using hadronic projectiles and neutron rms radii derived from theoretical calculations. For that reason we have adopted a semi-phenomenological approach that covers a broad range of possible neutron density distributions.

It is shown below that the feature of neutron density distributions which is most effective in determining strong interaction level shifts and widths in pionic atoms is the radial extent, as represented, for example, by $`r_n`$, the neutron density rms radius. Other features such as the detailed shape of the distribution have only minor effect. For that reason we chose the rms radius as the prime parameter in the present study. Since $`r_p`$, the rms radius for the proton density distribution, is considered to be known, we focus attention on values of the difference $`r_nr_p`$. In the analysis of Ref.  the nuclear densities were either of the single particle (SP) variety or of the macroscopic type  where the values of $`r_nr_p`$ were taken from the relativistic mean field (RMF) results of Lalazissis et al.  for each of the nuclei involved. Smaller values of $`r_nr_p`$ were also used in Ref. . In the present work we vary, among other things, the values of $`r_nr_p`$ and we therefore look for a simple parameterization that will be easy to relate to the RMF results. Choosing 43 stable nuclei along the periodic table, the calculated RMF values  of $`r_nr_p`$ are described very well by the expression

$$r_nr_p=\alpha \frac{NZ}{A}+\beta .$$
(9)
Figure 1 shows the results of this linear fit in the neutron-proton asymmetry $`(NZ)/A`$, where we have arbitrarily assigned an error of $`\pm `$0.03 fm to each value of $`r_nr_p`$, as a representative value of discrepancies between calculated and experimental charge rms radii. The RMF fit parameters obtained are $`\alpha `$=1.51$`\pm `$0.07 fm, $`\beta `$=$``$0.03$`\pm `$0.01 fm.

An expression of the form (9) has also been used recently in connection with antiprotonic atoms , with parameter values of $`\alpha `$=1.01$`\pm `$0.15 fm, $`\beta `$=$``$0.04$`\pm `$0.03 fm. Note that antiprotonic atoms are sensitive to nuclear densities well below 10% of the central density and as such it is not very reliable to infer from these data quantities which relate to the bulk of the density. Nevertheless we consider these results as an experimental indication on values of $`r_nr_p`$, particularly as they have been used  in analyses of pionic atoms. In the present work the neutron rms radii were then varied using Eq.(9) with $`\beta `$=$`0.035`$ fm and varying the parameter $`\alpha `$ between 0.3 and 2.7 fm. This way we have covered a very wide range of possible values of neutron rms radii.

In order to allow for possible differences in the shape of the neutron distribution, the ‘skin’ and the ‘halo’ forms of Ref.  were used, as well as an average between the two. Assuming a two-parameter Fermi distribution both for the proton (unfolded from the charge distribution) and for the neutron density distributions

$$\rho _{n,p}(r)=\frac{\rho _{0n,0p}}{1+\mathrm{exp}((rR_{n,p})/a_{n,p})},$$
(10)
then for each value of $`r_nr_p`$ in the skin form the same diffuseness parameter for the protons and the neutrons $`a_n=a_p`$ was used and the $`R_n`$ parameter was determined from the rms radius $`r_n`$. In the halo form the same radius parameter $`R_n=R_p`$ was assumed and $`a_n^\mathrm{h}`$ was determined from $`r_n`$. In the ‘average’ option the diffuseness parameter was set to be the average of the above two diffuseness parameters $`a_n^{\mathrm{ave}}=(a_p+a_n^\mathrm{h})/2`$ and the radius parameter $`R_n`$ was determined from the rms radius $`r_n`$. In this way we have used three shapes of the neutron distribution for each value of its rms radius all along the periodic table.

## IV Results and discussion

### A Dependence on neutron distribution

For the study of the dependence on the neutron density distribution using the parameterization of Eq. (9) above we first chose a data base of 100 points, from 20Ne to 238U, including the recently published results for the deeply bound 1$`s`$ states in 115,119,123Sn  and those for the 1$`s`$ and 2$`p`$ states in 205Pb . Ne was chosen as the lightest element because it is the first element in the periodic table for which the RMF results  apply and because we use the two-parameter Fermi parameterization for the densities. Fits to the data were made by varying six of the parameters of the potential, as discussed is Sect. II.

Figure 2 shows $`\chi ^2`$ values for the three types of neutron density, namely, skin, halo and average, as a function of $`\alpha `$, the coefficient of the asymmetry parameter $`(NZ)/A`$ which determines the value of $`r_nr_p`$, see Eq. (9). Also marked by asterisks are the points associated with the particular values of $`\alpha `$=1.01 fm and $`\alpha `$=1.51 fm, which refer to the ‘$`\overline{p}`$’ and the RMF prescriptions, respectively. It is easy to conclude from this figure that in the context of global fits to pionic atom data the RMF radii and the skin shape are the favoured option. Note that the minimum of $`\chi ^2`$ implies $`\chi ^2`$/F, the $`\chi ^2`$ per degree of freedom, of 1.8, which represents a very good fit. Increasing the data base to 120 points by adding 20 points from 12C to 19F and repeating the fits, the overall picture remains essentially the same, but the minima become less sharp and imply slightly larger values of $`\chi ^2`$/F. For these additional nuclei (except 19F) a modified harmonic oscillator form  was used for the densities.

The above results may be interpreted as indicating that on the average the difference $`r_nr_p`$ is described rather well by the RMF model ($`\alpha `$=1.5 fm). However, when other experimental data had been analyzed in terms of the difference between rms radii , quite often smaller values were obtained for $`r_nr_p`$, corresponding to $`\alpha `$ between 1 and 1.5 in Eq. (9). Note that Brueckner-Hartree-Fock calculations generally yield values of $`r_nr_p`$ corresponding to $`\alpha `$ 1.0 . We therefore conclude that on the average, rms radii of neutron density distributions are between the values inferred from analyses of antiprotonic atoms ($`\alpha `$ =1.0)  and the predictions of RMF calculations .

Turning to the question of the dependence of extracted parameters of the pion-nucleus optical potential on the neutron distribution used in the analysis, Fig. 3 shows the values of $`b_1`$ which correspond to the $`\chi ^2`$ results displayed in Fig. 2. It is evident that there is a fairly strong dependence of the extracted $`b_1`$ on the prescription chosen for $`r_nr_p`$, whereas the dependence on the shape of the density distribution, as indicated by the scatter of points for each value of $`\alpha `$, is small. A similar plot of the results of fits to the 120 data points mentioned above is almost indistinguishable from the present figure. Adopting an average value of $`\alpha `$= 1.25$`\pm `$0.25 fm as discussed above, one can estimate the additional uncertainty of the extracted potential parameters, due to neutron densities, with the help of plots like the one shown in Fig. 3.

### B Reduced data bases

Several recent publications  have claimed exceedingly small uncertainties for the parameter $`b_1`$ when it is extracted from data bases which are very small in comparison with the 100-120 points of data used above. Such claims are based on adopting several linearizations of terms in the optical potential Eq. (4), notably the Seki-Masutani correlation , which are at variance with the results of global fits to large scale data bases, see Refs.  and the following subsection. We have, therefore, repeated the analysis using progressively reduced data bases in order to examine in a consistent manner the dependence of the uncertainties on the data, including also the uncertainty in the neutron density ditribution. The results of this survey are summarized in Tables I and II. The results in Table I are based on the use of the ‘RMF skin’ neutron densities and the results in Table II are based on the use of the ‘$`\overline{p}`$ skin’ densities. The ‘global 2’ and ‘global 3’ data sets are the ones mentioned in the previous subsection (‘global 1’ was used in Ref. ). Three ‘partial’ data sets are included in the tables, as follows: (i)1$`s`$ states in light pionic atoms from 12C to 28Si; (ii)1$`s`$ states in light $`N=Z`$ pionic atoms plus the deeply bound 1$`s`$ states in 115,119,123Sn  and 205Pb ; and (iii)only the deeply bound states of (ii). Within each table it is seen that the parameter values are consistent with each other as the size of the data base is changed, but the uncertainties increase as the data base becomes smaller. Note that Re$`B_0`$ is different from zero, within errors, only for the global fits. Claims based on ‘partial’ data sets that no Re$`B_0`$ is required to fit the data therefore do not hold for pionic atoms in general. In Tables I and II, for the fit to only deeply bound states, we have set Re$`B_0`$=0 since otherwise the convergence proved unreliable with so few data points.

Turning to the role played by the neutron density distributions, it is evident from Tables I and II that the dependence of the extracted parameters on the rms radius of the neutron density distribution must be taken into account in the evaluation of the final results. Table III summarizes the final values of $`b_1`$, for $`\alpha `$=1.25$`\pm `$0.25 fm (Eq. (9)) as discussed above. Since one deals here with an average behaviour of $`r_nr_p`$ along the periodic table, the effects of the uncertainty $`\mathrm{\Delta }\alpha =\pm 0.25`$ fm are included in quadrature. Again all the different results for $`b_1`$ are consistent with each other, but only those due to global fits to 100-120 data points may be regarded as significantly different from the free $`\pi N`$ value  of $`b_1^{\mathrm{free}}=0.0885_{0.0021}^{+0.0010}m_\pi ^1`$. Claiming such disagreement by analysing reduced data sets can only be described as unfounded.

### C The Seki-Masutani linearization

Having demonstrated the dependence of the uncertainties of the derived values of $`b_1`$ on the size of the data base, we turn now to another problem which has affected recent derivations of this parameter , namely, the use of the Seki-Masutani (SM) linearization . In the SM scheme the Re$`B_0`$ term which is quadratic in the densities is linearized with the help of an ‘effective’ density such that the two terms with $`b_0`$ and Re$`B_0`$ are lumped together into a single term $`b_0^{}\rho (r)`$. The existence of correlations between the parameters $`b_0`$ and Re$`B_0`$ has been known for a long time  and it can be seen also in Tables I and II. However, the same tables show also that when large data bases are being used, then the value of each one of these two parameters is reasonably well determined, although not to a very good accuracy. In order to explicitly show the consequences of using the SM linearization we have repeated the fits reported in Table I with Re$`B_0`$ set to zero. The results are summarized in Table IV where it is seen that for the large data bases the values of $`\chi ^2`$ have increased significantly. Note that the ‘natural’ unit for measuring these increases is $`\chi ^2`$/F, the $`\chi ^2`$ per degree of freedom. For the smaller data sets the increases are not statistically significant, in full agreement with the uncertainties for Re$`B_0`$ quoted in Tables I and II. It is concluded that the SM correlation is valid, as may be expected for any non-linear problem over some range of parameters, but when the quality of fit $`\chi ^2`$ criteria are taken into consideration then the use of the SM linearization is unjustified and in the case of large data bases it biases the results.

Another interesting conclusion obtained from Table IV is that when Re$`B_0`$ is held fixed the values of $`b_1`$ and their uncertainties hardly change compared to the unrestricted analysis. The values of $`b_0`$ naturally change in order to compensate for the corresponding changes in Re$`B_0`$ but the uncertainties of $`b_0`$ go down by almost an order of magnitude, which is unphysical. Since the derivation of $`b_1`$ by Geissel et al.  is based on “$`b_0^{}b_1`$ constraints” ($`b_0^{}`$ closely related to our $`b_0`$), we conclude that the uncertainties of $`b_1`$ obtained in such limited analysis of deeply bound states, which is based on the SM linearization, are bound to be artificially small.

## V The role of deeply bound states

In order to understand the role played by the deeply bound states in pionic atoms of 115,119,123Sn and of 205Pb in global fits, we compare values of $`\chi ^2`$/N, the total $`\chi ^2`$ per point obtained in the fits, with the corresponding values for the contribution of the deeply bound states. For the ‘global 2’ fit, using RMF-skin densities (Table I), we have $`\chi ^2`$/N=2.0 whereas the contribution of the deeply bound states amounts to 1.0 per point for these states. For the ‘global 3’ fit the numbers are 1.7 and 0.8, respectively. Removing the deeply bound states from the ‘global 2’ data and repeating the fit, we find $`\chi ^2`$/N=2.0 and then using the resulting parameters to predict the shifts and widths of the deeply bound states, we get $`\chi ^2`$/N=1.2 for these states. Therefore in all three examples the $`\chi ^2`$/N for the deeply bound states are smaller than the average, signifying that the deeply bound states do not depart from the picture obtained for normal pionic atoms along the periodic table. Limiting the discussion to only 1$`s`$ states over the periodic table and confronting in Table III the results for light nuclei only with the results for light $`N=Z`$ nuclei plus the deeply bound states (the latter subset was advocated in Ref.), again it is concluded that the deeply bound states do not provide new information on the isovector $`s`$-wave amplitude $`b_1`$. This observation is not entirely new and it was shown already in  that the deeply bound states in 207Pb were in line with the expectations based on normal pionic atom data. A possible explanation was also indicated in Ref. , namely, “the same mechanism which causes the deeply bound states to be narrow also masks the deep interior of nuclei where new effects could possibly be observed”.

The fact that the deeply bound pionic atom states have failed to provide new information on the pion-nucleus interaction can be understood from simple arguments of overlap between the atomic wavefunction and the nucleus . Figure 4 displays absolute values squared of the atomic radial wavefunctions multiplied by the nuclear density for a normal 1$`s`$ state (in 20Ne) and for a deeply bound 1$`s`$ state (in 208Pb). It is seen that the Coulomb wavefunction would have indeed penetrated deeply into the heavy Pb nucleus, but due to the strong interaction it is repelled such that its overlap with the nucleus is sufficiently small to make the width of the state relatively narrow and thus making the state observable. In fact, the deeply bound 1$`s`$ wavefunction does not overlap with inner regions of the nucleus more so than a normal 1$`s`$ wavefunction does. This shows clearly why deeply bound states do not play any special role in the determination of pionic atom potentials. Expanding the statement made in Ref.  we note that the same mechanism which causes the deeply bound states to be narrow and observable, namely, the strong repulsion  of the wavefunction out of the nucleus, also masks the nuclear interior such that the penetration of the deeply bound pionic atom wavefunction is not dramatically enhanced compared to the normal states.

## VI Summary and conclusions

In this work we have presented a comprehensive analysis of uncertainties for the parameters of the pion-nucleus optical potential which are extracted from experimental strong interaction effects in pionic atoms. Emphasis was placed on the parameter $`b_1`$, the isovector $`s`$-wave $`\pi N`$ amplitude, which has been discussed widely in recent years in connection with possible medium modifications that could also explain the long-standing ‘anomalous’ repulsion. In addition to errors which are inherent in the $`\chi ^2`$ fit procedure, we have examined errors which are introduced because of the insufficient knowledge of neutron density distributions in nuclei, both in shape and in radial extent. The uncertainties involved in the use of partial data sets were compared with those obtained from global analyses of large amount of data covering all the periodic table. It was shown that only in the latter case the uncertainties are sufficiently small that a discrepancy between the value $`b_1=0.108\pm 0.007m_\pi ^1`$ derived from pionic atoms and the corresponding value for the free $`\pi N`$ interaction can be established. In particular, the deeply bound states do not provide sufficiently extensive and varied data base in order to constrain $`b_1`$ to the required precision. The failure of the deeply bound states to provide meaningful new information on the pion-nucleus interaction has been demonstrated and a simple explanation in terms of overlap arguments was offered.

Finally, we mention Weise’s suggestion for a density dependence of $`b_1`$ motivated by a partial restoration of chiral symmetry in dense matter . Its application was shown  to remove most of the anomaly when applied to large data sets. Updating the fit termed ‘W’ in  for the present analysis of the 120 data points and including also the uncertainties due to the neutron densities, we find the following low-density limit values: $`b_0=0.009\pm 0.007m_\pi ^1`$ and $`b_1=0.088\pm 0.005m_\pi ^1`$, in full agreement with the free $`\pi N`$ values  of $`b_0^{\mathrm{free}}=0.0001_{0.0021}^{+0.0009}m_\pi ^1`$ and $`b_1^{\mathrm{free}}=0.0885_{0.0021}^{+0.0010}m_\pi ^1`$. The parameter Re$`B_0`$ assumes then the value $`0.02\pm 0.03m_\pi ^4`$, which is most acceptable (see Ref. ). Although successful on a phenomenological level, this approach has been deemed controversial. For example, Chanfray et al.  have recently discussed the additional medium effect of a $`\sigma `$ field on the renormalization of the pion field in dense matter, finding a considerably smaller renormalization effect than that suggested by Weise. On the other hand, Kolomeitsev et al.  have argued that Weise’s fairly substantial renormalization effect can be derived by considering the energy dependence of $`V_{\mathrm{opt}}`$ within a systematic chiral perturbation expansion and enforcing consistently gauge invariance. This latter suggestion will be addressed in a forthcoming paper.

We wish to acknowledge the assistance of G. Yaari in the calculation of nuclear densities. This research was partially supported by the Israel Science Foundation.