B G Giraud
email: giraud@dsm-mail.saclay.cea.fr

Lon-Chang Liu
email: liu@lanl.gov

Retrieving Nuclear Information from Protons Propagating through a Thick Target

Keywords:

The multiple scattering of high-energy particles in a thick target is formulated in an impact parameter representation. A formalism similar but not identical to that of Molière is obtained. We show that calculations of particle beam broadening due to multiple Coulomb scattering alone can be given in closed form. The focus of this study is on whether or not the broadening of the Coulomb angular distribution prevents the retrieval of nuclear-interaction information from measuring the angular distributions of charged particles scattered from a thick target. For this purpose, we study multiple scatterings with both the nuclear and Coulomb interactions included and we do not make a small-angle expansion. Conditions for retrieving nuclear information from high-energy protons propagating through a block of material are obtained.

Introduction and basic formalism

Understanding and calculation of the broadening of a particle beam when it propagates through a block of materials are important not only to make multiple-scattering corrections for cross-section measurements in physics experiments but also to many applications such as radiography by means of high-energy protons. Many different theories of the multiple scattering of electrons by thick targets have been formulated in the past sixty years [1]- [7]. The theory of Molière [1] [7] has received extensive attention because it gives the best agreement with data concerning the broadening of Coulomb angular distribution, arising from the multiple scattering of charged particles from the atoms in thick targets. An excellent and succinct derivation of the Molière theory was given by Mott and Massey [START_REF] Mott | The Theory of Atomic Collisions[END_REF].

As one deals with a large (almost astronomical) number of scatterers in a thick target, the numerical aspect of the calculation becomes extremely demanding. One can obtain good results provided that careful approximations are carried out. Over the years, successful parametrizations of the broadening of angular distributions due to Coulomb multiple scattering have been established [START_REF][END_REF]. For hadronic projectiles, such as protons, nuclear interactions also contribute to multiple scattering. However, this latter aspect has not yet received sufficient attention in the literature. In this work, our focus is, therefore, on effects of multiple scattering on angular distributions in the Coulomb-nuclear interference region and in the region where the nuclear interaction dominates. We formulate the multiple-scattering problem in such a way that not only in the case with Coulomb multiple scattering alone can analytical evaluations of the beam broadening become possible, but also the calculations are greatly facilitated when both nuclear and Coulomb interactions are taken into consideration.

For a very high energy (e.g., ≥ 20 GeV) proton scattering from a single nucleus, the Coulomb cross section decreases rapidly with the increase of scattering angles in such a way that the Coulomb cross section is negligible with respect to the nuclear cross section already at scattering angles as small as several milliradians. However, the Coulomb peak is rapidly broadened by proton-atomic nucleus multiple scatterings. Clearly, for the purpose of extracting the forward amplitude of the basic hadron-nucleus strong interaction, one should use targets as thin as possible and then employ Molière-type theory to correct for the Coulomb contribution from the measurement [10]. On the other hand, in many practical applications the thickness of the "target" is often fixed by specific needs, which is by no means thin. It becomes, therefore, interesting to know how much nuclear information can still be learned from hadrons scattered from a thick target. Certainly, the feasibility of learning nuclear information can bring added value to probing materials with hadronic beams. In other words, will the broadening of the Coulomb angular distribution render impossible any study of the nuclear cross sections? One naturally expects that the survival of nuclear information, if any, depends on the target thickness, i.e., on the number of nuclear scatterers that a proton encounters in a block of material. We use our formalism to examine this question.

In this section, after deriving the basic multiple-scattering formalism, we discuss the important bearing of an ability to retrieve nuclear information on many applications. In Section II, we show the broadening of angular distributions by Coulomb multiple scattering in closed form. The broadening of angular distributions by combined Coulomb and nuclear interactions are studied by means of semi-analytical models in regions of small momentum transfers (Section III) as well as large momentum transfers (Section IV). We find that it is possible to retrieve nuclear information from protons scattering from a thick target. Conclusions and suggestions are presented in Section V.

It is well known that high-energy elastic scattering is basically forward peaked, which allows to a very good approximation to neglect the longitudinal momentum transfer. Glauber [START_REF] Glauber | High-Enery Physics and Nuclear Structure[END_REF] has shown that, in an impact-parameter representation, every function O(q) of the tranverse momentum transfer q, whether O is an amplitude or a cross-section, can be parametrized in terms of a profile function φ(b) defined in the impact plane,

O(q) = (2π) -2 d b exp(i q • b) φ(b) = (2π) -1 ∞ 0 b db J 0 (qb) φ(b), (1)
where azimuthal symmetry is assumed, J 0 is the Bessel function of the first kind, and q and b are the moduli of the transverse momentum transfer q and impact parameter b, re-spectively. Nonessential complications (e.g. spins, etc) are here understood. Conversely, the profile function results from the inverse Fourier transform and is given by

φ(b) = d q exp(-i q • b) O(q) = 2π ∞ 0 q dq J 0 (qb) O(q). (2)
Without loss of generality, we consider a target which consists of one kind of nuclei.

For a thin target of thickness t and atomic density ρ the probability that a beam particle undergoes a scattering is p tot 1 = t ρ σ tot 1 , where σ tot 1 is the total cross section [START_REF] Joachain | Quantum Collision Theory[END_REF]. (The subscript 1 denotes the single scattering.) The transmission probability is, therefore, given by p trans

1 = 1 -p tot 1 .
From the definition of the differential cross section σ 1 ( q), one obtains the sum rule σ tot 1 = d q σ 1 ( q). The scattering probability density p 1 ( q) is related to the differential cross section by [START_REF] Joachain | Quantum Collision Theory[END_REF], [START_REF] Rodberg | Introduction to the Quantum Theory of Scattering[END_REF] p 1 ( q) = t ρ σ 1 (q). It is the probability that a particle experiences scattering in the direction q. The sum rule of σ 1 ( q) leads to the sum rule p tot 1 = d q p 1 ( q). Most often, only the modulus q counts, σ 1 ( q) = σ 1 (q). Hence,

σ tot 1 = π d(q 2 ) σ 1 (q).
In thick targets the beam can bounce forward from many nuclei and/or electronic clouds and/or different atoms. These multiple scatterings are incoherent because the scatterers are separated far apart with respect to the ranges of the screened Coulomb and nuclear interactions so that the scattering waves are already in the asymptotic region before the next collision occurs. Furthermore, the target is not crystalline on a macroscopic scale, and thus the distance between scatterers are largely random. One must, therefore, add probabilities (not amplitudes) coming from individual scatterings. If one splits the thick target with thickness T into a large number N of thin targets each with thickness T /N, then the total multistep probability for a particle to be transmitted without any scattering is

P trans M ≡ P 0 = (1 -σ tot 1 T ρ/N) N .
The differential probability density for just one scattering in this situation with many thin targets, each with thickness T /N, is

P 1 ( q) = N σ 1 ( q) (T /N) ρ 1 -σ tot 1 T ρ/N N -1 . (3)
This represents a combination of scattering from any one single layer and transmission through the remaining (N -1) layers without scattering. The factor N in front of the right hand side accounts for the N layers, obviously.

For double scattering one must count pairs of layers and fold two single-scattering probability densities p 1 . Hence,

P 2 ( q) = N(N -1) 2 T ρ N 2 σ 2 ( q) 1 - σ tot 1 T ρ N N -2 , (4)
where

σ 2 ( q) = d q ′ σ 1 ( q ′ ) σ 1 ( q -q ′ ). ( 5
)
For triple scattering, an identical argument induces the result

P 3 ( q) =   N 3   T ρ N 3 σ 3 ( q) 1 - σ tot 1 T ρ N N -3 , (6)
where

σ 3 ( q) = d q ′ d q ′′ σ 1 ( q ′ ) σ 1 ( q ′′ -q ′ ) σ 1 ( q -q ′′ ) (7)
is a double convolution. Again, the first factor is the counting of all triplets of layers.

It is useful at this stage to introduce the profiles

φ 1 (b) = d q exp(-i q • b) σ 1 ( q) = 2 π ∞ 0 q dq J 0 (qb) σ 1 (q), (8)
with φ 1 (0) = σ tot 1 , and, ∀n,

φ n (b) = d q exp(-i q • b) σ n ( q) = [φ 1 (b)] n . (9)
If we define Φ 1 as the profile of P 1 , then we see that,

Φ 1 (b) = N T ρ N φ 1 (b) 1 -σ tot 1 T ρ/N N -1 , (10)
and, more generally, the profile of P n reads,

Φ n =   N n   T ρ φ 1 N n 1 - σ tot 1 T ρ N N -n . (11)
The "total multistep profile" due to P M = N n=1 P n can be written as,

Φ M (b) = 1 - σ tot 1 T ρ N + T ρ φ 1 (b) N N -P 0 = 1 + ν [φ 1 (b)/σ tot 1 -1] N N -(1 -ν/N) N , (12)
where we have introduced the parameter ν = σ tot 1 T ρ. Like p tot 1 , ν is dimensionless. However, it is not a probability for a thick target. Indeed, because the mean free path, Λ, of a beam particle is 1/(ρσ tot 1 ) , therefore, ν = T /Λ represents the average number of collisions of a beam particle when it passes through a target of thickness T . Hence ν can be a very large number.

In the limit, N → ∞, we obtain,

Φ M (b) = exp ν φ 1 (b) σ tot 1 -1 -exp(-ν) . (13)
Here, we emphasize the nonlinear action of ν on Φ M (b). The second term in Eq.(13) comes from the limit of P 0 , it also shows that ν has the meaning of a beam decay rate in a target of thickness T . Note that this latter term is neglected in the Molière theory [1] as well as in its reformulation by Bethe [7] and by Mott and Massey [START_REF] Mott | The Theory of Atomic Collisions[END_REF]. Hence, we find a form similar but not identical to that of the Molière theory for the final, multistep probability density P M :

P M (q) = (2π) -2 exp(-ν) d b exp(i q • b) exp ν φ 1 (b) σ tot 1 -1 . (14)
Again, we note that P M (q) depends nonlinearly on ν. It is appropriate to underline the importance of the above-mentioned extra term in the present formulation. If T is small rather than large, T = t, then ν is also small. The presence of this extra term reduces P M to p 1 = t ρ σ 1 , as should be for thin targets. On the other hand,

φ 1 (b) → 0 as b → ∞.
Consequently, the first term in the curly bracket has a limit, exp[ ] → 1, that exactly compensates the second term in the curly bracket, ensuring the convergence of the integration for P M (q).

It is interesting to note from Eq. ( 14) that

P tot M ≡ d q P M ( q) = 1 -exp(-ν) = 1 -P M t , (15)
where P M t is the total transmission probability. This last equation is a sum rule for the multiple-scattering probability. Multiple-scattering differential cross sections can be related to probability densities by the general relation

P M (q) = T ρ σ M (q) , (16)
in the same way σ 1 ( q) is to p 1 ( q). In summary, three steps thus occur in this formalism: (i)

Fourier transform thin target data σ 1 into their profile φ 1 , see Eq. ( 8); (ii) Find ν = T ρ σ tot 1 and exponentiate T ρ φ 1 , see Eq. ( 13); (iii) Fourier transform Φ M back into a probability distribution P M , see Eq. ( 14).

It is reasonable to parametrize the single-scattering distribution as σ 1 = σ 1c +σ 1n with σ 1c being the Coulomb cross sections and σ 1n the sum of cross sections of nuclear scattering and nuclear-Coulomb interference. The separation of σ 1 into σ 1c and σ 1n induces the same for the profile function:

φ 1 = φ 1c + φ 1n . The relation exp(βφ 1c ) exp(βφ 1n ) -1 = [exp(βφ 1c ) -1] + exp(βφ 1c ) [exp(βφ 1n ) -1] (17)
with β = ν/σ tot 1 then leads to a split of P M as the sum of two probability densities,

P M c = exp(-ν) (2π) 2 d b exp(i q • b) exp ν φ 1c (b) σ tot 1 -1 (18)
and

P M n = exp(-ν) (2π) 2 d b exp(i q • b) exp ν φ 1c (b) σ tot 1 exp ν φ 1n (b) σ tot 1 -1 . ( 19
)
This allows a perturbative consideration of nuclear effects at those angles where Coulomb scattering dominates. Note that the exponent in Eq. ( 18) contains a denominator σ tot 1 and not σ tot 1c . Hence P M c is proportional to a pure Coulomb process with an effective value of ν, namely ν ef f = (σ tot 1c /σ tot 1 ) ν.

Alternately, at angles where nuclear scattering might dominate, the roles of φ 1c and φ 1n can be interchanged to generate similar formulae, namely,

P ′ M n = exp(-ν) (2π) 2 d b exp(i q • b) exp ν φ 1n (b) σ tot 1 -1 (20)
and

P ′ M c = exp(-ν) (2π) 2 d b exp(i q • b) exp ν φ 1n (b) σ tot 1 exp ν φ 1c (b) σ tot 1 -1 . (21)
One can parametrize the screened Coulomb interactions as σ 1c (q) = α C α (q 2 + κ 2 α ) -α . The powers α, screening momenta κ α , and normalizations C α are mainly functions of the charge Z of each individual nucleus. But σ 1n will depend on both Z and the mass number A. The global normalization of σ M will also depend on the target thickness or the parameter ν.

Hence, the theory is essentially driven by three parameters of a thick target, namely, Z, A and ν. Experimental measurements of σ M might, conversely, permit a determination of such three parameters when the nuclear nature of the target is unknown a priori, as is most often the case for radiographic studies where a recovery of A, besides Z, would be precious. Success will occur, however, only if multiple scattering does not spoil the information carried by A.

This question is the main concern of the following sections.

Broadening of angular distributions by multiple scatterings

The mean-square width q 2 1 of the distribution σ 1 represents a useful observable for the broadening of the cross-section distribution σ 1 and can be defined by an integral

q 2 1 = (σ tot 1 ) -1 d q q 2 σ 1 ( q) = 2 π (σ tot 1 ) -1 ∞ 0 dq q 3 σ 1 (q), (22)
if it converges. The use of Eq.( 8) and elementary properties of the Fourier transform allows us to write

q 2 1 = -lim b → 0 1 b d db b d db φ 1 (b) φ 1 (0) , (23)
where the operator between brackets ( ) comes from a two-dimensional Laplacian in cylindrical coordinate space. The same procedure gives the mean-square width of σ M ( q) as

q 2 M = -lim b → 0 1 b d db b d db Φ M (b) Φ M (0) , (24)
where Φ M (0) = (1e -ν ). Assume, for the sake of the argument, that σ 1 (q) is a Gaussian,

σ 1 (q) = C exp[-q 2 /(2κ 2 )], (25)
where C is a suitable normalization and the parameter 1/κ is the interaction range. For instance, if one discusses screened-Coulomb interactions, then 1/κ is of the scale of an atomic radius. Then one obtains σ tot 1 = 2 π C κ 2 , q 2 1 = 2κ 2 , and

φ 1 (b) = 2 π C κ 2 exp(-b 2 κ 2 /2). (26)
From Eqs. ( 13) and (24) one further obtains

Φ M (b) = exp ν e -b 2 κ 2 /2 -1 -exp[-ν] ( 27
)
and

q 2 M = ν q 2 1 1 -exp(-ν) . ( 28
)
Since ν is large in general, the denominator is ≃ 1. Consequently, the multiple scattering has broadened the mean-square width by a factor ν, as might be expected from a Brownian motion in the transverse-momentum space. The multiplication of q 2 1 by ν also occurs if we start from a "polynomial Gaussian distribution" (q/κ) 2n exp[-q 2 /(2κ 2 )]. This growth rate is very general and can be viewed as one more version of the central limit theorem. As additional evidence, one finds that if σ 1 (q) is of the functional form σ 1 (q) = C(q 2 +κ 2 ) -n with an exponent n > 2, then q 2 1 = κ 2 /(n -2) and again q 2 M = ν q 2 1 /(1e -ν ) ≃ ν q 2 1 .

In the following, we illustrate the broadening of the cross-section distribution in the case of a screened Coulomb scattering. We fit σ 1 at small angles by a few terms of the sum

σ 1c (q) = m>2 C m (q 2 + κ 2 m ) m , (29)
where m can be half-integers as well as integers, and C m and κ m are fitting parameters. It follows that

σ tot 1c = π m C m (m -1) κ 2(m-1) . (30)
For definiteness, we take two terms with m = 5/2 and 4, namely,

σ 1c (q) = C 5/2 (q 2 + κ 2 5/2 ) 5/2 + C 4 (q 2 + κ 2 2 ) 4 . (31)
Hence,

σ tot 1c = 2πC 5/2 3κ 3 5/2 + πC 4 3κ 6 4 . (32)
Dividing both sides by σ tot 1c , we obtain

1 = 2πC 5/2 3κ 3 5/2 σ tot 1c + πC 4 3κ 6 4 σ tot 1c ≡ a 5/2 + a 4 . (33)
Eq. (33) shows that both a 5/2 and a 4 are dimensionless numbers between 0 and 1. An advantage of using Eq. ( 29) is that its Fourier transform gives the profile function φ 1c in terms of analytical functions which can be easily analyzed, i.e.,

φ 1c (b) = 2πC 5/2 (1 + κ 5/2 b) exp(-κ 5/2 b) 3 κ 3 5/2 + πC 4 b 3 K 3 (κ 4 b) 24 κ 3 4 , ( 34
)
where K 3 is the modified Bessel function of the second kind. One verifies easily that φ 1c (0) = 2πC 5/2 /(3κ 3 5/2 ) + πC 4 /(3κ 6 4 ) = σ tot 1c . The "Coulomb" multistep profile then reads

Φ M c = exp ν a 5/2 (1 + κ 5/2 b)e -κ 5/2 b + a 4 b 3 8κ 3 4 K 3 (κ 4 b) -1 -e -ν . (35)
When b → ∞, it is easy to verify that Eq. ( 34) induces exponential decreases with ranges κ -1 5/2 and κ -1 4 . Numerical integrals with such integrands converge well. The final integral for the Coulomb cross section then reads

σ M c (q) = 1 2 π T ρ ∞ 0 db b J 0 (qb) Φ M c (b). ( 36
)
Let κ be an average between the two momenta κ 5/2 and κ 4 , which are both atomic scales. It is now convenient to scale momenta and lengths as q = κQ and b = B/κ. The dimensionless Q will be a few units or a few tens, if one wants to describe scattering angles moderately larger than the Coulomb peak. We also expect that the values of B contributing to the integral

σ M c (Q) = exp(-ν) 2πκ 2 T ρ ∞ 0 dB B J 0 (QB) × exp νa 5/2 1 + κ 5/2 κ B e -κ 5/2 B/κ + νC 4 B 3 8κ 3 K 3 κ 4 κ B -1 (37)
should be mainly between 0 and several units. In atomic units, all parameters κ, κ 5/2 /κ and κ 4 /κ are of order 1. It remains to estimate the dimensionless magnitudes of ν a 5/2 and ν a 4 /8.

From Eq.(33), a 5/2 and a 4 are moderate fractions of 1. It is thus the large number ν that drives the integrand.

It is also convenient to write

σ 1 (q) = κ -4 σ 1 (Q) , φ 1 (b) = κ -2 φ 1 (B) (38)
with σ 1 (Q) and φ 1 (B) being dimensionless. Because we work with systems of atomic scale, we further set κ to be 1, meaning that our primary scale is "atomic". In this scale, all lengths and momenta will, respectively, be given in units of atomic radius and its inverse.

To show the shrinking of profiles by multiple scatterings we plot φ 1c and Φ M c in Fig. 1 as the "crosses" and solid curves when σ 1c = (Q 2 + 1) -5/2 , and, respectively, as the "circles" and dashed curves when σ 1c = (Q 2 + 1) 3 Loss of nuclear information at small momentum transfers

In this section, we introduce a semi-realistic model for σ 1 (Q) which contains "nuclear" information. We will investigate (a) changes of normalizations and not just shrinking or dilation of shapes of σ(Q), and (b) how nuclear information may become lost. We will, therefore, illustrate the blurring of signal through the study of various relevant quantities, such as P M , P M c , P M n , P ′ M c , P ′ M n . We also use analytical models to ensure that the blurring comes from physics and is not a result of numerical imprecision. A good analytical model must satisfy the following constraints: (i) positivity of the sum of σ 1c and σ 1n ; (ii) big contrasts between maxima and minima; (iii) analyticity in both the momentum and the impact parameter representations; and (iv) significant differences between the atomic and the nuclear scales for profiles.

Let

σ 1c (Q) = (1 + Q 2 ) -2 (1 + Q 2 /100) -1 . ( 39
)
This is qualitatively realistic, because the factor, (1 + Q 2 ) -2 , represents a screened Coulomb scattering. The additional, artificial factor, (1 + Q 2 /100) -1 , is here just for the convergence of Q 2 1c . Then we further use, with σ 1n (0) = 0, the following semi-realistic σ 1n :

σ 1n (Q) = 11 Q 2 /65060976287632087746874800000000 × exp(1/2 -Q 2 /50) 24684880296681586800 e 12 25 (3600 -169 Q 2 + Q 4 ) 2 -175731507577476 (Q 2 -144) 2 (491485925 -20593402 Q 2 + 215573 Q 4 ) +40919125 e 119 50 (Q 2 -25) 2 (3973881778272 -103796385841 Q 2 + 463199137 Q 4 ) (40)
The quality of the model with respect to the requirements (ii) and (iv) are evidenced by Fig.

3.

As it is allowed by the split of σ 1 into a "Coulomb" part (σ 1c ) and a "nuclear" part (σ 1n ), our σ 1n can be positive or negative, as long as σ 1 remains positive. Our σ 1n was fine tuned to create four clear "nuclear" signals, namely two maxima of σ 1 near Q = 8 and 18, and, as signatures of interferences, two sharp minima at Q = 5 and 12. Furthermore, we adjusted its parameters so that the maxima do not exceed ∼ 1% of the forward peak of σ 1c . Note also that our model σ 1n has only two maxima and, thus, carries no "nuclear" information for Q > 40. This is designed to track whether or not the maxima, if they survive the blurring of angular distribution by multiple scatterings, would migrate towards larger values of Q. The log 10 σ 1 and log 10 σ 1c of our toy model are shown as functions σ(Q) in Fig.

It is trivial to deduce σ 1n visually.

The corresponding profiles read, in closed forms,

φ 1c (B) = 200 π -K 0 (B) 9801 + K 0 (10 B) 9801 + B K 1 (B) 198 , (41)
These profile functions are shown in Fig. 4. The width of φ 1n is significantly smaller than that of φ 1c , as one should expect when comparing a "nuclear" profile to an "atomic" one. A geometrical ratio of widths might be ∼ 10 -4 or even ∼ 10 -5 , but the model ratio we choose, between ∼ 1/5 and ∼ 1/10, is sufficient for a pedagogical study and much more convenient numerically. This choice of "data" gives, after a numerical implementation of Eq. ( 14), the total multistep probability distributions shown in Fig. 5. The solid curve is the same as that in Fig. 3, namely, log 10 σ 1 . The dashed, linked-crosses, and linked-circles curves represent log 10 P M for ν = 4, 9 and 16, respectively. The result is striking, on two counts: (i) the forward peak is more and more damped, the distributions extending more and more towards larger momenta, and (ii) the "nuclear information", whether minima or maxima, becomes rapidly blurred beyond recognition. Furthermore, the broadening of distributions does not seem to push much residual information towards larger momenta. The broadening process is also confirmed by the behavior of the component P M c of P M , shown in Fig. 6. In our model σ tot 1 ≃ 5.7 and σ tot 1c ≃ 3.0. We chose a large nuclear contribution, σ tot 1n ≃ 2.7, in order to emphasize nuclear effects. However, at Q < 40, even this exaggerated nuclear information did not survive multiple scatterings.

In Fig. 7, we show the various probability distributions P M n (Q). We note again that multiple scatterings wash away nuclear information. A similar feature is also seen in the P ′ M n given in Fig. 8. Besides the damping and information loss which are evident from Figs. 7 and8, we may stress a feature of Fig. 8, namely the transformation of "negative cross sections" into positive ones after multiple scattering. In order to create interferences, it was necessary, at the stage of making a model for σ 1n , to create negative values interfering with σ 1c . As has already been pointed out, this is allowed as long as σ 1 remains positive; there is a degree of freedom in modeling σ 1n . The solid curve in Fig. 8 shows log 10 |σ 1n |. One sees four arches, the first and the tiny third ones meaning negative values. Such "negative" arches disappear in the dashed curves representing P ′ M n . This disappearance justifies the use of models where σ 1n can be not everywhere positive as long as σ 1n + σ 1c is everywhere positive, as was discussed after Eq. ( 40).

An advantage of our use of special analytical forms for the cross sections, Eqs. (39,40), is that such forms induce analytical profiles, Eqs. (41,42), which in turn allow analytical forms for the multistep profiles, Eqs. ( 13) and (18-21). Values of Q 2 can then be easily obtained from the use of Eq. (24). The rates of broadening as functions of ν can also be readily calculated. Fig. 9 shows how, at values of ν smaller by several orders of magnitudes than those estimated from geometric cross sections, the square-widths Q 2 of P M n , P M c already increase linearly with ν. We have also noted a similar behavior of the widths of P ′ M n and P ′ M c . In summary, the signature of nuclear information (diffractive oscillations in the differential cross section) in the region of small momentum transfers is washed away by the broadening of the angular distribution. This happens even with our model that has exaggerated nuclear cross sections. In the next section, we examine if there exist momentum-transfer regions where the multiple scattering of the proton does not completely blurr nuclear signals. 4 Conditions for observing nuclear information Let r ≡ q 1/2 /q min1 be the ratio of the half-width of the Coulomb peak to the momentum transfer at which the first minimum due to nuclear diffraction is observed for a thin target. At high energies, both q 1/2 and q min1 occur at very small angles. Consequently, r = θ 1/2 /θ min1 with the θ's being the respective scattering angles corresponding to q 1/2 and q min1 . From Eq. (28), it is reasonable to expect a rule, q 2 M ≃ ν q 2 1 ≃ ν q 2 1/2 (= νr 2 q 2 min1 ), hence that there is a critical value ν crit ≃ r -2 , above which nuclear signals will be obliterated by the broadening of the Coulomb peak. In other words, nuclear signals can only be observed at q ≫ q min1 for ν > ν crit .

In Fig. 10, we show the elastic scattering differential cross sections of protons scattered from a thin 208 Pb target at ∼20 GeV, which we have calculated by using the method of optical model of Ref. [START_REF] Glauber | [END_REF] with a screened Coulomb interaction. The calculated cross sections exhibit the main characteristics of high-energy proton-nucleus scattering, namely, a narrow forward Coulomb peak and the diffractive oscillations at larger angles. Here, the diffractive pattern constitutes the nuclear signal. One notes that the first diffractive minimum lies at about ∼ 6 milliradians. We have noted from our calculation that θ 1/2 for the Coulomb peak is of order ∼ 0.003 milliradians. It is therefore reasonable to assume that, for high-energy proton scattering from nuclei, r is of order ∼ 10 -3 or less in general. At most one might consider r of order ∼ 10 -2 . Accordingly, although the geometric size of a nucleus is typically ∼ 10 4 -10 5 smaller than that of its atom, the range of "nuclear information profiles" at high scattering energies may be taken ∼ 10 2 to ∼ 10 3 smaller than the range of the screened atomic profile, and possibly much smaller. While the model used in the previous section where ∼ 0.5 < r <∼ 0.1 is pedagogically justified, an analysis with a smaller r is in order. We first introduce a model in which r = 10 -2 ; the profile function φ 1 (B) is the sum of a "Coulomb" term,

φ 1c (B) = B K 1 (B) (43)
and a "nuclear" term,

φ 1n (B) = 4 × 10 -6 1 + exp[800 (B -1/100)] . (44)
The profile φ 1c gives a bare Coulomb cross section of the form σ 1c (Q) ∝ (1 + Q 2 ) -2 and the Woods-Saxon profile φ 1n makes, in practice, a window with range r = 1/100 indeed. The coefficient 800 in its exponent creates a "smoothed" Heaviside function. Both profiles are normalized so that σ tot 1c = 1, and σ tot 1n /σ tot 1c = 4 × 10 -6 except for a negligible factor 1 + e -8 . This cross-section ratio is quite compatible with the r 2 suggested by Fig. 10. Hence, the set of parameters given in Eqs. ( 43) and ( 44) is more realistic than that used in the previous section.

The result for various angular cross sections σ M (Q; ν), compared with the single scattering The ν crit ∼ r -2 rule is also seen in the previous section, where the use of ∼ 0.5 < r <∼ 0.1 induces the loss of nuclear signal as early as ν >∼ 10. To further verify this rule, we use the same φ 1c but use instead r = 0.5 × 10 -2 in φ 1n , namely,

φ 1n (B) = 4 × 10 -6 1 + exp[1600 (B -1/200)] . ( 45
)
The results are shown in Fig. 12. As we can see, the observation of nuclear signals is much improved; in agreement with the rule that r = 1/200 elevates ν crit to a higher value, ∼ 4×10 4 . As a last test of the ν crit ∼ r -2 rule, we keep φ 1c (B) = BK 1 (B) and let

φ 1n (B) = 2 × 10 -8 1 + exp[8000 (B -1/1000)] , (46)
which is an obvious r = 1/1000 case. The results are shown in Fig. 13, where it is clear that, as expected, ν crit occurs between 10 -5 and 10 -6 .

Q Because ν = T ρ σ tot 1 ≃ T ρ σ tot 1c , the existence of a ν crit induces a critical target thickness T crit such that the retrieval of nuclear signal is possible for target thickness T sufficiently less than T crit ; namely,

T < T crit ≃ ν crit ρ σ tot 1c . (47)
For p-208 Pb elastic scattering at 20 GeV, σ tot 1c ≃ 6.7 × 10 8 mb = 6.7 × 10 -19 cm 2 . The density d and the atomic mass number A of lead are 11.3 g/cm 3 and 208, respectively. . In so far as Q = 1 corresponds to θ ∼ 0.003 milliradians, the survivor oscillation seen for ν = 10 5 in Fig. 13 between Q ∼ 6000 and Q ∼ 12000, compatible with the expected period ∼ 2π/r, would demand experimental measurements at angles of order a few dozens of milliradians at most.

The proton-nucleus Coulomb cross section

σ tot 1c is ∝ Z 2 R 2 e ∼ Z 4/3
, where Z is the target charge and R e is the root-mean-square radius of electric charge distribution in an atom with R e ∼ a 0 /Z 1/3 and a 0 being the first Bohr radius [START_REF] Gottfried | Quantum Mechanics[END_REF]. Hence, one can estimate σ tot 1c;pA for proton scattering from a given nucleus A at 20 GeV by using σ tot 1c;pA ∼ λ σ tot 1c;p-P b with the scaling factor λ = (82/Z) -4/3 . Hence, for a same ν one has T (pA) = T (pP b)ρ P b /(ρ A λ).

Results for a sample of atomic nuclei at ν = 10 5 are given in Table I. Of course, the price one pays in studying the nuclear signals that survive the multiple-scattering broadening is that one has to measure the angular distribution with good energy resolution at large proton scattering angles. In the case of 20 GeV incoming protons, the angles are about tens of milliradians, where the magnitudes of the cross sections are quite small. However, such measurements should be feasible with the currently available technology.

The main mathematical and physical statement proposed by this work about multiple scattering consists in folding probabilities rather than scattering amplitudes. This is justified by the incoherence which is expected between the different scatterers of a thick, non crystalline target. Simultaneously, an eikonal approximation, justified by the very high energy of the beam, allows a familiar impact parameter representation with profiles. Furthermore, small-momentum expansions [1] are not employed in the formulation. As a consequence of such initial statements, a Poisson process is found, leading to an elementary formalism of convolutions and exponentiations in a context of Fourier-Bessel transforms. This Poisson process is nothing but a random walk in transfer momentum space.

The central limit theorem is at work and the details of nuclear oscillations and interferences between Coulomb and nuclear scattering are blurred very fast as soon as the parameter ν = T ρ σ tot 1 , a measure of the number of collisions, exceeds a critical value of order r -2 . Here r is the ratio of the range of the nuclear profile to that of the atomic profile.

Below this critical value of ν, and at moderate and large momentum transfers (at the cost of very small elastic cross sections in the latter case) our conclusion is that some nuclear information remains observable. Such information is contained in oscillations of the multistep angular cross section σ M (Q) with periods ∼ 2π/r, oscillations that are similar to those of the Bessel function, J 1 (r Q), which typically represents pure nuclear diffractive scattering.

Our model analysis shows that the characteristic distances between successive cross-section maxima and minima in the angular distribution remain essentially unchanged while each of these oscillations dampens as ν increases.

From the point of view of retrieving nuclear signals from protons traveling through a thick target, for which no sufficient attention was given in the literature, our work is more of a general feasibility study rather than a specific numerical evaluation. We have made use of analytical and semi-analytical models to bring out the basic features of the underlying physics. We believe that the positive feasibility concluded from this study will sustain tests in detailed numerical applications.

One problem which has not been solved in the present work, however, is to find an estimate of the ν dependence of such periods ∼ 2π/r. Our numerical evidence suggests that the dependence is not very strong, despite all the causes for a broadening of the signal, but our models and calculations lack the precision needed to tabulate such periods into functions of ν. This effort is under consideration for an extension of the present work.

In summary, below ν crit , which is of order 1/r 2 with r being the ratio of the range of nuclear profile to the range of atomic profile, the nuclear signals can be retrieved from proton scattering from a thick target of thickness T < T crit . We suggest a conservative upper bound, namely, T ≤ 0.1 T crit for practical considerations.

- 4 .Fig. 1 .Fig. 2 .

 412 Fig. 1. Crosses: φ 1c (B) for σ 1c (Q) = (Q 2 + 1) -5/2 . Solid curve: the corresponding multistep Φ M c (B) if ν = 4. Circles and dashed curve: φ 1c (B) and Φ M c (B) (with ν = 4) for σ 1c = (Q 2 + 1) -4 . All profiles normalized to 1 at B = 0.

Fig. 3 .

 3 Fig. 3. log 10 σ 1 (solid curve) and log 10 σ 1c (dashed curve) as functions of Q. and, with u = -e 12/25 and v = -e 119/50 ,

Fig. 4 .

 4 Fig. 4. log 10 φ 1 (solid curve), log 10 φ 1n (dotted curve), and log 10 φ 1c (dashed curve) as functions of B.

Fig. 5 .

 5 Fig.5. The dependences on Q of log 10 σ 1 of single scattering (solid curve) and log 10 P M of multiple scattering when ν = 4 (dashed curve), 9 (linked-crosses curve), and 16 (linked-circles curve), respectively.

Fig. 6 .

 6 Fig. 6. Solid curve: log 10 σ 1c (Q). Dashed, linked-crosses and linked-circles curves: log 10 P M c (Q) for ν = 4, 9 and 16, respectively.

Fig. 7 .

 7 Fig. 7. Dashed, linked-crosses and linked-circles curves: log 10 P M n (Q) for ν = 4, 9 and 16, respectively.

Fig. 8 .

 8 Fig. 8. Solid curve: log 10 |σ 1n (Q)| for a nuclear signal. Dashed, linked-crosses and linked-circles curves: probability distributions log 10 P ′ M n (Q) when ν = 4, 9 and 16, respectively.

Fig. 9 .

 9 Fig.9. The dependences on log 10 ν of log 10 Q 2 for P M (Q) (solid curve), P M c (Q) (dashed curve) and P M n (Q) (dotted curve). Note that the slopes ≃ 1 when ν >> 10.

Fig. 10 .

 10 Fig. 10. Differential cross sections of p-208 Pb elastic scattering at 20 GeV.

σ 1 (Fig. 11 .

 111 Fig. 11. The dependences on Q of log 10 σ 1 (solid curve) and of log 10 σ M for ν = 2000 (dashed curve) and ν = 10000 (dotted curve). All with r = 10 -2 .

Fig. 12 .

 12 Fig. 12. The dependences on Q of log 10 σ 1 (solid curve) and of Y = log 10 σ M for ν = 2000 (dashed curve) and ν = 10000 (dotted curve). All with r = 0.5 × 10 -2 .

Fig. 13 .

 13 Fig.[START_REF] Rodberg | Introduction to the Quantum Theory of Scattering[END_REF]. The dependences on Q of log 10 σ 1 (solid curve) and of log 10 σ M for ν = 10 5 (dashed curve) and ν = 10 6 (dotted curve). All with r = 10 -3 .

  Hence, ρ = (d/A)N Avog = 3.27 × 10 22 cm -3 , with the Avogadro number N Avog = 6.02 × 10 23 [1/mole]. Because from Fig.10it is likely that r ∼ 10 -3 and our analysis indicates that

Table 1

 1 Target thickness T corresponding to ν = 10 5 for 20-GeV protons.

	Atom d [g/cm 3 ] Z	λ	A ρ [10 22 cm -3 ] T [cm]
	Pb	11.3	82 1.000 208	3.27	4.6
	Cu	8.9	29 0.250 64	8.39	7
	Al	2.7	13 0.086 27	6.02	29
	Mg	1.74	12 0.077 24	4.37	44
	Be	1.85	4 0.018	9	12.38	67

ν crit ∼ 10 6 , then Eq. (47) gives T crit ≃ 46 cm. Hence, a nuclear signal can be retrieved at ν ≤ 10 5 , which corresponds to T ≤ 4.6 = 0.1 T crit

[cm]


Adhitya Kamakshidasan

Spatial Comparison of Cricketers

Keywords:

There has been an increasing demand from the cricketing community to introduce newer metrics to analyze the game. Data from ball tracking and prediction have urged statisticians to re-look at existing comparative measures. Using methods from topological data analysis, we introduce a new technique to compare cricketers using spatial features. We use data from IPL seasons (2012 -2017) to compare our results with an existing ranking scheme.

INTRODUCTION

When players are hit with injuries in the Indian Premier League (IPL), teams generally scout for like-for-like replacements without changing team composition characteristics. Such replacements are often chosen by using aggregate benchmarks like Strike Rate and Economy, which fail to capture spatial information.

Spatial information associated with a player's performance is usually represented by a pitch-map (see Figure 2). Though such visualizations are keystone in understanding performance over a single match, their interpretation over multiple matches is limited, owing to the high number of deliveries. Cricket analysts are commonly interested in understanding strong and weak scoring regions of a player across the pitch. For example, experts often suggest that Suresh Raina struggles playing short pitched deliveries, despite being the highest aggregate run scorer across all IPL seasons.

Topology-based techniques are ideal for understanding such expert insights due to two main reasons --they can efficiently identify locations of interesting features over a spatial domain; and their ability to look at the overall distribution of values. From a cricketing perspective, such techniques helps in identifying the abilities of a player across different lines and lengths of deliveries. In a more broader sense, we propose that these techniques can also be used to compare and cluster players of different traits, for injury replacement strategies by team managements.

In this work, using spatio-temporal data and methods from topological data analysis, we present an approach to compare cricketers.

SPATIAL COMPARISON

Irrespective of the format, cricket is always played on a rectangular pitch that is 22 yards long and 10 feet wide. Every delivery in a match, is bowled on the same pitch and has various attributes associated with it. Some attributes include delivery type, pitch location, shot played, runs scored/conceded, players involved and scoring region. We believe that a player's utility can be directly determined based on their behaviour across the pitch. A key step in identifying such behaviour would be to locate regions corresponding to significant aspects [START_REF] Gudmundsson | Spatio-temporal analysis of team sports[END_REF].

A scalar function f : D → R maps points in a spatial domain D to real values. In this work, we are interested in the spatial region corresponding to a pitch, which is represented by a planar domain.

* e-mail: adhitya.kamakshidasan@inria.fr A function value is defined over each point in this planar domain with the goal to capture the spatial utility of a player. An example of a scalar function, that is used in this work, is the number of runs scored off each delivery. We choose this function for illustration, since it intuitively captures the scoring capability of a player over various locations across the pitch. A high function value at a given location implies easier scoring, thereby enhancing/decreasing the utility of batsman/bowler.

In order to efficiently compute the topological features of the scalar function f , it is represented as a piecewise linear (PL) function f : K → R. The planar domain D of the function is represented by a 2D regular grid K. The function is defined on the vertices of the grid and linearly interpolated across each box. To account for the high density of a point cloud representing a player, the set of points are first interpolated to a regular grid using scattered data interpolation techniques and then a scalar function is computed for each player.

The critical points of a smooth real-valued function are exactly where the gradient of a scalar field becomes zero. The critical points of PL function are classified based on the behaviour of the function within a local neighbourhood. An important measure that is often associated with critical points is its persistence. Persistence diagrams encode the persistence of features as points in a plane.

Understanding player performance across regions of the pitch, is akin to understanding the persistent features of super-level components of the scalar field. Pitch regions having high persistence indicate easier scoring ability. Here, the specific high and low function values of a scalar field are not too important, but rather their distribution is key for comparing the behaviour of two players. We capture this, by finding similarity between persistence diagrams.

For two scalar fields, the similarity of persistence diagrams is generally measured by using a distance function like the Wasserstein distance or the Bottleneck distance. Low persistent features can be considered to be topological noise in the scalar field [START_REF] Tierny | The Topology ToolKit[END_REF]. In our case, two players can be compared to one another by their respective persistence diagrams after removal of low persistent features.

We also compare the players using Earthmovers distance. This distance reflects the minimal amount of work that must be performed to transform one distribution into the other by moving "distribution mass" around. For this, grayscaled pixel values of the scalar field are used to construct histograms for similarity detection between players.

METHOD

Hawk-Eye is a commercial vendor that captures real-time spatiotemporal data for assisting officials with decision making processes. We demonstrate our method using Hawk-Eye data scrapped from various sources on the internet. The pair (x, y) specifies the pitching location of a delivery on a cricket pitch, where y is the horizontal coordinate and x is the vertical coordinate relative to an origin at the centre of off-stump for a right-handed batsman. The positive y-axis points to the right from the bowler's perspective and the positive x-axis points towards the non-strikers stumps.

Our analysis includes ball-by-ball data from all IPL matches across 2012 and 2017. In this work, we look only at pitching locations and runs scored.

For a given interpolation method and a distance measure, the scalar field of one player can be compared to all other players in the consideration set, to obtain a unique rank descriptor for the player. The rank descriptor of a player showcases similarity/dis-similarity between all other players. For a given player, the descriptor would rank the player most similar (himself) as numero uno, and rank the most dissimilar player as equivalent to the number of players in the consideration set. This can be extended to obtain a rank matrix for all players (see Figure 4) wherein each row represents the rank descriptor for a player. For the construction of the rank descriptor, either Bottleneck/Wasserstein distance measures are computed between the persistence diagrams or an Earthmovers distance measure is computed between the respective scalar fields. We compare our results with Cricmetric's ranking system. Cricmetric provides a non-parametric way of ranking players by computing an eigenfactor score (EFscore). A similar rank matrix is constructed for Cricmetric by averaging EFscores across seasons followed by the L 1 norm. One common way to compare ranking schemes is by using the Pearson Correlation Coefficient (PCC). This measures the linear relationship between two data-sets. Like most other correlation coefficients, this one also varies between -1 and +1 with 0 implying no correlation. We compare each player's rank descriptor with that of Cricmetric's by using PCC. Table 1 reports the median correlation values for each ranking scheme across players in the comparison set.

RESULTS

Our method identifies six different clusters of players from Figure 4. We asked a cricket expert to help us identify and group several characteristics of these players to interpret our results. The player types found along the diagonal of the rank matrix, are listed from top-left to bottom-right as follows -"game closers", "slow-starting quickscorers", "hit-or-miss players", "multiple role players", "stroke makers", and "extremely aggressive starters".

From Table 1 we can clearly see that the correlation values with Cricmetric are not high. Since there is no "ground truth" to ranking systems, we cannot argue that our method performs better/worse than another one. However, the table indicates a shortcoming in that, our method is highly susceptible to the interpolation technique and the distance measure used for comparison.

Figure 1 :

 1 Figure 1: Critical points of a scalar field f defined on a PL 2-manifold (left) and its persistence diagram D( f ) (right). In the persistence diagram, each pair of critical points are represented by a vertical bar (yellow) and their persistence is given by height of the bar.

Figure 2 :

 2 Figure 2: Sunil Narine's bowling across six IPL seasons. Pitch-map visualization for runs conceded function (first row). Delaunay triangulation of pitch locations with z-coordinate as scalar function value (second row). Scalar field after Shepard's interpolation and removal of topological noise respectively (third, fourth rows).

Figure 3 :

 3 Figure 3: Grayscale pixel distribution of Sunil Narine (left) and Lasith Malinga (right) used by Earthmovers distance for player comparison.

Figure 4 :

 4 Figure 4: Rank matrix for batsmen obtained using Shepard's interpolation and Earthmovers distance. Each batsman on the left is ranked with a batsman on the top. Similar batsmen are highlighted in green.

Table 1 :

 1 Median correlation values -Spatial Ranking vs Cricmetric

	Distance Measure	Interpolation Method Sibson's Spline Kriging Shepard's
	Earthmovers	0.279	0.301	0.153	0.229
	Bottleneck	0.16	0.094	0.094	0.028
	Wasserstein -1	-0.185	0.112	0.078	-0.168


B G Giraud
email: giraud@dsm-mail.saclay.cea.fr

A Weiguny
email: weiguny@uni-muenster.de

Coordinates, modes and maps for the density functional

Special bases of orthogonal polynomials are defined, that are suited to expansions of density and potential perturbations under strict particle number conservation. Particle-hole expansions of the density response to an arbitrary perturbation by an external field can be inverted to generate a mapping between density and potential. Information is obtained for derivatives of the Hohenberg-Kohn functional in density space. A truncation of such an information in subspaces spanned by a few modes is possible. Numerical examples illustrate these algorithms.

I. INTRODUCTION

The well-known Hohenberg-Kohn (HK) density functional [1] and its finite temperature generalization by Mermin [2] suffer from the absence of constructive algorithms after their respective existence theorems. The Thomas-Fermi (TF) approach, however, and related developments such as [3], [4], have gone a long way into creating functionals with practical physical values. For reviews on the effectiveness of the detailed forms of the functional found empirically, see for instance [START_REF] Dreizler | Density Functional Theory[END_REF] and [START_REF] Perdew | A Primer in Density Functional Theory[END_REF]. For applications of Skyrme forces to nuclear densities, see for instance [7] and [8].

Standard perturbation theories (particle hole hierarchy of excitations, configuration mixing, generator coordinates, etc.), extrapolating from well understood mean field theories, give a constructive approach to the intricacies of a true ground state (GS), at the well known heavy cost of calculations with many degrees of freedom. But such theories proved to be practical, because suitable truncations were found that restricted calculations to few modes, collective or not, subspaces with fewer degrees of freedom. The purpose of the present note is to attempt answering a similar question in the space of densities rather than the space of wave functions: are there possible truncations, is there a possibility to restrict the functional to a set of few density modes?

For this we visit again the fundamentals of the HK functional F [ρ] in a systematic approach, based upon the following chain of arguments, i) given the full Hamiltonian, H = T + V + U, with a fixed kinetic operator T = i t i , a fixed two-body potential operator V = i>j v ij , and a variable one-body potential operator U = i u i , assume a non degenerate, square normalized GS Ψ, with its corresponding eigenvalue E, density ρ and functional F [ρ] ≡ Ψ|(T + V )|Ψ = E -Ψ|U |Ψ ; find the functional derivatives δρ/δu and, considering first F as a functional of u rather than ρ, find δF/δu, ii) expand such functional derivatives into suitable bases, to describe them by convenient matrices and vectors, iii) then invert δρ/δu to know δu/δρ, iv) furthermore obtain δF/δρ by eliminating δu between δF/δu and δu/δρ; further information about F might be obtained by integrating δF/δρ, or by comparing with phenomenological approaches, such as gradient expansions, v) at each stage, try a compression of the information, by a truncation of the theory to a few "density modes". A preliminary question is in order, however: can this formal program be carried if particle number is conserved in the mean only, as occurs with Lagrange multiplier techniques? According to [START_REF] Dreizler | Density Functional Theory[END_REF], the chemical potential, as a function of a continuous particle number, shows derivative discontinuities. We thus find it safer, in this paper, to stick to "slices" of the functional, those for fixed, integer particle numbers. We also restrict our considerations to pure eigenstates of H, at zero temperature.

In Section II, we carry the program in the absence of V ; the trivially soluble situation of independent fermions allows us to easily describe the mapping, ρ ↔ u, for both infinitesimal and finite variations of u. In Section III we reinstate V, but use the Hartree-Fock (HF) approximation to still obtain ρ without excessive technical complications. Section IV is dedicated to a better understanding of the "tangent" mapping δρ/δu [9] when full correlations are present. Then Section V introduces, via a new family of orthogonal polynomials, candidates for density and potential space modes, that might allow a compacted description of the functional and its ρ ↔ u mapping. An investigation of the relevance of such modes is provided numerically, for a toy model of independent fermions. The numerical investigation is continued in Section VI, by means of a second toy model, with now correlated fermions. Section VII contains a discussion and a conclusion.

II. PARTICLE-HOLE EXPANSIONS FOR INDEPENDENT FERMIONS

Let Z be the particle number for a finite Fermion system. For simplicity we ignore discrete labels such as spins and isospins. Set v = 0, temporarily, the case of two-body forces being discussed later. The one-body potential function u is taken here as a local potential u(r). The Hamiltonian then boils down to H = Z i=1 (t i + u i ) and its GS, assumed to be non degenerate, is a Slater determinant. The HK functional, in this special case, reduces to the kinetic expectation value, Ψ|T |Ψ .

Consider a perturbation δu of the local potential u. We find it practical to expand it in an orthonormal basis of functions w α (r), namely δu(r) = ∞ α=1 w α (r) δu α . According to the HK theorem [1], this basis must be orthogonal to a "flat potential component" w 0 (r) = 1. This is satisfied if we find a basis such that, ∀α > 0, dr w α (r) = 0. It is then understood that the index α will run from 1 to ∞. For obvious practical reasons, however, the expansion will sooner or later be truncated at some finite rank of the basis.

The perturbation δu induces a perturbation δρ of the density of the GS Ψ of H, and we find it convenient to expand δρ in the same orthonormal basis {w β (r)}, namely δρ(r) = ∞ β=1 w β (r) δρ β . The fact that this basis satisfies the constraint, ∀β > 0, dr w β (r) = 0, is very useful because δρ does not change the particle number, namely δρ automatically satisfies the condition dr δρ(r) = 0. It is very convenient that the constraint of particle number conservation and that of "non constant potential variation" allow the same choice for our forthcoming basis {w α }.

A trivial particle-hole argument then provides that perturbation δΨ induced by δu. Let the "hole index" i = 1, ..., Z and "particle index" I (running from Z + 1 to ∞) denote occupied and empty orbitals, respectively, with the corresponding single particle energies η i and η I and orthonormal wave functions ψ i and ψ I . For the sake of simplicity in the following, such orbital wave functions ψ(r) are assumed to be real in the coordinate representation. Each filled orbital picks a variation δψ i (r) = I ψ I (r) I|δu|i /(η iη I ). Hence,

δρ(r) = 2 iI ψ i (r)ψ I (r) I|δu|i η i -η I . (1)
This reads, when δu and δρ are expanded,

δρ β = 2 iIα D βiI I|w α |i η i -η I δu α , (2)
where D denotes the projection of a particle-hole product of orbitals upon the basis {w α },

D βiI = dr w β (r) ψ i (r) ψ I (r). (3)
Note, incidentally, that particle-hole orthonormality ensures that, ∀iI, dr ψ i (r)ψ I (r) = 0. Hence, functions w α expanded in the basis of particle-hole products ψ i ψ I , a basis to be orthonormalized, automatically fulfill the requested condition for δu and δρ. Furthermore, positivity of the density is guaranteed as variations δρ in Eq. ( 1) are based on variations of the wave function, in particular by particle-hole admixtures to the GS determinant Ψ. Define the matrix,

N βα = 2 iI D βiI I|w α |i η i -η I . ( 4
)
Notice also that the perturbation matrix element, I|w α |i , coming from δu, is nothing but an integral of a three term product, amounting to D αiI . Notice finally that the energy denominators correspond to a propagator G = Q(E -QHQ) -1 Q, if Q is the particle-hole space projector. This operator is diagonal in the particle-hole space, obviously. Then, in a condensed notation, N = 2 D G D, where the tilde denotes transposition between index α and pair index iI. Consider now δρ and δu as just vectors with components δρ β and δu α , respectively, sooner or later truncated. Since Eq. ( 2) reads δρ = N δu, the HK theorem states that, under the usual condition of non degeneracy for the GS of H, an inversion is possible. Namely for any δρ which leaves ρ + δρ in the manifold of actual densities, there exists a unique u + δu, provided δu does not add a constant component to u. Under such precautions, the infinite matrix N can be inverted, and the same can be expected under "reasonable" truncations of N . Accordingly, while Eq. ( 2) provides the functional derivative δρ/δu, one obtains the functional derivative δu/δρ,

δu α = β N -1 αβ δρ β . (5)
Now, that variation δE of the GS energy induced by δu is trivial. It just reads,

δE = Ψ|δU |Ψ , (6)
because of the stationarity of the GS energy with respect to δΨ. Accordingly, for the functional under study,

δF = δE -δ Ψ|U |Ψ = -δΨ|U |Ψ -Ψ|U |δΨ = -2 Ψ|U |δΨ = -2 Ψ|U G δU |Ψ , (7)
hence,

δF = -2 iI i|u|I (η i -η I ) -1 I|δu|i , (8)
and finally, with proper expansions,

δF = -2 iI i|u|I (η i -η I ) -1 α D αiI β N -1 αβ δρ β . (9)
It may be convenient here to set from u a column vector U with components

U iI = dr u(r) ψ i (r) ψ I (r), (10)
hence the functional derivative δF/δρ reads, in a condensed, matrix and vector notation,

δF = -Ũ G D D G D -1 δρ. (11)
This simplifies if one observes that, because of orthogonality between particle and hole orbitals, the product ψ i (r) ψ I (r) can be expanded in the w-basis as,

ψ i (r) ψ I (r) = β D iIβ w β (r). (12)
Accordingly,

U iI = dr u(r) β D iIβ w β (r) = β u β D iI β , (13)
with u β the components of u in our special basis. Combining Eqs. ( 11) and ( 13) results in

δF = - αβγ u β D G D βγ D G D -1 γα δρ α = - α u α δρ α . (14)
This avoids the transition between different bases through the matrix D. One thus recovers the trivial result, δF = dr u(r) δρ(r), but it must be kept in mind that, here, u has become a functional of ρ. Whether this simplification is made or not, this makes a set of numerical, non linear, coupled, partial differential equations relating F and ρ.

The non linearity comes in particular from the orbitals and single particle energies which occur in the definition of N . We stress again that the vector, Ũ D -1 , just makes an "α" representation of u, converted from its particle-hole representation U.

A comment about Legendre transforms is here in order [9]. According to the Hellmann-Feynman theorem, δE/δu = ρ. But then, F ≡ Edr u(r) ρ(r) is nothing but the Legendre transform of E and the primary degree of freedom is not u any more, but ρ. Note that the reasoning remains if V is reinstated. In all cases, u is recovered from δF/δρ = -u.

For Eq. ( 11), and its generalization if two-body forces are present, to become a tool to obtain information about F, dynamical models are obviously necessary. These are the subjects of several of the forthcoming sections.

III. TWO-BODY FORCES AND HARTREE-FOCK MODEL

In this section, we stay with Z fermions, but reinstate in the Hamiltonian the two-body interaction v ij with the operator V = Z i>j=1 v ij . The HK functional is Ψ|(T + V )|Ψ . While a Slater determinant Ψ was available as the true GS of a simpler H = T + U in the previous section, we cannot usually obtain the true GS with two-body forces present in a full H = T + V + U. Thus, in this section, we tolerate for Ψ the HF ground state of H, with energy E 0 , and furthermore assume that this HF approximation does not create degeneracies between distinct Ψ's. Under this precaution of uniqueness, there exists an extension of the HK theorem. Indeed, in the space of Slater determinants, let Ψ and Ψ ′ be the HF GSs of H = T +V +U and H ′ = T +V +U ′ , respectively, and let ρ(r) and ρ ′ (r) be their respective densities. The two Hamiltonians differ by their (local) one-body operators U = Z i=1 u(r i ) and

U ′ = Z i=1 u ′ (r i ) only. Their HF GS energies E 0 = Ψ|H|Ψ and E ′ 0 = Ψ ′ |H ′ |Ψ ′ ,
non degenerate, may be equal or distinct. Now, if ρ and ρ ′ were equal, then the usual HK arguments, namely

E ′ 0 -E 0 < Ψ|(H ′ -H)|Ψ = dr [u ′ (r) -u(r)] ρ(r) and E 0 -E ′ 0 < Ψ ′ |(H -H ′ )|Ψ ′ = dr [u(r) -u ′ (r)] ρ(r)
, necessarily lead to ρ = ρ ′ , by contradiction. It can be stressed here that, again because of the stationarity of the energy with respect to variations of Ψ, we can still take advantage of Eq. ( 6) for the variation of the energy induced by a variation δu. This reads, with notations already used in the previous section,

δE 0 = Ψ|δU |Ψ . ( 15
)
The same holds every time we approximate the GS by means of the Rayleigh-Ritz principle in a restricted space of wave functions. As a general consequence, we obtain again Eq. ( 7), namely, δF = -2 Ψ|U |δΨ , for every such variational approximation of Ψ. That variation of δΨ induced by δu is slightly more complicated, in the HF case, than in the trivial case of Section II where H = T + U. Indeed, each filled orbital is driven by the perturbed HF equation,

[η i + δη i -u(r) -δu(r)] [ψ i (r) + δψ i (r)] + h2 2m ∆ r [ψ i (r) + δψ i (r)] = Z j=1 dr ′ v(r -r ′ ) [ψ j (r ′ ) + δψ j (r ′ )] 2 [ψ i (r) + δψ i (r)] - Z j=1 dr ′ v(r -r ′ ) [ψ j (r ′ ) + δψ j (r ′ )] [ψ i (r ′ ) + δψ i (r ′ )] [ψ j (r) + δψ j (r)] . (16)
The non locality of the HF mean field, because of antisymmetrization, is written in an explicit way in the right-hand side above, in the coordinate representation. An equivalent form of this perturbed HF Eq. ( 16) is obtained if we retain its first order terms only and consider the particle-hole infinitesimal components δc iI , again assumed here to be real numbers,

(η i -η I ) δc iI -I|δu|i = jJ [ IJ|v|ij + Ij|v|iJ ] δc jJ . (17)
Thus ψ i (r) becomes ψ i (r) + I ψ I (r) δc iI . Notice that δη i drops out from the calculation, as should be expected.

Then define in particle-hole space the symmetric matrix, with antisymmetrized matrix elements of v,

A (iI)(jJ) = (η i -η I )δ ij δ IJ -IJ|v|ij -Ij|v|iJ . ( 18
)
Here δ is a Kronecker symbol and we must use pairwise indices (iI) when defining the inverse A -1 to be used; this A -1 generalizes the propagator used in Eq. ( 4), and thus,

δc iI = jJ A -1 (iI)(jJ) J|δu|j . ( 19
)
This leads to the variation δρ, and the analog of Eq. ( 2) reads,

δρ β = 2 (iI)(jJ)α D βiI A -1 (iI)(jJ) J|w α |j δu α , (20)
where the overlap matrix D is the same as defined in the previous section. Actually, this boils down to the even simpler formula, in matrix and vector notations,

δρ = 2 D A -1 D δu, (21)
where the tilde again denotes transposition of that connection D between the particle-hole products ψ i (r)ψ I (r) and their rearrangement into an orthonormal basis {w α (r)}. Note, incidentally, that, if v = 0, the matrix F ≡ 2 D A -1 D boils down to the matrix N , which is obviously negative semidefinite. We even expect N to be negative definite. The same is expected for D A -1 D. The stability of our HF solutions is assumed as long, at least, as v is a weak enough interaction, and this "definiteness" of F is intuitively most likely.

In the following, we shall also need the inverse of F . The final result for the variation of F reads,

δF = -Ũ A -1 D D A -1 D -1 δρ, (22)
and, like in Section II, this expression, in a transparent notation, reads δF =α u α δρ α . For obvious reasons of numerical convergence, the number of needed w α states must be large enough to overlap a sufficient number of particle-hole components of δρ. But as will be found in the coming numerical applications, a surprisingly small number of w α states might sometimes be sufficient.

IV. SYMMETRY OF THE DENSITY-POTENTIAL MAPPING IN GENERAL

With the exact ground energy E and exact GS Ψ of a full H = T + V + U, and Q = 1 -|Ψ Ψ| the projector out of Ψ, the Brillouin-Wigner perturbation theory gives the exact result for first order functional derivatives,

|δΨ = Q E -QHQ δU |Ψ = ∞ n=1 |Ψ n Ψ n |δU |Ψ E -E n . (23)
We assume here that a resolution of the identity with real numbers and reasonable truncations, convenient for numerics, are available. The sum over excited states Ψ n includes integrals over the continuum, if necessary. Let us single out the first of our identical particles and integrate out all the other ones, to define the following transition densities, Θ n (r) = dr 2 dr 3 ... dr Z Ψ n (r, r 2 , r 3 , ..., r Z ) Ψ(r, r 2 , r 3 , ..., r Z ).

(

) 24
Notice that, from its very definition, Θ n integrates out to 0, namely dr Θ n (r) = 0. Hence Θ n can be represented in the w-basis without any loss of information. Since δU = Z i=1 δu(r i ) is a symmetric operator, it is clear that Eq. ( 23) also reads,

|δΨ = ∞ n=1 |Ψ n Z dr δu(r) Θ n (r) E -E n = n α |Ψ n Z dr w α (r) Θ n (r) E -E n δu α , (25)
where we have again expanded δu in the basis {w}. There pops out a matrix,

D α n = Z dr w α (r) Θ n (r), (26)
as a generalization of the matrix D α iI . Now, by definition, the density of the GS is,

ρ(r) = Z dr 2 dr 3 ... dr Z [Ψ(r, r 2 , r 3 , ..., r Z )] 2 , (27)
and its variation is, δρ(r) = 2Z dr 2 dr 3 ... dr Z Ψ(r, r 2 , r 3 , ..., r Z ) δΨ(r, r 2 , r 3 , ..., r Z ).

This becomes, if one replaces δΨ by its expression, Eq. (25),

δρ(r) = 2Z n α Θ n (r) D α n E -E n δu α . (29)
An expansion of δρ in the {w} basis gives its coordinates,

δρ β = 2 n α D β n 1 E -E n D α n δu α , (30)
hence, in an obvious notation, a symmetric "flexibility" matrix F = 2 D G D connecting δu and δρ. In hindsight, the symmetry of F (and of its approximations under the Rayleigh-Ritz variational principle) is straightforward. Indeed, since u = -δF/δρ, then δu α /δρ β =δ 2 F/(δρ α δρ β ). All denominators E -E n being negative definite, the negative definite nature of this exact F is also transparent. We conclude this section on the general case with explicit expressions for δF/δρ and δ 2 F/(δρδρ ′ ). With Eqs. ( 25), (26) the general form for δF, see Eq. ( 7), reads

δF = -2 Ψ|U |δΨ = -2 nα Ψ|U |Ψ n (E -E n ) -1 D αn δu α = -2 nα U n (E -E n ) -1 D αn δu α . ( 31
)
Here the numbers

U n = dr Ψ| Z i=1 u(r i ) |Ψ n = Z dr u(r) Θ n (r), (32)
are now the components of u in the space of transition densities, generalizing Eq. ( 10) for states Ψ, Ψ n containing correlations. Upon inverting Eq. ( 30) we find, as a generalization of Eq. ( 22),

δF = -Ũ G D D G D -1 δρ. (33)
This simplifies if we expand

Θ n (r) = Z -1 β D nβ w β (r). (34)
Then

U n = dr u(r) β D nβ w β (r) = β u β D nβ , (35)
with again the components of u in our special w-basis, u β = dr u(r) w β (r). Accordingly,

Ũ G D γ = βn u β D βn G nn Dnγ = 1 2 β u β F βγ , (36)
hence finally, as expected,

δF = - αβγ u β F βγ F -1 γα δρ α = - α u α δρ α , (37)
as before in sections II and III. Similarly, the second derivative of F is found directly from the inverse of Eq. (30),

δ 2 F/(δρ α δρ β ) = -F -1 αβ = -2DG D -1 αβ . (38)
Hence, from Eqs. (33) and (38),

δF/δρ α = 2 β Ũ G D β δ 2 F/(δρ β δρ α ), (39)
another useful equation to test phenomenological functionals F [ρ], by calculating quantities such as ŨG D from microscopic wave functions and energies for simple systems.

V. ONE DIMENSIONAL TOY MODEL, SPECIAL POLYNOMIALS Assume that r is just one dimensional, running from -∞ to ∞. Define t = -d 2 /(2dr 2 ), with a nucleon mass m = 1, h = 1 and p = -id/dr. In the present section, we are first interested in the Hamiltonian H 0 = Z i=1 (p 2 i /2 + r 2 i /2); it is not a bad approximation to most shell model Hamiltonians, whether one considers one-body potentials only or HF solutions to problems with two-body potentials as well. A trivial scaling of coordinates and momenta allows us to reduce to the case, ω = 1, any situation, Z i (p 2 i /2 + ω 2 r 2 i /2), where the physical spring constant would be different. Then we shall consider the functional F (ρ) = Ψ|H 0 |Ψ for a family of additional one-body potentials u, with the corresponding GS density ρ(r) of

H = Z i=1 [ t i + r 2 i /2 + u(r i ) ]. (40)
Set temporarily u = 0, namely consider H 0 and its GS density ρ 0 (r). Since ω = 1, which will be understood from now on, both initial particle and hole orbitals ψ k (r) are just trivial products ϕ k of a Hermite polynomial, a common Gaussian and a suitable normalization, ϕ k (r) = π -1 4 e -1 2 r 2 P k (r). For the sake of illustration, we list here the first five Hermite polynomials, with their coefficients adjusted for orthonormalization,

P 0 = 1, P 1 = 2 1 2 r, P 2 = 2 -1 2 (2r 2 -1), P 3 = 3 -1 2 r (2r 2 -3), P 4 = 2 -1 6 -1 2 (4r 4 -12r 2 + 3). ( 41
)
Particle-hole products, ϕ i (r) ϕ I (r), make, in turn, just polynomials again, now multiplied by e -r 2 . To build our basis, {w α (r)}, it is tempting to orthonormalize the set {ϕ i ϕ I } containing that Gaussian, e -r 2 , and recover forms 2

1 4 ϕ k (r √
2), with Hermite polynomials again, compressed by the obvious √ 2 for their argument r, because of the new factor, e -r 2 . This is correct for odd parity functions. But, for even parity ones, the constraint, dr w α (r) = 0, would be violated. Hence, out of each even function, 2

1 4 ϕ 2k (r √ 2)
, k > 0, we subtract a term proportional to 2

1 4 ϕ 0 (r √
2), letting the subtraction cancel the integral, ∞ -∞ . (Alternately, we considered all elementary functions r 2k e -r 2 .) Then we use a Gram-Schmidt process to reorthonormalize such subtracted states. Notice that the subtraction cancels out the polynomial state of degree zero, and therefore the transformation from Hermite polynomials to this new set of orthonormal polynomials is not unitary, but only isometric, with "defect index" 1. In other words, our basis has codimension 1. This is also clear from the degree 2 of the lowest member of the new even basis. For an illustration, we list the first four even states obtained,

w 2 (r) = 2 2 1 4 (2 r 2 -1)/ √ 3 π -1 4 e -r
√ 35 π -1 4 e -r 2 . ( 42
)
For the sake of comparison with Hermite polynomials, which rather go with a factor e -1 2 r 2 , we perform the transformation, r → r/ √ 2 on Eqs. ( 42) and multiply the results by a factor 2 -1 4 to retain their (ortho)normalization. Discarding norm coefficients from the resulting polynomials we get,

Q 2 = r 2 -1, Q 4 = 2 r 4 -7 r 2 + 1, Q 6 = 4 r 6 -32 r 4 + 47 r 2 -11, Q 8 = 8 r 8 -116 r 6 + 438 r 4 -453 r 2 + 39. ( 43
)
But, as already noticed, products ϕ i ϕ I carry a factor e -r 2 and we find it natural, in the following, to stick to those polynomials trivially derived from Eqs. (42),

Γ 2 (r) = 2 (2 r 2 -1)/ √ 3 , Γ 4 (r) = 2
and so on. We generated such polynomials up to degree 100 and will send them to interested readers. Such polynomials are orthonormal under the metric weight, e -2r 2 2/π. They must be completed by odd Hermite polynomials, P 2k+1 (r √ 2), suitably adjusted for the same metric. Hence, for instance,

Γ 1 (r) = 2 r , (45a) Γ 3 (r) = r (4r 2 -3) 2/3 , (45b) Γ 5 (r) = r (16r 4 -40r 2 + 15)/ √ 30 , (45c)
Γ 7 (r) = r (64r 6 -336r 4 + 420r 2 -105)/(6 √ 35) . ( 45d
)
More technicalities on such polynomials Γ and related polynomials can be found in [10].

It is then trivial to calculate both even and odd blocks, respectively, of the initial matrix D, see Eq. ( 3), according to the parity of the subscript of w and its associated polynomial Γ. With due normalizations, this reads,

D βiI = ∞ -∞ dr Γ β (r) e -r 2 (2/π) 1 4 P i (r) e -1 2 r 2 π -1 4 P I (r) e -1 2 r 2 π -1 4 . ( 46
)
For instance, if the hole label is restricted to i = 0, and the particle label I runs from 1 to 3, the sets of non vanishing odd, respectively even, matrix elements boil down to,

D 101 = 2 -3 4 π -1 4 , D 103 = - √ 3/ 4 (2π) 1/4 , D 303 = 2 -7 4 π -1 4 , D 202 = √ 3 2 -7 4 π -1 4 . (47)
We found it useful to precalculate and store such initial matrix elements D for the particle index I running up to 100 and the α index running up to 100 also. This fastens generic calculations of D when u becomes finite, as one represents (t + r 2 /2 + u) by a matrix on the oscillator basis, diagonalizes it with eigenvalues η k and orthonormal eigenvectors X ℓk , and finally expands orbitals of both holes and particles as ψ k (r) = ℓ X ℓk ϕ ℓ (r) in the same basis.

In [10] we set Z = 4, considered u to be an infinitesimal δu in the neighborhood of u = 0, then calculated and diagonalized the functional derivative N = δρ/δu. The eigenvectors of N defined density and potential infinitesimal perturbations having the same shapes. Now we set again Z = 4, but are rather interested in cases where u is finite. We are concerned in particular with the mapping between u and ρ, in that representation provided by the "modes" w α . Truncations at a maximum degree N are necessary. The finite expansion,

u(r) = N α=1 u α w α (r), (48)
defines those processed perturbations u. Given u, it is trivial to diagonalize H with a good numerical accuracy and obtain ρ. Then it is easy to obtain "coordinates in density space,"

ρ α = ∞ -∞ dr w α (r) [ρ(r) -ρ 0 (r)]. (49)
The harmonic potential, r 2 /2, serves here as the origin in potential space, and the origin in density space is the corresponding density ρ 0 . We show in Figures 1 and2, respectively, a grid of values {u 2 , u 4 } in potential space and its image grid of density coordinates {ρ 2 , ρ 4 }. Dots at grid corners help matching the object and the image. In this calculation, all coordinates u α have been set to vanish, except u 2 and u 4 , but it must be stressed that the resulting density variation, ρρ 0 , has non vanishing coordinates ρ 6 , ρ 8 , ... besides ρ 2 and ρ 4 . Such additional coordinates are small, but not very small, as shown by the grid for ρ 6 , ρ 8 in Figure 3. Qualitatively, if u contains one mode w α only, then ρ β tends to decrease when |α -β| increases. But this is likely to be valid for small enough u's only, in a linear response regime. Curvature effects, evidenced by Figs. 2 and3, must be expected further.

Of interest are plots of F in ρ-space. If u has two components u 2 , u 4 only, assume that F is a function of ρ 2 , ρ 4 only. Then a gradient ∇F (ρ) can be observed directly. In Figure 4, the 3D plot of F shows slight deviations from a traditional paraboloid. This is even more visible in Figure 5, showing the vector field {u 2 , u 4 }(ρ 2 , ρ 4 ), namely -∇F. The field, read from Figs. 2 and1, focuses towards the origin in ρ-space, but with clear distortions. We know that the field has a vanishing curl; it can be integrated back into F.

It is also trivial to create an approximate F in the following way: i) assume indeed that F depends only on ρ 2 and ρ 4 for this toy model, ii) take a few exact (numerical, actually) values of F at random points taken from the partner grids shown in Figs. 1 and2, iii) set a simple parametric ansatz such as,

VI. CORRELATIONS, FROM ANOTHER TOY MODEL

In the previous sections we skirted around the difficulty of obtaining a GS with true correlations. Now we shall mix several Slater determinants, each made of Z = 4 harmonic oscillator, one dimensional orbitals taken from h 0 = 1 2 (p 2 + r 2 ). The Hamiltonian is a complete one,

H = 4 i=1 1 2 - d 2 dr 2 i + r 2 i + u 2 w 2 (r i ) + u 4 w 4 (r i ) -V a 4 i>j=1 [δ(r i -r j -R a ) + δ(r i -r j + R a )] . (52)
This contact, finite range attraction between particles is expected to create a reasonable amount of correlations and was chosen to allow an easy precalculation and tabulation of matrix elements of v. Such matrix elements are again understood to be antisymmetrized. The first Slater determinant, Φ 0 , in the mixture contains the lowest Z orbitals of the harmonic oscillator. It is expected to make the dominant component of the configuration mixture Ψ, as we shall keep u 2 , u 4 within the grid seen in Fig. 1, and also the strength V a moderate. Let ξ i and ϕ j , i, j = 1, ..., Z be the orbitals of two Slater determinants Ξ and Φ, respectively. Define the cofactors C ij and double cofactors C ikjl of the determinant of scalar products ξ i |ϕ j . Such cofactors are very simple in the present orthogonal basis of orbitals, obviously. Then the matrix elements needed for the Hamiltonian matrix and the calculation of the HK functional read,

Ξ|(H 0 + U )|Φ = ik C ik ξ i |(h 0 + u)|ϕ k , Ξ|V |Φ = 1 4 ijkl C ikjl ξ i ξ j |v|ϕ k ϕ l . (53)
With V a = 5 and R a = 1 the grids shown by Figures 7 and8 come from a calculation with a single particle basis made of the first 10 harmonic oscillator orbitals and a corresponding Slater basis of 61 states made of Φ 0 and all positive parity one-particle-one-hole and two-particle-two-hole determinants built upon Φ 0 . Because of V the centers of the density grids are not at the origin defined by Φ 0 , obviously. Indeed, for u 2 = u 4 = 0, this calculation gives {ρ 2 , ρ 4 , ...ρ 14 } = {-.14, -.36, -.43, -.06, .04, -.01, .005} as the coordinates of the ground state density shift ρρ 0 and the ground state energy is E = -6.7, significantly down from Φ 0 |H 0 |Φ 0 = 8. Had we taken as the origin in density space the density of the HF solution for u = 0, different drifts of grid centers would have been observed. Such new drifts are likely to make better signals of true correlations in Ψ. For the sake of comparison between the first and the second toy models we kept ρ 0 as a reference, but we have a direct access to the amount of true correlations: it is easy here to calculate the density matrix ρ, the diagonal of which gives ρ. Here, for u = 0, the trace of ρ -ρ2 is of order 5%, a reasonable amount. At that grid corner, u 2 = u 4 = -2.2, the trace even reaches 8%. It can be concluded that this second toy model does create correlations.

While ρ 2 and ρ 4 both vary by ∼ .3 across their grid and ρ 10 , ρ 12 ... can be neglected, the variation of ρ 6 and ρ 8 across the grid is of order ∼ .1, which is not so small. For the sake of comparison with Fig. 4 we now show in Figure 9 a plot of the HK functional in terms of {ρ 2 , ρ 4 } again, but it will not be not forgotten in our agenda to find those two orthonormal combinations {σ 2 , σ 4 } of ρ 2 , ρ 4 , ρ 6 and ρ 8 which allow the "flattest" projection of the grid. It is clear that, in the spirit of Eqs. (50) and (51), the best parametrization of F should now be in terms of {σ 2 , σ 4 }. In any case, from Fig. 9, the minimum of F occurs at {ρ 2 , ρ 4 } = {-.14, -.36}, the grid center, as should be. Fig. 9 also shows deviations from a paraboloid clearly stronger than those of Fig. 4.

VII. SUMMARY, DISCUSSION AND CONCLUSION

Particle number conservation is essential to that key element in the proof of the HK theorem, the one-to-one correspondence between density and external potential. However, variations of the potential should not be trivial: they should differ from constants. We treated both constraints of i) matter conservation and ii) non triviality of potential variations on the same footing: a vanishing average for both the potential and the density variations.

This constraint of vanishing average was implemented by means of a new family of orthogonal polynomials; hence appeared a set of modes in both the density and the potential spaces. We proved, numerically with toy models, that such modes might have a physical meaning, on two counts, i) converse linear responses δρ/δu, δu/δρ might be reasonably simple when described in terms of such modes, and ii) the HK functional itself might be practically truncated into projections into subspaces spanned by a few modes.

We have not discussed in this paper the constraints of positivity of the density, but it is clear that, within an algebra of polynomials such as ours, positivity conditions are not too difficult to implement. We have not discussed either more subtle constraints related to the Sobolev nature of the topological spaces available for densities. For this question, we refer to [11], [12], [13]. It can be stressed again that the compatibility of truncations of densities into a finite number of "polynomial modes" with such fine constraints can easily be tested.

An ultimate goal would be to create a constructive theory of the HK functional. The functional differential equation, δF/δρ = -u[ρ], cannot be integrated in the density space as long as u[ρ] is not known accurately enough. Because of our detour through many-body perturbation theories we are clearly far from the goal, but this work gives a frame in which the task should become easier. The detour might allow a compression of the needed information through, for instance, the parametrization of a limited set of matrix elements of our matrices D for mean field approximations or D in the correlated cases. Our main results are i) the existence of those special, orthogonal constrained polynomials Γ α and associated modes w α which design convenient sets of coordinates and convenient parametrizations of F, δF/δρ, δ 2 F/(δρ δρ ′ ), ... etc., ii) the explicit relation we showed between these modes and the traditional perturbation theories used in the many-body problem and iii) the likely possibility of truncated descriptions and accurate parametrizations.

Actually, in the context of extended systems, the idea of density waves as important modes of the system has always been present. There remains to be seen, obviously, if our modes can be generalized to infinite systems. It also remains to be seen whether, for finite or infinite systems, our truncations are always justified, whether infinite resummations are possible, whether collective degrees of freedom are present, or absent, because of our new representation. Also, because of our choice of a Gaussian weight for the new family {Γ}, the present results are better meaningful if restricted to nuclei. A generalization to atoms and molecules obviously demands other weights for the constrained polynomials.

The usual perturbation theories have their hierarchy of modes in many-body space, most often a hierarchy of particlehole components. Our approach, tuned to the one-body nature of the density functional, replaces the particle-holes by other modes, in a transparent way.

It is a pleasure to thank Y. Abe, R. Balian and B. Eynard for stimulating discussions.

FIG. 1 .FIG. 2 .FIG. 3 .FIG. 4 .

 1234 FIG.1. Grid of parameters u2, u4 for the potential u = u2 w2 + u4 w4 used in the toy model.

FIG. 6 .

 6 FIG.6. Contours for Fapp.

FIG. 7 .FIG. 8 .FIG. 9 .

 789 FIG. 7. Second toy model: influence of the two-body force V on the ρ2, ρ4 image of the grid of Fig. 1. Compare with Fig. 2.

comes from fitting 26 values for Z = 4. Figure 6 shows several resulting contours, the smallest of which locates the minimum of F app very slightly only away from the origin, that is the true minimum by the very construction of the toy model. The value of the functional at the minimum turns out to be 8.0005165, instead of strictly 8. This toy numerical exercise demonstrates the possibility of contracting the description of the functional to few degrees of freedom.