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

K Katō
email: kato@nucl.sci.hokudai.ac.jp

A Ohnishi
email: ohnishi@nucl.sci.hokudai.ac.jp

S M A Rombouts
email: stefan.rombouts@rug.ac.be

Existence of Density Functionals for Excited States and Resonances

We show how every bound state of a finite system of identical fermions, whether a ground state or an excited one, defines a density functional. Degeneracies created by a symmetry group can be trivially lifted by a pseudo-Zeeman effect. When complex scaling can be used to regularize a resonance into a square integrable state, a DF also exists.

The aim of density functional (DF) theory is to contruct a functional that provides the energy expectation value for a correlated many-body state as a function of the one-body density, such that minimization of the DF leads to the exact ground state (GS) density. Since the existence theorem proven for GSs by Hohenberg and Kohn (HK) [1], its extension by Mermin [2] to equilibrium at finite temperatures, and the further development by Kohn and Sham (KS) [3] of an equivalent, effective, independent particle problem, a considerable amount of work has been dedicated to generalizations such as spin DFs [4], functionals taking into account the symmetries of the Hamiltonian [5], calculations of excited state densities [6,7], treatments of degeneracies or symmetries of excited states [8,9] and quasiparticles [10]. For the reader interested in an even more complete reading about both basic questions and applications, we refer to [11]- [19].

DFs for resonant states have received much less attention. We want to study this problem here. First we will address two related issues, namely that of a unified theory for ground and excited states and that of a theory for non degenerate and degenerate ones. A generalized existence theorem can be constructed by modifying the Hamiltonian in such a way that the spectrum is shuffled but the eigenstates are left unchanged, and by making a systematic use of the Legendre transform (LT) for a detailed analysis of the density.

A reminder of the HK proof is useful here. Consider a finite number A of identical fermions, with a † r and a r their creation and annihilation operators at position r, and the physical Hamiltonian,

H = T + V + U, where T = A i=1 t i , V = A i>j=1 v ij and U = A i=1 u i are
the kinetic, two-body interaction and one-body potential energies, respectively. For simplicity, we consider such fermions as spinless and isospinless and work at zero temperature. Both v and u may be either local or non local. Next, embed the system into an additional one-body, external field, W = A i=1 w i , to observe its (non linear!) response. The Hamiltonian becomes

K = H + W. It is understood that w is local, r|w| r ′ = w(r) δ( r -r ′ )
, although a DF theory with non-local potentials exists [20]. The usual Rayleigh-Ritz variational principle, where |ψ is just an A-particle, antisymmetric, square normalized, otherwise unrestricted wave function, applied to F M = min ψ F, with F = ψ|K|ψ , generates ψ min , the exact GS of K, with the exact eigenvalue F M . The minimum is assumed to be non degenerate, smooth, reached. Clearly, ψ min and F M are parametrized by w. An infinitesimal variation δw triggers an infinitesimal displacement δψ min , with δF M = ψ min |δW|ψ min . There is no first order contribution from δψ min . Define the one-body density matrix in coordinate representation, n( r, r ′ ) = ψ min |a † r a r ′ |ψ min . Its diagonal, ρ( r) = n( r, r), is the usual density deduced from |ψ min | 2 by integrating out all particles but one. Since δF M = d r ρ( r) δw( r), then δF M /δw( r) = ρ( r). Freeze t, v and u and consider F M as a functional of w alone. The HK process then consists in a Legendre transform of F M , based upon this essential result, δF M /δw = ρ. This LT involves two steps: i) subtract from F M the functional product of w and δF M /δw, i.e. the integral d r w( r) ρ( r), leaving F M = ψ min |H|ψ min ; then ii) set ρ, the "conjugate variable of w", as the primary variable rather than w; hence see F M as a functional of ρ. Step ii) is made possible by the one-to-one (1 ↔ 1) map between w and ρ, under precautions such as the exclusion of trivial variations δw that modify w by a constant only, see for instance [16] and [19]. The 1 ↔ 1 map is proven by the usual argument ad absurdum [1]: if distinct potentials w and w ′ generated ψ min and ψ ′ min (distinct!) with the same ρ, then two contradictory, strict inequalities would occur, d r

[ w( r) -w ′ ( r) ] ρ( r) < F M -F ′ M , and, d r [ w( r) -w ′ ( r) ] ρ( r) > F M -F ′
M . An inverse LT returns from F M to F M , because δF M /δρ = -w. Finally, the GS eigenvalue E 0 of H obtains as E 0 = min ρ F M [ρ]; the GS wave function ψ 0 of H is the wave function ψ min when w vanishes; that density providing the minimum of F M is the density of ψ 0 .

Consider now any excited bound eigenstate ψ n of H, with its eigenvalue E n . Then, trivially, ψ n is a GS of the semipositive definite operator (H -E n ) 2 . Since E n is not known a priori, consider rather an approximate value E n , obtained by any usual technique (configuration mixing, generator coordinates, etc.) and assume that E n is closer to E n than to any other eigenvalue E p . Then ψ n is a GS of (H -E n ) 2 . The possible degeneracy degree of this GS is the same whether one considers

H or (H -E n ) 2 . Introduce now K = (H -E n ) 2 + W.
If there is no degeneracy of either ψ n or its continuation as a functional of w, then the HK argument holds as well for K as it does for K. Hence a trivial existence proof for a DF concerning ψ n . But most often, ψ n belongs to a degenerate multiplet. Degeneracies are almost always due to an explicitly known symmetry group of H. Notice however that the external potential w does not need to show the same symmetry; hence, in general for K, there is no degeneracy of its GS; a unique ψ min emerges to minimize the expectation value of K. However, for that subset of zero measure in the space of potentials where w shows the symmetry responsible for the degeneracy, and in particular for the limit w → 0, precautions are necessary. Consider therefore an (or several) additional label(s) g sorting out the members ψ ng of the multiplet corresponding to that eigenvalue

(E n -E n ) 2 of (H -E n ) 2 .
There is always an operator G related to the symmetry group, or a chain of operators G j in the reduction of the group by a chain of subgroups, which commute with H and can be chosen to define g. For simplicity, assume that one needs to consider one G only. Then define g as an eigenvalue of G and assume, obviously, that the spectrum of G is not degenerate, to avoid a reduction chain of subgroups. It is obvious that, given some positive constant C, and given any chosen γ among the values of g, there is no degeneracy for the GS of (H 2 , even if w has the symmetry. When several labels become necessary with a subgroup chain reduction, it is trivial to use a sum j C j (G j -γ j ) 2 of "pusher" terms. Finally a DF results, now from the HK argument with K. We stress here that pusher terms, because they commute with H, do not change the eigenstates of either H nor (H -E n ) 2 . Only their eigenvalues are sorted out and reorganized. Note that the pusher expectation value vanishes for ψ nγ . Naturally, when w is finite, eigenstates of K differ from those of K, but what counts is the information given by the DF when w vanishes. A simplification, avoiding cumbersome square operators H 2 , is worth noticing. Consider the operator, K = H + C (G -γ) 2 + W. At the limit where w vanishes, there is always a choice of a positive constant C which makes the lowest state with quantum number γ become the GS. This leads to a more restricted density functional that is of interest for the study of an yrast line.

-E n ) 2 + C (G -γ) 2 . Nor is there a degeneracy of the GS of K = (H-E n ) 2 +C (G-γ) 2 +W = K+C (G-γ)
That DF, F M [ρ], based upon K, provides the expectation value,

F M [ρ] = ψ min | [ (H -E n ) 2 + C (G - γ) 2 ]
|ψ min , where ψ min , square normalized to unity, is also constrained by the facts that ψ min |a † r a r |ψ min = ρ( r) and K|ψ min = ε|ψ min for the eigenvalue ε = F M . It may be interesting to find a DF that provides the expectation value of H itself. This can be done by taking the derivative of F M [ρ] with respect to E n , at constant ρ. We suppose that this derivative exists, which is the case for a discrete spectrum at least. With the notation | ψ = d|ψ /d E n , and using the fact that ψ min |W| ψmin + ψmin |W|ψ min = w( r) dρ( r)/d E n d r = 0, one can write:

dFM [ρ] d En = 2 ψ min | E n -H |ψ min + ψmin |(ε -W)|ψ min + ψ min |(ε -W)| ψmin = 2 ψ min | E n -H |ψ min . (1)
Therefore we can define a new DF, = 0 for ψ min = ψ nγ . Hence the functional F D [ρ] is stationary at the exact density ρ = ρ nγ . It is not expected to be minimal at ρ nγ , however, unless the resulting eigenstate corresponds to the absolute GS when w vanishes.

F D [ρ] = E n - dF M [ρ] 2 d E n , (2)
Resonances may be defined as special eigenstates of H if one uses an argument à la Gamow, allowing some radial Jacobi coordinate r ≥ 0 to show a diverging, exponential increase of the resonance wave function at infinity of the form exp(ipr), where the channel momentum p is complex and ℑp < 0. It is well known that those eigenvalues E n describing resonances are complex numbers, with ℑE n < 0. There have been extensive discussions in the literature about the physical, or lack of, meaning of such non normalizable wave functions and about the wave packets which might be used to replace them, [21,22,23,24]. The point of view we adopt in this note is based upon the Complex Scaling Method (CSM) [25,26,27,28]: a modest modification of H transforms narrow resonances into square integrable states; then there is no difference between the diagonalization for bound states and that for resonances. The cost of the CSM, however, is a loss of hermiticity: the CSM Hamiltonian H ′ is non hermitian, somewhat similar to an optical Hamiltonian [25,26,27,28].

Given the ket eigenstate equation, (H ′ -E n )|ψ n = 0, where |ψ n is now a square integrable resonance wave function, we can consider the hermitian conjugate equation,

ψ n |(H ′ † -E * n ) = 0.
Clearly, ψ n is a GS, as both a ket and a bra, of the hermitian and semipositive definite operator,

Q exact = (H ′ † -E * n ) (H ′ -E n )
, with eigenvalue 0. Applying the same argument as before, but now to Q exact instead of (H -E n ) 2 , demonstrates the existence of a DF around the targeted resonant state.

In practice we do not know E n exactly. Given a sufficiently close estimate . For E n = E n the functional will be stationary at the density of the exact resonant state. While provinding a proof of existence, the construction of the exact functional for H ′ requires the knowledge of the exact eigenvalue E n . This might be an inconvenient limitation but fortunately calculations of numbers such as E n are usually much easier and much more precise than calculations of wave functions ψ n and/or of their densities.

E n of E n , an approximate GS eigenvalue |E n -E n | 2 occurs for Q apprx = (H ′ † - E * n ) (H ′ -E n ), at first order with respect to ∆Q = Q apprx -Q exact . Since ψ n is not a ket eigenstate of H ′ † = H ′ -2i ℑH ′ ,
If the resonance has good quantum numbers (QN)s inducing degeneracies, the same pusher terms as those which have been discussed above can be added to create a unique GS, from the operator,

Q exact + C (G -γ) 2 .
The HK argument, implemented with the full operator,

K ′ = (H ′ † -E * n ) (H ′ -E n ) + C (G -γ) 2
+ W, then proves that DFs exist for those resonances regularized by the CSM. Notice, however, that a simplified theory, with an "yrast suited" operator K ′ , linear with respect to H ′ , is not available here, since the restoration of hermiticity forces a product H ′ † H ′ upon our formalism.

We now consider a special case of rather wide interest in nuclear and atomic physics. i) Good parity of eigenstates of H 0 = T + V or H = H 0 + U when u is restricted to be even, is assumed in the following. Hence our eigendensities, quadratic with respect to the states, have positive parities. ii) We also assume that the number of fermions is even. iii) The QNs in which we are interested in this Section are the integer angular momentum L and magnetic label M of an eigenstate Ψ LM of H, where it is understood that the two-body v and one-body u interactions conserve angular momentum. When w is switched on and is not rotationally invariant, eigenstates of K, K, or K may still tolerate such labels LM by continuity. First, consider w = 0. The density ρ LM comes from the product Ψ * LM Ψ LM , but it does not transform under rotations as an {LM } tensor. Rather, it is convenient to define "auxiliary densities",

σ λ0 ( r) = L M=-L (-) L-M L -M L M | λ 0 ρ LM ( r),
where L -M L M | λ 0 is a usual Clebsch-Gordan coefficient. Each function σ λ0 ( r) now behaves under rota-tions as a {λ0} tensor. It can therefore be written as the product of a spherical harmonic and a radial form factor, σ λ0 ( r) = Y λ0 (r) τ λ (r) = (2λ + 1)/4π L λ (cos β) τ λ (r), where L λ is a Legendre polynomial and the angle β is the usual polar angle, counted from the z-axis. Conversely,

ρ LM ( r) = 2L λ=0 (-) L-M L -M LM |λ0 Y λ0 (r) τ λ (r). (3)
This provides a "Fourier analysis" of ρ LM in angular space. The density is parametrized by scalar form factors, τ λ . Since L is here an integer and furthermore ρ L-M and ρ LM are equal, and since Clebsch-Gordan coefficients have the symmetry property

L M L ′ M ′ | λ M ′′ = (-) L+L ′ -λ L ′ M ′ L M | λ M ′′ , then necessarily τ λ = 0 if λ is odd.
There are thus (L + 1) scalar functions, τ 0 , τ 2 , ..., τ 2L , to parametrize (L + 1) distinct densities ρ L0 , ρ L1 , ..., ρ LL . Because of the quadratic nature of the density observable, the even label λ for angular "modulation" of ρ runs from zero to twice L, with a "2L cut-off"; a signature, necessary if not sufficient, for an "L-density". Reinstate now w as the LT conjugate of ρ LM . It makes sense to study situations where w is restricted to expansions with (L + 1) arbitrary scalar form factors, w( r) = 2L even λ=0 Y λ0 (r) w λ (r). With inessential factors such as (-) L-M L -M L M | λ 0 omitted for simplicity in the following, every pair {r τ λ , r w λ } is conjugate. An eigendensity of K, K, K may have an infinite number of multipole form factors, but, with such restricted potentials w, only τ 0 , τ 2 , ..., τ 2L are chosen by the LT relating F M and F M .

It can make even more sense to restrict w to one multipole only, w( r) = Y λ0 (r) w λ (r), with λ = 0, or 2, ... or 2L, to study each multipole of ρ separately. For simplicity we now use the easier version of the theory, with that operator K which is suited to the yrast line. Add therefore to H a pusher term Z LM leaving intact the eigenstates, namely

Z LM = B [ L• L-L(L+1)] 2 +C (L z -M ) 2 .
Hence

K LMλ = H + Z LM + W λ = T + V + U + Z LM + W λ .
Here the subscript λ specifies that w is reduced to one multipole only. Then L is the total angular momentum operator and L z is its third component. This operator This note offers theorems for the existence of exact DFs for every excited bound state, and even narrow resonances, and every set of good QNs used in nuclear, atomic and molecular physics. Furthermore, the densities used as arguments of our DFs do not need to be fully three-dimensional ones; they can be radial form factors of multipole components of the states under study. Our existence theorems, though, suffer from the usual plague of the field: constructive algorithms are missing and empirical approaches will have to be designed. What is the corresponding (KS) theory [3]? In its usual form, the task of calculating the kinetic energy part of the DF is actually left to the solution of a Schrödinger equation, and this can be trivially generalized to any one-body part. Published studies of the KS formalism are actually dedicated to calculations of the functional derivative, δV xc /δρ( r), of the exchange and correlation part of the DF, coming from the two-body part V of the DF. Our present use of modified Hamiltonians, or even squares of H, introduces two-body operators, but also three-and four-body operators. For the versions where no squares of H occur, see the yrast suited operator K and Section 3, the nature of the three-and four-body terms, typically coming from ( L. L) 2 , is not forbidding, because of obvious factorization properties. Hence a KS theory might be realizable for such simplified versions. With squared Hamiltonians, however, a KS theory seems out of reach at present. A systematic analysis of solvable models on a basis of "modes" [19], however, may help to extrapolate such models into practical rules. For the discussion of differentiability, representability and fine topological properties of the w-and ρ-spaces, we refer again to [16]. Up to our understanding of the topology of the variational spaces, flat or curved [29], of general use in nuclear, atomic and molecular theory, the validity domain of our existence theorems is quite large. We have not used the time dependent formalism, although much progress has been made in deriving excitation energies from it [30]. A generalization of our arguments to finite temperatures seems plausible, however, and insofar as inverse temperature may be viewed as an imaginary time, a generalization to a time dependent theory is not excluded.

Z
B.G.G. thanks T. Duguet and J. Meyer for stimulating discussions, and also the hospitality of the Hokkaido University and Kyoto Yukawa Institute for part of this work. S.M.A.R. thanks D. Van Neck for interesting discussions.

  such that F D [ρ] = ψ min |H|ψ min and F D [ρ nγ ] = E n for the density ρ nγ of the eigenstate ψ nγ of H at energy E n . Furthermore one finds that δFD δρ [ρ nγ ] = 0, because δ ψmin| δρ H|ψ min + ψ min |H δ|ψmin δρ

  it is also perturbed at first order in ∆Q. Still one can copy the construction for F D [ρ], see Section 3, if one interprets the operator d/d E * n as d/dℜ E n + id/dℑ E n . The resulting functional F D [ρ] is linear in H ′

  LM moves the eigenvalues of H so that the lowest eigenstate of H with quantum numbers {LM } becomes the GS of H+Z LM . The commutator [H, Z LM ] vanishes indeed, and given A, t, v and u, there are always positive, large enough values for B and C that reshuffle the spectrum such that the lowest {LM } eigenstate Ψ LM becomes the GS of H + Z LM under this Zeeman-like effect. We stress again that Z LM changes nothing in the eigenfunctions, eigendensities, etc., of all our Hamiltonians if w is rotationally invariant. Furthermore, angular momentum numbers remain approximately valid for eigenstates of K LMλ if w is weak, and the same numbers might still make sense as labels by continuity when stronger deformations occur. Then the usual ad absurdum argument generates a map w λ ↔ τ λ , where τ λ (r) is the form factor of the λ-multipole component of the GS density for K LMλ , leading to an exact DF, for every {LM } lowest state and every even λ between 0 and 2L. A generalization to operators K LMλ , involving (H -E n ) 2 , is trivial.


B G Giraud
email: giraud@spht.saclay.cea.fr

Constrained orthogonal polynomials

We define sets of orthogonal polynomials satisfying the additional constraint of a vanishing average. These are of interest, for example, for the study of the Hohenberg-Kohn functional for electronic or nucleonic densities and for the study of density fluctuations in centrifuges. We give explicit properties of such polynomial sets, generalizing Laguerre and Legendre polynomials. The nature of the dimension 1 subspace completing such sets is described. A numerical example illustrates the use of such polynomials.

I. INTRODUCTION

Generalizations Γ n of Hermite polynomials H n were recently [1] proposed to describe, for instance, density perturbations constrained by a condition of matter conservation. Because of the constraint, such polynomials cannot form a complete set, but span a subspace well suited to specific applications. In particular, the polynomials Γ n used in [1] were motivated by the consideration in nuclear physics of the Hohenberg-Kohn functional [2] and similar functionals along the Thomas-Fermi method [3,4]. Indeed, in such approaches, the ground state of a quantum system is shown to be a functional of its density ρ(r), and there is a special connection between ρ(r) and the mean field u(r) driving the system. It was thus convenient to expand variations of ρ in a basis {w m (r)} of particle number conserving components, δρ(r) = m δρ m w m (r), with the term-by-term constraint, ∀m, dr w m (r) = 0. This spares, in the formalism, the often cumbersome use of a Lagrange multiplier. Simultaneously, it was convenient to expand variations of u in a basis orthogonal to the flat potential, because, trivially, a flat δu, as just a change in energy reference, cannot influence the density. The same basis can thus be used for δu(r) = n δu n w n (r), since the very same condition, dr w n (r) = 0, induces orthogonality to a constant δu. Because of the nuclear physics context of [1], harmonic oscillators shell models were considered and the basis contained a Gaussian factor, e -1 2 r 2 . The same functional approaches [2][3][4] are also of a general use in atomic and molecular physics, where Gaussian weights would be clumsy and radial properties are best fitted with simple exponential weights [START_REF] Messiah | Mécanique Quantique[END_REF]. Furthermore, in [1], the discussion was restricted to one dimensional problems. In the present note, we want to include two and three dimensional situations. We shall thus use weights of the form e -1 2 r , with 0 ≤ r < ∞, but integrals will carry a factor r ν , with ν a positive exponent, suitable for dimension d. This will lead to generalizations of Laguerre polynomials.

This note is also concerned with compact domains, of the form 0 ≤ r ≤ 1 for instance. This might correspond for instance to expansions of density fluctuations in cylindrical vessels used for chemical processes, where mass conservation is also in order, or maybe in centrifuges. Radial integrals with factors r and r 2 in both the constraint and orthogonalization conditions will lead to generalizations of Legendre polynomials.

For any positive weight µ(r), and any dimension d, a constraint of vanishing average, dr r ν µ(r) Γ n (r) = 0, is incompatible with a polynomial Γ of order n = 0. Therefore, in the following, the order hierarchy for the constrained polynomials runs from n = 1 to ∞, while that for the traditional polynomials runs from 0 to ∞. We study in some generality the "Laguerre" case in Section II. In turn, the "Legendre" case is the subject of Section III. A brief Section IV discusses possible applications to the study of density fluctuations in centrifuges. Section V answers a question which was omitted in [1], that of the nature of the projector onto the subspace spanned by the constrained states and the nature of the codimension of this subspace. A numerical application is provided in Section VI. A discussion and conclusion make Section VII.

II. MODIFICATION OF LAGUERRE POLYNOMIALS BY A CONSTRAINT OF ZERO AVERAGE

In this Section we consider basis states carrying a weight e -1 2 r , in the form

w n (r) = e -1 2 r G d n (r), where G d n is a polynomial. It is clear that G d
0 cannot be a finite, non vanishing constant if the constraint, ∞ 0 dr r d-1 e -1 2 r G d 0 (r) = 0, must be implemented. Hence set integer labels m ≥ 1 and n ≥ 1 and define polynomials G d n by the conditions,

∞ 0 dr r d-1 e -r G d m (r) G d n (r) = g d n δ mn , ∞ 0 dr r d-1 e -1 2 r G d n (r) = 0 , (1)
1 where δ mn is the usual Kronecker symbol and the positive numbers g d n are normalizations, to be defined later. It is elementary to generate such polynomials numerically, in two steps by brute force, namely i) first create "trivial seeds" of the form, s d n (r) = r n -r n d , where the subtraction of the average, r n d = 2 n (d -1 + n)!/(d -1)!, ensures that each trivial seed fulfills the constraint, then ii) orthogonalize such seeds by a Gram-Schmidt algorithm. The first polynomials read,

G 1 1 = r -2, G 1 2 = r 2 -5r + 2, G 1 3 = r 3 -10r 2 + 20r -8, G 1 4 = r 4 -17r 3 + 78r 2 -108r + 24, (2a)
G 2 1 = r -4, G 2 2 = r 2 -8r + 8, G 2 3 = r 3 -14r 2 + 44r -32, G 2 4 = r 4 -22r 3 + 138r 2 -288r + 144, (2b) G 3 1 = r -6, G 3 2 = r 2 -11r + 18, G 3 3 = r 3 -18r 2 + 78r -84, G 3 4 = r 4 -27r 3 + 216r 2 -606r + 468. ( 2c
)
All these are defined to be "monic", namely the coefficient of r n is always 1. For an illustration we show in Figure 1 the new polynomials G Rather using the Gram-Schmidt method, we find it easier, and more elegant, to generate the polynomials G d n , starting from the initial table, Eqs. (2a,2b,2c), by means of the following recursion formula,

G d n (r) = (r -d) G d n-1 (r) -2 r G d ′ n-1 (r) + (n + d -1) (n -2) G d n-2 (r), (3)
where the prime denotes the derivative with respect to r. Its simple structure can be proven analytically as follows: i) Let us first create some kind of a "less trivial seed" at order n, assuming the polynomial G d n-1 is known. For this, try r G d n-1 . By partial integration, we see that,

∞ 0 dr r d-1 e -1 2 r r G d n-1 (r) = 2 ∞ 0 dr e -1 2 r r d G d n-1 (r) ′ , (4)
where again a prime means derivation with respect to r. Thus ii) It turns out that such coefficients vanish if m < n -2. Indeed, an integration of σ d n against G d m , weighted by r d-1 e -r , gives, by partial integration of the

σ d n ≡ r G d n-1 -2 r G d ′ n-1 -2 d G d n-
G d ′ n-1 term, ∞ 0 dr e -r r d-1 σ d n (r) G d m (r) ≡ ∞ 0 dr r d-1 e -r (r -2d) G d n-1 (r) -2r G d ′ n-1 (r) G d m (r) = ∞ 0 dr e -r r d-1 G d n-1 (r) (r -2d) G d m (r) + 2 ∞ 0 dr G d n-1 (r) e -r r d G d m (r) ′ = ∞ 0 dr e -r r d-1 G d n-1 (r) -σ d m+1 (r) -2d G d m (r) . (5)
In the bracket [ ] in the last right-hand side of Eq. ( 5) the seed σ d m+1 has order m + 1 and, by definition, G d m is of order m. By definition also, G d n-1 , of order n -1, is orthogonal to all those polynomials of lower order, that are compatible with the constraint. This integral, Eq. ( 5), thus vanishes as long as m + 1 < n -1. It can be concluded that the difference, σ d n -G d n , contains only two contributions, namely those from G d n-2 and G d n-1 . Explicit forms for their coefficients are obtained by elementary manipulations, leading to Eq. ( 3). Elementary manipulations also give,

2 r G d ′′ n -(r -2d) G d ′ n + n G d n = (n -1) (n + d) G d n-1 . (6)
Here, in the same way as a prime means first derivative with respect to r, we used double primes for second derivatives.

Finally the normalization of the polynomials is obtained easily as,

g d n ≡ ∞ 0 dr e -r r d-1 [G d n (r)] 2 = (n -1)! (n + d)! . (7)

III. MODIFICATION OF LEGENDRE POLYNOMIALS BY A CONSTRAINT OF ZERO AVERAGE

Legendre polynomials, and their associates and generalizations (Gegenbauer, Chebyshev, Jacobi) are defined with respect to the [-1, 1] segment. Exceptionally in the literature, one finds shifted Legendre polynomials, adjusted to the [0, 1] segment. We are here interested in applications to radial densities in cylinders, or the small circles of toruses, or spheres. Hence we shall use 0 ≤ r ≤ 1. It is clear that the case, d = 1, does not make an original problem, since Legendre polynomials, whether translated and/or scaled or not, already average to 0 as soon as their order n is ≥ 1. We keep the case, d = 1, for the sake only of completeness and in this Section we consider d = 1, 2, 3, with a geometry factor r d-1 . The weight is µ(r) = 1, hence our states are described by just a polynomial G d n of order n. It is again obvious that G d 0 cannot be a non vanishing constant if the constraint,

1 0 dr r d-1 G d 0 (r) = 0, is implemented. Hence set m ≥ 1, n ≥ 1, and define polynomials G d n from conditions, 1 0 dr r d-1 G d m (r) G d n (r) = γ d n δ mn , 1 0 dr r d-1 G d n (r) = 0 , (8)
where the normalizations γ d n are again to be defined later. It is obvious that the shifted (and shrunk) Legendre polynomials L n (2r -1), n ≥ 1, satisfy both constraint and orthogonality relations for d = 1, because they are orthogonal to any constant polynomial, of order 0. The polynomials G 1 n = L n (2r -1) thus make nothing new. We turn therefore to d = 2 and d = 3, with a brute force construction as in the previous Section. But the defining conditions, Eqs. (8), show a difference with Eqs. (1): both orthogonality and constraint conditions now carry the same weight, namely µ 2 = µ, while in the previous case, Eqs. (1), there were different weights, because of the exponentials e -r and e -1 2 r . A similar difference between µ 2 and µ happened in the "Hermite" case, naturally. Hence now, in this Legendre case, we can Gram-Schmidt orthogonalize even more trivial seeds r n , without subtractions, and accept those orthogonal polynomials with order m ≥ 1. The table of first results reads,

G 1 1 = 2r -1, G 1 2 = 6r 2 -6r + 1, G 1 3 = (2r -1)(10r 2 -10r + 1), G 1 4 = 70r 4 -140r 3 + 90r 2 -20r + 1, (9a) G 2 1 = 3r -2, G 2 2 = 10r 2 -12r + 3, G 2 3 = 35r 3 -60r 2 + 30r -4, G 2 4 = 126r 4 -280r 3 + 210r 2 -60r + 5, (9b) G 3 1 = 4r -3, G 3 2 = 15r 2 -20r + 6, G 3 3 = 56r 3 -105r 2 + 60r -10, G 3 4 = 210r 4 -504r 3 + 420r 2 -140r + 15. (9c)
Easy, but slightly tedious manipulations validate the following recursion relations,

n G 1 n = (2n -1) (2r -1) G 1 n-1 -(n -1) G 1 n-2 , (10a)
(n + 1) (2n -1) G 2 n = 2 [ (4n 2 -1)r -2n 2 ] G 2 n-1 -(n -1) (2n + 1) G 2 n-2 , (10b) n 2 (n + 2) G 3 n = (2n + 1) [ 2n(n + 1)r -(n 2 + n + 1) ] G 3 n-1 -(n -1) (n + 1) 2 G 3 n-2 . ( 10c
)
and the differential equation,

r (r -1) G d ′′ n + [ (d + 1) r -d ] G d ′ n -n (n + d) G d n = 0. (11)
Finally the normalization of the polynomials reads,

γ d n ≡ ∞ 0 dr r d-1 [G d n (r)] 2 = 1/(2n + d) . (12)
We show in Figure 2 the plots of G d n for n = 1, 2 and d = 1, 2, 3.

0.2 0.4 0.6 0.8 1 r -2 2 4 6 CstLeg FIG. 2. Modified Legendre polynomials G 1 1 , G 1 2 (full lines), G 2 1 , G 2 2 (long dashes), G 3 1 , G 3 2 (dashes).

IV. POLYNOMIALS FOR CENTRIFUGES

The case of centrifuges is worth a short comment. As soon as the matter under centrifugation is compressible, the density becomes much larger at the outer edge, r = 1, than at the rotation axis, r = 0. Let h be the height of the centrifuge. Assume, for the sake for the argument, that one studies fluctuations about a reference density of the form, ρ(r) = ρ c e Kr 2 , where the parameter K contains all informations about the angular velocity, compressibility, etc. of the process. The factor, ρ c = M h 1 0 dr r ρ(r)

-1 = M h -1 2 K e K -1 -1
, ensures the conservation of the mass M included in the vessel. If a cause for fluctuations of ρ is an instability of K, the first order for density change is,

∂ρ ∂K (r) = 2 K r 2 e K -1 + e K -Ke K -1 [e K -1] 2 e Kr 2 , 1 0 dr r ∂ρ ∂K (r) = 0, ( 13
)
namely a polynomial of order 2 multiplied by e Kr 2 . Higher derivatives with respect to K will generate similar, even order polynomials, with the same property, 1 0 dr r ∂ n ρ/∂K n (r) = 0. An orthogonalization, under a metric ∝ e 2Kr 2 , might be useful. This new set of polynomials will depend on K, however, since r is already scaled to a radius 1 for the cylinder and thus K cannot be scaled away. Because of this K dependence we do not elaborate further on this issue. For a large list of ad hoc polynomials and integration weights, see [START_REF] Boujot | Thèse de 3e cycle, Calcul de moments appliqués à l'évaluation numérique d'intégrales[END_REF].

V. PROJECTOR ON THE CONSTRAINED SUBSPACE

For the sake of the discussion and short notations, set first d = 1, µ(r) = e -1 2 r , and temporarily include normalization to unity factors into both the Laguerre polynomials L n and the constrained G 1 n . This summarizes as,

∞ 0 dr [µ(r)] 2 L m (r) L n (r) = δ mn , ∞ 0 dr [µ(r)] 2 G 1 m (r) G 1 n (r) = δ mn , ∞ 0 dr µ(r) G 1 n (r) = 0, (14)
|σ N = N m=0 z m 2 -1 2 N n=0 z n |z n , z n = ∞ 0 dr r|z n , (15)
defines the proper "subtractor" operator |σ N σ N |. Indeed, from

( Q N -|σ N σ N | ) |τ = N n=0 τ n |z n - N m=0 z m 2 -1 N n=0 z n |z n N p=0 z p τ p , (16) one obtains
∞ 0 dr r| (Q N -|σ N σ N |) |τ = N n=0 τ n z n - N m=0 z m 2 -1 N n=0 z n z n N p=0 z p τ p = 0. (17)
Hence Q N -|σ N σ N | is the projector P N . Incidentally, the Laguerre result for σ N is very simple, because z m = 2, ∀m. But the ansatz for σ N , Eq.( 15), generalizes to all cases. For instance, with Hermite polynomials, odd orders already satisfy the constraint when integrated from -∞ to ∞, naturally, and thus do not contribute to σ N . Even orders contribute, and it is easy to verify, upon integrating from -∞ to ∞ again, that z 2p 2 = π 1 2 2 1-p (2p -1)!!/p!. It may be pointed out that the condition, dr µ(r) f (r) = 0, for functions f orthogonalized, like our polynomials, by a metric [µ(r)] 2 , might be interpreted as an orthogonality condition, dr f (r) [µ(r)] 2 g(r) = 0, with g(r) = [µ(r)] -1 . This makes g a candidate for the subtractor form factor σ. This is of some interest for the centrifuge case, where a state function such as, for instance, e -Kr 2 , remains finite when 0 ≤ r ≤ 1. But there is little need to stress that, when the support of µ extends to ∞, then µ -1 does not belong to the Hilbert space and cannot be used for σ.

More interesting is the limiting process, N → ∞, as illustrated by Figures 2345. Figs. 2 and3 show the shapes, in terms of r, of 2|P N |r and 10|P N |r , respectively, when the projectors are made of the modified Laguerre polynomials G 1 n . The build up of an approximate δ-function when N increases from N = 50 (short dashes) to N = 100 (long ones) and N = 150 (full lines) is transparent, although the convergence is faster when peaks are closer to the origin, compare Figs. 2 and3. The slower convergence in Fig. 3 is due to the cut-off imposed by exponential weights as long as N is finite. Given N, there is a "box effect", the range of the box being of order ∼ N. A similar build up is observed for our other families of constrained polynomials, with slightly different details of minor importance such as, for instance, a box range ∼ √ N for the Hermite case. The box effect is even more transparent in Figs. 4 and5, which show the shapes of subtractors 10|σ N σ N |r and 0|σ N σ N |r deduced from constrained polynomials of the Laguerre (Fig. 4) and Hermite (Fig. 5) type, respectively. (For graphical convenience, the polynomials Γ 1 n and H n used for the Hermite case, Fig. 5, are tuned to a weight e -r 2 rather than e -1 2 r 2 , but this detail is not critical.) It seems safe to predict that, given an effective length Λ(N ) for the box, the wiggles of the subtractor will smooth out when N → ∞ and that only a background ∼ -1/Λ(N ) will then remain.

VI. ILLUSTRATIVE EXAMPLE: TRAJECTORIES IN DENSITY SPACE

We return here to the toy model discussed in [1] and the corresponding, modified Hermite polynomials. The model consists of Z non interacting, spinless fermions, driven by a one dimensional harmonic oscillator H 0 = 1 2 (-d 2 /dr 2 +r 2 ). The ground state density from the Z lowest orbitals reads, ρ(r) = Z i=1 [ψ i (r)] 2 . Let i = 1, .., Z and I = Z + 1, ... ∞ label "hole" and "particle" orbitals, respectively. Add a perturbation δu(r) to the initial potential r 2 /2. The first order variation of the density is,

δρ(r) = 2 iI ψ i (r)ψ I (r) I|δu|i i -I . ( 18
)
If we expand δu and δρ in that basis {w n } provided by the new polynomials, the formula, Eq. ( 18), becomes,

δρ m = 2 iI n D m iI 1 i -I D n iI δu n , D n iI ≡ dr w n (r) ψ i (r)ψ I (r), (19)
where D denotes both a particle-hole matrix element of a potential perturbation and the projection of a particlehole product of orbitals upon the basis {w n }. In [1] we briefly studied the eigenvalues and eigenvectors of this symmetric matrix, F = D (E 0 -H 0 ) -1 D, where (E 0 -H 0 ) -1 is a short notation to account for the denominators and the particle-hole summation, and the tilde indicates transposition. It is clear that the invertible F represents the functional derivative δρ m /δu n and is suited for infinitesimal perturbations. We shall now take advantage of the representation provided by {w n } to study finite trajectories ρ(u).

For this, we consider a variable Hamiltonian, H m (λ) = H 0 + λ w m (r), made of the initial harmonic oscillator, but with a finite perturbation ∆u along one "mode" w m . It is trivial to diagonalize H m (λ) with an excellent numerical accuracy and thus obtain, given Z, the ground state density ρ(r, λ). Then it is trivial to expand the finite variation, ∆ρ = ρ(r, λ) -ρ(r, 0), in the basis {w n }. This defines coordinates ∆ρ n (λ; m) for trajectories, parametrized by the intensity of the chosen mode m for ∆u. In Figures 7 and8 we show, with Z = 4, results from H 4 = H 0 + λ 4 2 (2π) -1 4 15 -1 2 (8r 4 -14r 2 + 1) e -r 2 and H 6 = H 0 +λ 6 (2π) -1 4 105 -1 2 (32r 6 -128r 4 +94r 2 -11) e -r 2 , respectively. The case, H 2 = H 0 +λ 2 2 (2/π)

1 4 3 -1 2 (2r 2 -1) e -r 2 ,
makes almost a harmonic oscillator and is probably of academic interest only; anyhow we verified that its confirms the results with H 4 and H 6 . We use a basis {w n } containing a factor e -r 2 rather than e -1 2 r 2 to better match the same factor e -r 2 created by products of harmonic oscillator orbitals in the calculation of matrix elements z p |∆u|z q , but this technicality is not important for the physics.

The main result to be observed seems to be the lack of "collectivity" for such modes and for such elementary Hamiltonians. Indeed, for λ 4 = 2, the first five coordinates of ∆ρ read {0.016, -0.267, -0.055, 0.023, 0.018}, with a strong dominance of ∆ρ 4 , while for λ 6 = 2, these read {-0.013, -0.041, -0.376, 0.008, 0.040}, with a strong dominance of ∆ρ 6 . To clarify Figs. 7 and8, we had indeed to blow up each ∆ρ n by a factor 2 |n-m| , where m is the index of the driver mode in potential space. Other modes than m = 4 and m = 6 show the same property: in the density space, a trajectory driven by ∆u = λw m stays close to the same w m axis in that density space, although curvatures effects, while somewhat modest, are not absent. Such non linearity, slight curvatures are seen in Figs. 78, and also in Figure 9, where the three ∆ρ 4 , ∆ρ 6 , ∆ρ 8 sets of data shown by Fig. 7 are converted into a parametric plot for a trajectory. For graphical purposes again, ∆ρ 6 and ∆ρ 8 are blown up 4 times to create Fig. 9. It can be concluded, temporarily, that the "flexibility" matrix F is not too far from being diagonal in the {w n } basis, or in other words, that the w n modes indicate an approximately natural hierarchy in both the potential and the density spaces.

A subsidiary question pops up: that of the positivity of ρ. Indeed, while the space of potentials is basically a linear space, with arbitrary signs for u(r) when the position r changes, densities ρ(r) must remain positive for every r. This creates severe constraints for any linear parametrization of ∆ρ in terms of the basis {w n }. In our toy model, it turns out that ρ(r, 0) = π -1 2 (8r 6 -12r 4 + 18r + 9) e -r 2 /6. Hence, if we truncate ∆ρ to have two components only, w 2 and w 4 for instance, then ρ is the product of e -r 2 and a polynomial P(r), 6 π 1 2 P(r) = 8r 6 -12r 4 + 18r 2 + 9 + ∆ρ 2 12 (2π)

1 4 3 -1 2 (2r 2 -1) + ∆ρ 4 12 (π/2) 1 4 15 -1 2 (8r 4 -14r 2 + 1). ( 20
)
Rescale out inessential factors, for a simpler polynomial, P = 8r 6 -12r 4 +18r 2 +9+∆R 2 (2r 2 -1)+∆R 4 (8r 4 -14r 2 +1). Eliminate r between P and d P/dr. The resultant R(∆R 2 , ∆R 4 ), when it vanishes, gives the border of the convex domain of parameters ∆R 2 , ∆R 4 where P remains positive definite. This domain contains the origin, because of ρ(r, 0). The precise form of R is a little cumbersome and does not need to be published here. But the resulting border is shown in Figure 10. Generalizations to more ∆ρ parameters are obvious, with more cumbersome resultants R.

VII. DISCUSSION AND CONCLUSION

The subject of orthogonal polynomials has been so treated and overtreated that any claim to novelty must contain much more than a change of the integration measure. We took therefore a different approach, motivated by a law of physics and/or chemistry, matter conservation. This means a constraint of a vanishing average for the states described by weighted polynomials.

For a support [0, ∞[ and a simple exponential weight such as e -1 2 r , a non trivial generalization of Laguerre polynomials occurs. This extends the generalization of Hermite polynomials described in [1] with the support ] -∞, ∞[ and Gaussian weights such as e -1 2 r 2 . We also took care of cylindrical and spherical geometries, by replacing dr with dr r and dr r 2 , respectively. The new sets of constrained polynomials are clearly sensitive to the geometry.

For finite supports such as [0, 1] and constant weights, the constraint is already satisfied by the usual brand of orthogonal polynomials as soon as their order is ≥ 1. In that sense, we did not find significantly original generalizations of Legendre polynomials, although we generated polynomials fitted to the cylindrical and spherical geometries. The cause of the failure is transparent: when the weight µ(r) is a constant, there is no difference between the orthogonality metric µ 2 and the constraint weight µ.

For each set of new polynomials we found a recursion relation and a differential equation. There seems to be a systematic property for those cases where the constraint generates truly original polynomials, namely when µ 2 = µ. In such cases, recursion and differentiation seem to be necessarily entangled. This does not happen for traditional orthogonal polynomials, indeed, and this "entanglement" may deserve some future attention.

Constrained polynomials expressing matter conservation in centrifuges do make an original set if the fluid under centrifugation is compressible; a non constant reference weight µ is indeed in order there. But the set depends on the precise form of µ via potentially many physical parameters. We found it difficult to design, through scaling, a sufficiently "universal" set. "Centrifuge polynomials" will have to be calculated specifically for each practical situation.

For those new polynomials generalizing the Hermite and Laguerre ones, we found a description of the subspace accounting for their defect of completeness. A codimension 1 is the consequence of the constraint, expressed at first by the obvious lack of a polynomial of order n = 0.

Finally the use of such polynomials was illustrated by a toy model for the Hohenberg-Kohn functional. A slightly surprising result was found: our polynomials, those of low order at least, define potential perturbations which are reflected by density perturbations having almost the same shapes. This occurs despite the delocalization created by the kinetic energy operator, hints at short ranges in effective interactions and validates the localization spirit of the Thomas-Fermi method. Whether such hints are good when the full zoology of the density functional is investigated is, obviously, an open question; for a review of the richness of the functional, we refer to [START_REF] Perdew | A Primer in Density Functional Theory[END_REF]. If long range forces are active, a significant amount of delocalization between the "potential cause" and the "density effect" is not excluded. It would be interesting indeed to discover collective degrees of freedom in this connection between potential and density. In any case, our main conclusion may be that the new polynomials provide, for the context of matter conservation, a discrete and full set of modes and coordinates, hence a systematic and constructive representation of phenomena.

It is a pleasure to thank Y. Abe, J.-P. Boujot, B. Eynard, C. Normand, R. Peschanski and A. Weiguny for stimulating discussions.

  Then the kets and bras defined by r|w n = w n |r = w n (r) = µ(r) G 1 n (r) and r|z n = z n |r = z n (r) = µ(r) L n (r) provide two "truncation" projectors, P N = N n=1 |w n w n | and Q N = N n=0 |z n z n |, available for subspaces where polynomial orders do not exceed N. Their respective ranks N and N + 1, and the embedding and commutation relation, [P N , Q N ] = P N , are obvious. Obvious also is the limit, lim N →∞ Q N = 1. The role of the rank one |σ N σ N | difference P N -Q N is to subtract from any test state, |τ = N n=0 τ n |z n , that part which violates the condition of vanishing average. We shall show that the elementary ansatz,

FIG. 3 .

 3 FIG. 3. Shapes of projectors made of polynomials G 1n . Full line, 2|P150|r , long dashes, 2|P100|r , short dashes 2|P50|r .

FIG. 4 .FIG. 5 .

 45 FIG. 4. Shapes of projectors made of polynomials G 1n . Full line, 10|P150|r , long dashes, 10|P100|r , short dashes 10|P50|r .

FIG. 6 .

 6 FIG. 6. Subtractors made of Γn. Shapes centered at r = 0. Stronger wiggles, shorter cut-off, dashed line, N = 50. Weaker wiggles, larger cut-off, full line, N = 100.

FIG. 7 .

 7 FIG. 7. Coordinates of the perturbation density ∆ρ created by a perturbing potential ∆u = λ4 w4. Full line: 2 ∆ρ2. Long dashes: ∆ρ4. Moderate dashes: 2 ∆ρ6. Short dashes: 4 ∆ρ8. Very short dashes: 8 ∆ρ10.

FIG. 8 .FIG

 8 FIG.8. Same as Fig.7, but now ∆u = λ6 w6. Full line: 4 ∆ρ2. Long dashes: 2 ∆ρ4. Moderate ones: ∆ρ6. Short ones: 2 ∆ρ8. Very short dashes: 4 ∆ρ10.

  FIG. 10. Domain of values of ∆R2 and ∆R4 acceptable for the positivity of the density of the toy model. The domain sits inside the full line curve and left of the straight line. It contains the origin.


B G Giraud

R Peschanski

ON POSITIVE FUNCTIONS WITH POSITIVE FOURIER TRANSFORMS

Keywords: numbers: 02.30.Gp, 02.30.Mv, 02.30.Nw, 12.38.-t

Using the basis of Hermite-Fourier functions (i.e. the quantum oscillator eigenstates) and the Sturm theorem, we derive constraints for a function and its Fourier transform to be both real and positive. We propose a constructive method based on the algebra of Hermite polynomials. Applications are extended to the 2-dimensional case (i.e. Fourier-Bessel transforms and the algebra of Laguerre polynomials) and to adding constraints on derivatives, such as monotonicity or convexity.

Introduction

Positivity conditions for the Fourier transform of a function occur in various domains of physics. One often asks:

• What are the constraints for a real function ψ(r) ensuring that its Fourier transform

ϕ(s) = 1 √ 2π +∞ -∞
dr e isr ψ(r)

be real and positive?

• Conversely, what are the properties of ϕ if ψ is real and positive?

• Finally, what are the constraints on Fourier partners such that both ψ and ϕ be real and positive? † giraud@dsm-mail.saclay.cea.fr ‡ pesch@dsm-mail.saclay.cea.fr (331) Physicists often need to work with concrete constructions. The practical construction of a basis of functions satisfying the above mentioned positivity properties remains, up to our knowledge, an open problem. Such questions are quite relevant in Physics. As practical examples, let us quote two wellknown cases. A Fourier transform relates [START_REF] Glauber | High-Enery Physics and Nuclear Structure[END_REF] two quantities, namely the cross section and the profile of a nucleus, which both ought to be positive. In particle physics, a 2-d Fourier-Bessel transform relates the color dipole distribution in transverse position space (derived from Quantum Chromodynamics) and the transverse momentum distribution of gluons probed during a deep-inelastic collision [START_REF] Kovchegov | [END_REF]. Such questions occur also in probability calculus, for the relation between probability distributions and characteristic functions [START_REF] Lévy | Fonctions caractéristiques positives [Positive characteristic functions[END_REF], in crystallography and in condensed matter physics, e.g. for the interpretation of patterns, etc. . . . . The problem is simplified if related to another one, that concerning the functions which are invariant [START_REF] Titchmarsh | Introduction the Theory of Fourier Integrals[END_REF] up to a phase factor 1 by Fourier transforms. Indeed, the most familiar examples of positive self-dual functions or distributions, thus trivially verifying the double positivity condition, which are at most scaled under Fourier transformation (FT), are the Gaussian and the Dirac comb.

Many special cases can be found, where positivity is conserved, such as, for instance, the continuous family of functions exp(-r ν ), where 0 < ν ≤ 2. Various sufficient conditions for positivity can be found in the literature, such as the convexity of ψ [START_REF] Lafforgue | private communication; for a general theorem[END_REF] but, up to our knowledge, no general constructive method has been presented.

The present note attempts to give general positivity criteria, in a constructive way, by taking advantage of a representation under which the FT is essentially "transparent". Our method combines the advantages of selfduality properties with those allowed by an algebra of polynomials, where positivity means absence of real roots, hence reasonably simple conditions for the polynomial coefficients. For this sake, in the 1-d case, we select a basis made of convenient eigenstates of the FT, the Hermite-Fourier functions, i.e. the harmonic oscillator eigenstates. The method extends to the 2-d case, or Fourier-Bessel transform, by replacing Hermite by Laguerre polynomials.

There are general theorems about the characterization of Fourier transforms of positive functions [START_REF] Scharnhorst | For a comprehensive list of references[END_REF]. Let us quote in the first place the Bochner theorem and its generalizations [START_REF] Gel'fand | Generalized Functions[END_REF] which state that the Fourier transform of a positive function is positive-definite. But positive definiteness in the sense of such theorems does not imply plain positivity 2 . Hence our problem 1 They are called "self-Fourier" in [START_REF] Caola | [END_REF] if the phase factor is 1 and "generalized self-Fourier" [6] (or "dual") for other phases. 2 Positive definiteness means that for any real numbers x1, . . . x k and complex numbers ξ1, . . . ξ k , one has

P k,j ϕ(x k -xj) ξjξk ≥ 0.
actually could be rephrased [START_REF] Teugels | [END_REF] as "build positive-definite functions that are positive". Our formalism is the subject of Section 2. Numerical, illustrative examples are given in Section 3. Then Section 4 extends our algorithms to the 2-d problem. Hermite polynomials are replaced by Laguerre polynomials, but the algebra remains essentially the same. A brief discussion, conclusion and outlook are offered in Section 5.

Basic formalism

Consider the harmonic oscillator Hamiltonian, 1 2 (p 2 + r 2 ), and its eigenwavefunctions

u n (r) = π -1/4 e -1/2r 2 H n (r) . (2)
Here, we set H n to be a square normalized Hermite polynomial, (-) n e r 2 d n /dr n e -r 2 / √ 2 n n!, with a positive coefficient for its highest power term. For the sake of clarity, we list the first polynomials as,

H 0 = 1, H 1 = √ 2 r, H 2 = (2r 2 -1)/ √ 2, H 3 = (2r 3 -3r)/ √ 3
and their recursion relation

a n+1 H n+1 = 2r a n H n -2n a n-1 H n-1 , (3)
where

a n = √ 2 n n! .
It is known that the FT of such states brings only a phase

1 √ 2π ∞ -∞ dr e isr u n (r) = i n u n (s) , (4)
and thus such states give generalized self-dual functions with phase i n . Let one expand ψ in the oscillator basis, ψ(r) = N n=0 ψ n u n (r), with a truncation at some degree N. Then all odd components ψ 2p+1 must vanish if ϕ must be real, and the even rest splits, under FT, into an invariant part and a part with its sign reversed, namely

ϕ(s) = π -1/4 e -1/2 s 2 [P + (s) -P -(s)] , P + (s) = [N/4] p=0 ψ 4p H 4p (s) , P -(s) = [(N -2)/4] p=0 ψ 4p+2 H 4p+2 (s) , (5)
where the usual symbols [N/4] and [(N -2)/4] mean, respectively, the entire parts of N/4 and (N -2)/4.

Notice that, when all components ψ n vanish except ψ 0 , then both ψ and ϕ are positive, because H 0 = 1. Hence, one may, starting from this special point in the functional space of functions, investigate those domains of parameters ψ n where the polynomials P = P + + P -and Q = P + -P - have no real root. Notice that only even powers of r and s are involved. It will therefore be convenient to use auxiliary variables ρ = r 2 and σ = s 2 , and the domain of interest for the parameters ψ n will correspond to the absence of real positive roots for both ρ and σ.

The second ingredient of our approach is the well-known Sturm theorem [START_REF] Sturm | Mémoire sur la résolution des équations numériques[END_REF] which gives a way to characterize and localize the real roots of any given polynomial. The Sturm criterion can be expressed as follows:

"Given a polynomial P(x), its Sturm sequence S(x) ≡ {S 1 , S 2 , . . . S m , . . . S j≤N } is the set of polynomials

S 1 = P, S 2 = dP dx , . . . S m = -S m-2 + S m-2 S m-1 S m-1 . . . , (6)
where [ ] designates the polynomial quotient3 . To know the number of distinct roots between x = a and x = b, count the number N (a) of sign changes in S(a) and, similarly, count N (b). Then the number of roots is

|N (b) -N (a)| ."
The domain borders where the root number, |N (+∞) -N (0)| , changes have to do with cancellations of the resultant R between P and dP/dx. The cancellation of R corresponds to collisions between conjugate complex roots becoming real roots and conversely. Because of the demanded positivity of ρ and σ, the borders have also to do with sign changes of P(0) and Q(0), meaning real roots ρ and σ going through 0. All such technicalities are taken care of by the Sturm criterion, which furthermore allows the labeling of each domain by its precise number of real roots. This will be implemented here in an explicit way, analytically as much as possible, then numerically and graphically, for a few cases of a general illustrative value. For this, we will plot the shapes of domains labeled by the values of the Sturm criterion. Hence the combination of the self-dual properties of the quantum oscillator basis and of the Sturm theorem allows a constructive method for a systematic investigation of positivity conditions for a 1-d Fourier transform. This basis has the potential to represent any function in the Hilbert space L 2 , but the present study concerns a finite set of components. An extension to infinite series of components may deserve other tools.

Positivity domains

In this section we will apply our method to the case of a basis formed with 3 or 4 Hermite-Fourier functions. In Sec. 3.1 we consider the basis ψ 0 , ψ 4 , ψ 8 , then in 3.2 the basis ψ 0 , ψ 4 , ψ 8 , ψ 12 , which both lie in the subspace with eigenvalue 1, and in 3.3 the basis ψ 0 , ψ 2 , ψ 4 where ψ 2 is in the subspace with eigenvalue -1, furthermore, in 3.4 the influence of an additional constraint motivated by Physics, that of monotony for ϕ, and finally in 3.5 a comparison with the convexity constraint for ψ. Many other illustrations are possible, but these, Sec. 3.1-3.5, demonstrate the flexibility of our approach.

Mixture of three polynomials in the subspace with eigenvalue 1

Here we assume that ψ has only components ψ 0 , ψ 4 and ψ 8 , hence P reduces to P + in (5) and we can study

P = ψ 0 + ψ 4 2 √ 6 (4ρ 2 -12ρ+3)+ ψ 8 24 √ 70 (16ρ 4 -224ρ 3 +840ρ 2 -840ρ+105), (7)
where ρ = r 2 . One maintains ψ 8 > 0 and one reads P + (0) = ψ 0 + 3ψ

Because of the free scaling of P we normalize the polynomial so that its coefficients lie on a half sphere of unit radius

ψ 0 = cos α cos β , ψ 4 = sin α cos β , ψ 8 = sin β, -π < α ≤ π, 0 < β ≤ π/2 . ( 9
)
We thus show in figure 1 the domains where the number of roots increases from 0 to 4. The no root domain is the white triangle above α = β = 0 and slightly right from this point. In this domain, ψ (≡ ϕ) is self-Fourier and positive.

Mixture of four polynomials in the subspace with eigenvalue 1

Now we add to ψ a component ψ 12 , hence

P = ψ 0 + ψ 4 2 √ 6 (4ρ 2 -12ρ + 3) + ψ 8 24 √ 70 (16ρ 4 -224ρ 3 + 840ρ 2 -840ρ + 105) + ψ 12 1440 √ 231 (64ρ 6 -2112ρ 5 + 23760ρ 4 -110880ρ 3 + 207900ρ 2 -124740ρ + 10395). (10)
While borders corresponding to P + (0) = 0 obtain easily, the resultant R to be considered for other borders is unwieldy and is skipped here. Taking advantage of scaling we set:

ψ 0 = cos α cos β cos γ , ψ 4 = sin α cos β cos γ , ψ 8 = sin β cos γ , ψ 12 = sin γ , |α| ≤ π 2 , |β| ≤ π , 0 ≤ γ ≤ π 2 . ( 11
)
This choice of spherical coordinates was designed to ensure the positivity of ψ 12 , obviously, but also a dominance of ψ 0 near α = β = γ = 0. The dominance is clearly useful for the positivity of P. Then this S 3 sphere can be explored by various cuts according to fixed values of γ. The results are shown in figures 2-5, with γ = π/10 6 , π/15, 2π/15, π/6, π/5, 7π/30, 4π/15, 3π/10 , respectively. The color code for the number of roots is: 0 root, red; 1, yellow; 2, yellowish green; 3, bluish green; 4, blue; 5, dark purple; 6, pink (in the uncolored edition, the no root domain, if it exists, is that dark, small or tiny triangle slightly right of the map center). The red domain shrinks at first very slowly when γ increases, then faster when γ ≃ π/6. Beyond such an order of magnitude for γ, there is no red domain and the map becomes invaded by bigger and bigger pink patches, representing the dominance of the 6 positive, real roots of H 12 .

3.3. Two polynomials from subspace "1" mixed with one polynomial from subspace "-1"

If we consider a mixture of H 0 , H 2 and H 4 , the FT connects the two polynomials

P(ρ) = ψ 0 + ψ 2 2ρ -1 √ 2 + ψ 4 4ρ 2 -12ρ + 3 2 √ 6 , Q(σ) = ψ 0 -ψ 2 2σ -1 √ 2 + ψ 4 4σ 2 -12σ + 3 2 √ 6 . (12)
We study the positivity of each polynomial separately, then of both. Notice that the parametrization

ψ 0 = cos α cos β , ψ 2 = sin α cos β , ψ 4 = sin β , -π < α ≤ π, 0 < β ≤ π/2 , ( 13
)
reverses only the sign of ψ 2 if α becomes -α. This parity operation is seen in figure 6, the white domains of which correspond to the positivity of P and Q, respectively. The domain of simultaneous positivity for both is the white intersection domain in the left part of figure 7, with the expected symmetry. together with signatures for the signs of roots, such as P(0), etc. This can be done also in the "spherical representation". The white domains of figures 6, 7 are recovered.

Positivity with monotony

For some problems [START_REF] Kovchegov | [END_REF], it may be useful to request either ψ and/or ϕ to be monotonous functions in an interval such as [0, ∞]. We illustrate this in the case of an H 0 , H 2 , H 4 mixture, with the additional constraint dϕ dσ ∝ -12 cos α cos β+6 √ 2 cos β sin α(2σ-5)-sin β √ 6(4σ 2 -28σ+27) < 0 .

(15) The result appears in the right part of figure 7. The white domain, corresponding to such three simultaneous conditions of positivity and monotonicity, is a severe restriction of the white domain seen in the left part of figure 7.

Positivity from convexity

A practical condition for the positivity of ϕ is the convexity of ψ [START_REF] Lafforgue | private communication; for a general theorem[END_REF]. This has been useful in particular for the derivation of baryon mass inequalities in QCD [8]. Convexity is a sufficient, but not necessary condition. Indeed, ϕ ∼ ∞ 0 dr [1cos(sr)]/s 2 d 2 ψ/dr 2 > 0. We illustrate this convexity condition for a mixture of H 0 , H 4 , H 8 . Clearly, the presence of exp -r 2 /2 in front of a finite order polynomial, with the even parity of ψ and its derivability, are contradictory with "convexity everywhere"; a smooth, round maximum must occur at the origin. We thus study partial convexity conditions of the kind, "convexity between r c and +∞". A reasonable choice for the order of magnitude of r c is the position r c = 1 of the inflexion point of exp -r 2 /2 . The second derivative d 2 ψ/dr 2 belongs to the same algebra. We adjusted its Sturm criterion to various values of r c . Only small domains are found which ensure zero roots for d 2 ψ/dr 2 , because most mixtures of H 0 , H 4 , H 8 do oscillate. Figures 8 and9 show, in white again, with the parametrization by Eqs. ( 9), the survivor domain obtained if r c = 1, right part of figure 8, then r 2 c = 2/3, left part of figure 9, and r 2 c = 2/5, right part of figure 9, respectively. The domain shrinks in a smooth way when r 2 c decreases from 1 to 0.4 and disappears if r 2 c <∼ 0.4. It does not increase much when r c ≥ O(1). The left part of figure 8, a zoom of the left part of figure 6, shows that such partial convexity domains are already included in the positivity domain of ψ.

Positivity for the 2-dimensional Fourier transform

The Fourier-Bessel transform in which we are here interested reads,

ϕ(s) = ∞ 0 (r dr) J 0 (s r) ψ(r) . (16)
For the 2-d radial space, a complete basis of states results from substituting r 2 for r into Laguerre polynomials

∞ 0 (r dr) 2 e -r 2 L m (r 2 ) L n (r 2 ) = δ mn . (17)
For the sake of clarity, we list here the first four such normalized, "2-d radial" states,

{v 0 , v 1 , v 2 , v 3 } = √ 2 e -1 2 r 2 1, r 2 -1, r 4 -4r 2 + 2 2 , r 6
-9r 4 + 18r 2 -6 6 .

(18) One can verify that the states v n make eigenstates of the Fourier-Bessel transform,

∞ 0 (r dr) J 0 (s r) v n (r) = (-1) n v n (s) . (19)
Positivity conditions can again be implemented with the Sturm criterion.

For instance a mixture of v 0 , v 2 , from that subspace with eigenvalue 1, and v 1 , from that with eigenvalue -1, defines the following two reciprocal partners:

P(ρ) = ψ 0 + ψ 1 (ρ -1) + ψ 2 ρ 2 -4ρ + 2 2 , Q(σ) = ψ 0 -ψ 1 (σ -1) + ψ 2 σ 2 -4σ + 2 2 . ( 20
)
The left part of figure 10 shows the positivity domain for ϕ and its right part shows the joint positivity domain for ψ and ϕ. Similarities between polynomials involved in the present "mixed parity" case and those of case C of the previous section create topological similarities with figure 6 (right) and figure 7 (left) but numerical details do differ. Interestingly enough, a sufficient condition [START_REF] Lafforgue | private communication; for a general theorem[END_REF] for the positivity of ϕ in this 2-d situation is that the third derivative, d 3 ψ/dr 3 , be negative. But, as in the previous section, case E, the presence of a truncated number of basis terms in our expansion may reduce this negativity condition for d 3 ψ/dr 3 to a domain r ≥ r t only. We show in figure 11 the result for a mixture of v 0 , v 2 and v 4 if r 2 t = 1.4, a demanding situation resulting into the tiny domain in the right part of figure 11. The domain belongs to the positivity domain seen in the left part of the same figure.

Discussion and conclusion

To summarize our results, we have built a basis of functions verifying positivity together with their Fourier transform. The method is based upon algebras of Hermite polynomials (for 1-dimensional FT) or Laguerre polynomials in the variable r 2 (for 2-dimensional radial FT).

The Fourier transform has four eigenvalues 1, i, -1, -i, and thus four highly degenerate eigensubspaces. Two of such subspaces are compatible with real functions remaining real. To span each subspace, we used a basis made of Hermite-Fourier (or "Laguerre-Fourier") states. The orders (in r) of the polynomials have to be multiples of 4 if the eigenvalue is 1, and multiples of 4 plus 2 if the eigenvalue is -1. At the cost of a truncation of such bases to a maximum order N, the conditions for positivity, convexity, etc. thus reduce to simple manipulations of polynomial coefficients based on the Sturm theorem. In each truncation case, one can find suitable domains for the parameters which mix the various basis polynomials. Such domains have been illustrated by the figures shown in this paper.

Some qualitative considerations may be drawn about the solutions we found. In the left part of figure 12, we display two typical solutions with self-Fourier properties. In the right part, we display one solution connecting two distinct partners ψ and ϕ. Their shapes show distinctive features. One class of shapes, which are monotonic, seem to remain closer to the bare Gaussian, the building block of our method; it always belongs to the subsets we found. The other class, with oscillations, is different. One might ask if our truncations to a maximum order N limit the flexibility of our method to the vicinity of the Gaussian. It seems not to be so, as shown for instance by the very oscillating solution 4 in the left part of figure 12, which exhibits an approximate periodicity in an interval; it is reminiscent of the Dirac comb, which is, of course, outside the set of functions constructed with a finite number of polynomials. However, from our numerical experience, see in particular figures 2 and 3 and Eqs. ( 10) and ( 11), we may risk the conjecture that, when we allow a mixture including H 0 , H 4 , . . . H 4N and when the weight of H 4N is that maximum allowed by positivity, see the tiny red triangle in the right part of figure 3, then the corresponding function ψ = ϕ has a Dirac comb limit when N → ∞. The wavelength of this limit comb is also conjectured to be 1. Combs with different wavelengths may be obtained as limits, but they would not be self-Fourier. For each truncation at a given polynomial order 4N, the locus of such candidates for comb limits can be defined by projecting the combs into the subspace spanned by exp(-r 2 /2) H 0 , exp(-r 2 /2) H 4 , . . . exp(-r 2 /2) H 4N , obviously.

Our bases are flexible enough to reconstruct any function having positivity properties, but in some cases convergence might be slow. It is not excluded that other bases exist, which might be more convenient to speed up the convergence and make easier the search for positivity domains. Another open problem is that of positivity for periodic functions. Such questions are beyond the scope of the present paper.

It is a pleasure to thank R. Balian, R. Enberg, R. Lacaze, C. Marquet, P. Moussa, G. Soyez and A. Voros for stimulating discussions. R.P. thanks T. Lafforgue, (Lycée Blaise Pascal, Orsay), for his fruitful contributions. Thanks are also due to E.H. Lieb and K. Scharnhorst for calling our attention to related works, quoted in the references.

Fig. 1 .

 1 Fig. 1. Mixture of H 0 , H 4 and H 8 . White, no real positive ρ root; black, 4 roots; darkening grays, increasing root number.

Fig. 2 .

 2 Fig. 2. Mixture of H 0 , H 4 , H 8 and H 12 . Root number maps. Color code: 0 root, red; 1, yellow; 2, yellowish green; 3, bluish green; 4, blue; 5, dark purple; 6, pink. Left: cut of the parameter sphere S 3 when γ = π/10 6 . Right: cut for γ = π/15 (uncolored edition: the 0 root domain is the dark triangle-like domain near the center).

Fig. 3 .

 3 Fig. 3. Left: Same as Fig. 2, with γ = 2π/15. Right: γ = π/6. See how the no root triangle shrinks.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Left: γ = π/5; Right: γ = 7π/30. Absence of no root domain.

Fig. 6 .

 6 Fig. 6. Mixture of H 0 , H 2 and H 4 . White domain, 0 root. Grey, 1 root; black, 2 roots. Left: results for P; Right: Q.

Fig. 7 .

 7 Fig. 7. Positivity domains for the parity mixed case. Left: for both P and Q. Right: P, Q and monotony of ϕ.

Fig. 8 .

 8 Fig. 8. H 0 , H 4 , H 8 mixture. Left: zoom of the positivity domain. Right: convexity domain for r > r c = 1.

Fig. 9 .

 9 Fig. 9. Same mixture as Fig. 8. Convexity domains. Left: r 2 c = .67. Right: r 2 c = .4.

Fig. 10 .

 10 Fig. 10. Positivity domains for the "radial 2-d", parity mixed case. Left: ϕ. Right: both ψ and ϕ.

Fig. 11 .

 11 Fig. 11. Left: Positivity domain for the "radial 2-d", v 0 , v 2 , v 4 mixture. Right: third derivative negativity domain if r 2 t = 1.4.

Fig. 12 .

 12 Fig. 12. Left: two examples of self Fourier states (ϕ = ψ). Right: mixed parity case (ϕ = ψ).

4

  Its equation reads ψ = exp(-r 2 /2) h 0.566053 + 0.0488517(3 -12r 2 + 4r 4 ) +0.0011871(105 -840r 2 + 840r 4 -224r 6 + 16r 8 )

	+0.0000164538 "	10395 -124740r 2 + 207900r 4 -110880r 6 + 23760r 8
	-2112r 10 + 64r 12	"i	.

The sequence, made of polynomial remainders of the division of Sn-1 by Sn-2, with (-) signs, clearly stops at some j ≤ N.


J.-P Noel

V Egloff

M Kooli

R Gauchi

J.-M Portal

H.-P Charles

P Vivet

B Giraud

Computational SRAM Design Automation using Pushed-Rule Bitcells for Energy-Efficient Vector Processing

This paper presents a new methodology for automating the Computational SRAM (C-SRAM) design based on off-the-shelf memory compilers and a configurable RTL IP. The main goal is to drastically reduce the development effort compared to a full-custom design, while offering a flexibility of use and a high-yield production. The proposed C-SRAM architecture has been developed to process energy-efficient vector data coupled with a scalar processor, while limiting the data transfer on the system bus. The results obtained by post P&R simulations show that 2RW and 4RW C-SRAM configurations using the double pumping technique achieved the highest performance to process vectorized MAC operations compared to the others configurations. Moreover, it has been shown that the impact of the digital wrapper decoding and executing the instructions can be mitigated by increasing the memory cut size to represent less than 10% in area and 20% in power consumption.

I. INTRODUCTION

For edge AI applications processing a large amount of vector data, parallel computing is the straight forward solution. However, due to the transfer of large vector data (up to 512 bits) through the system bus, current commercial vector processors (e.g. SIMD) are no longer enough efficient to perform lowenergy operations. Indeed, it has been shown that the cost of data transfer between the memories and processing elements is much higher than the operation itself in high-performance computer architectures for advanced technology nodes [1][START_REF] Hennessy | Computer Architecture: A Quantitative Approach[END_REF]. In this context, In-and Near-Memory Computing (IMC/NMC) architectures defined by a limited data transfer on the system bus seem to be particularly well suited. In the recent literature, several Computational SRAM (C-SRAM) architectures, in which the computing is performed in the bitcell matrix (IMC) and/or at peripheral circuit level (NMC), have been proposed. Most of them, based on silicon results, are full-custom solutions, proposing specific peripheral circuits and/or bitcells complying with logic design rules. This has the advantage of maximizing the energy efficiency of computation, up to several tens of TOPS/W [START_REF] Zhang | In-Memory Computation of a Machine-Learning Classifier in a Standard 6T SRAM Array[END_REF][START_REF] Zhang | Recryptor: A Reconfigurable Cryptographic Cortex-M0 Processor With In-Memory and Near-Memory Computing for IoT Security[END_REF][START_REF] Wang | A Compute SRAM with Bit-Serial Integer/Floating-Point Operations for Programmable In-Memory Vector Acceleration[END_REF]. Nevertheless, this design methodology presents two major drawbacks: difficult silicon qualification and low adaptability to specific applications requirements. First, to qualify with a high-yield production a custom-designed C-SRAM IP, specific test circuits (e.g. BIST) need to be developed, resulting in major changes to the well-established validation flow. Then, to expand the operating range, several memory sizes need to be qualified, which also requires a major development effort. Finally, to meet EDA tool compatibility requirements in a time-limited design phase, automated generation of C-SRAM views (HDL, physics, timing and power...) need also to be developed. Finally, it is necessary to easily manage a large number of configurations (capacity, form factor, operations...), while guarantying a high-yield production.

In this paper, we propose a new methodology for automating the C-SRAM design based on off-the-shelf SRAM compilers for the storage part and a configurable digital wrapper (RTL IP) for the computing part. The main contributions of the paper are:  To evaluate the impact, for each selected memory type, the energy efficiency of a Multiply-Accumulate (MAC) operation performed in the digital wrapper.

 To quantify the additional area and energy cost of the computing part (RTL IP).

This work prefigures the algorithm that could be implemented in the automated decision-making system (ADMS) to select the optimal memory type according to the user's constraints (e.g. memory capacity, instruction frequency, required operations ...).

In order to take full advantage of the parallel computing based on the use of C-SRAM as a vector processing accelerator, the code compiler must also be optimized to execute as many vector-based kernels as possible of a given algorithm. This part is out of the scope of this paper.

The remainder of the paper is organized as follow: Section II demonstrates the benefits of using a C-SRAM as an energyefficient vector processing accelerator over conventional SIMD processors. Section III describes the proposed methodology for automating the C-SRAM design based on off-the-shelf SRAM compilers and a configurable RTL IP. Section IV presents the simulation results (post P&R) highlighting the correlation between the optimum energy efficiency to perform MAC operations and the number of ports in the selected memory. Finally, Section V summarizes the paper and presents the perspectives of applications.

II. COMPUTATIONAL SRAM: AN ENERGY-EFFICIENT VECTOR PROCESSING ACCELERATOR

The C-SRAM architecture presented in this paper is designed to be implemented with a scalar CPU as an energy-efficient vector processing accelerator thanks to the reduced transfer of vector data on the system bus [START_REF] Gauchi | Exploration of a Scalable Vector-based In-Memory Computing Architecture via a System-on-Chip Evaluation Framework[END_REF]. To highlight this assumption, Figure 1 provides a comparison between scalar, vector (e.g. a SIMD processor) and scalar/vector (based on C-SRAM). Compared to the scalar architecture, the vector architecture proposes parallel computations on vector data (up to 512-bit for a commercial SIMD processor).

However, it requires data transfers over the system bus towards the registers of the processor. This drawback can be solved by the scalar/vector architecture based on C-SRAM by limiting the data transfer, while keeping parallel computation on vector data. Most of vector operations could be thus directly performed in the C-SRAM (containing the data to be processed) after receiving the instructions from the CPU (up to 64-bit) instead of transferring the data vectors (128-bit or more) to the CPU registers [START_REF] Kooli | Smart Instruction Codes for In-Memory Computing Architectures Compatible with Standard SRAM Interfaces[END_REF]. As for conventional vector processors, the size of C-SRAM operators can be adjusted according to the element length (8-bit, 16-bit, 32-bit ...) of the vector data. The size of the latter are used to define the length of C-SRAM words (128-bit, 256-bit ...).

III. PROPOSED C-SRAM DESIGN METHODOLOGY

SRAM compiler re-use & digital wrapper

To take advantage of the configurability (size, form factor, power management, ECC ...) and the reliability (silicon qualification) of existing SRAM compilers, while limiting the development efforts, we propose to use a digital wrapper based on a configurable RTL IP (Fig. 2). The wrapper has to decode and execute the C-SRAM instructions sent by the CPU without transferring vector data through the system bus. In addition, this wrapper also enables the data vectorization coming from the system bus (up to 64-bit) to the C-SRAM (from 128-bit), and vice-versa. The customizable parameters are the size of the input (depending on the system bus) and the vector data (depending on the application). Then, several operations can be selected from a pre-defined list of a specific Instruction Set Architecture (not described in this paper). Finally, the management of the wrapper pipeline depends on the timing constraints imposed by the frequency difference between the instructions sent by the CPU and the memory accesses.

Another advantage to use existing SRAM compilers and wrap the generated instances with a glue logic is to enable a full flexibility in terms of memory partionning. Due to performance constraints (frequency, power, ...), the memory word length can be limited above the targeted data vector length (e.g. 512-bit). In this case, several memory instances with a limited word length (e.g. 32-or 64-bit) can be used to reach the targeted data vector length, while sharing the same digital wrapper to limit area and leakage power penalties [START_REF] Gauchi | Memory Sizing of a Scalable SRAM In-Memory Computing Tile Based Architecture[END_REF]. A similar partionning can be done to overcome the limitation of the number of words per SRAM instance, as shown in Fig. 3.

To simplify the development, we propose to divide the C-SRAM design in two parts: storage and computing (Fig. 4). The storage part is based on the use of SRAM compilers to generate the memory cuts. The computing part is based on the use of a configurable RTL IP to implement the digital wrapper. To select the appropriate memory type according to the computing features coming from the user's requests, we propose to use an ADMS. It is based on rules from the results presented in section IV. The algorithm on which it is based is not detailed in this paper.

The outputs of the proposed C-SRAM design flow enable the placement and routing of the digital wrapper with the conventional EDA tools (step 2 in Fig. 4). Nevertheless, the development of this design flow could be pushed one step further to generate an assembled hardware macro with all the EDA views necessary for its implementation (step 3 in Fig. 4).

Multi-port SRAM compiler possibilities based on pushed-rule foundry bitcells

Three types of SRAM bitcells are systematically densityoptimized by the founders: Single-Port (SP), Dual-Port (DP) and Two-Port (TP). These pushed-rule bitcells are optimally designed then qualified on silicon to meet the performance requirements of the process technology. Next, several SRAM compiler types are usually developed, each of them for a specific use. These compilers can be classified by the number of access ports, and more specifically by the number of read / write accesses in a same clock cycle. For example, a 1RW SRAM compiler means one memory access in one clock cycle, either to write or read a data at a specific address.

To increase the number of synchronous memory accesses, the use of DP or TP bitcells has been generalized. Typically, DP bitcells are used for video processing where multi-access memories (e.g. 2RW) optimize the data-parallel processing [START_REF] Nii | 2RW Dual-port SRAM Design Challenges in Advanced Technology Nodes[END_REF]. While TP bitcells are commonly used as register files where simultaneous read and write accesses (e.g. 1R1W) are needed for high-speed computing [START_REF] Kulkarni | 5.6 Mb/mm2 1R1W 8T SRAM Arrays Operating Down to 560 mV Utilizing Small-Signal Sensing With Charge Shared Bitline and Asymmetric Sense Amplifier in 14 nm FinFET CMOS Technology[END_REF]. In any case, the use of this type of bitcells results in a significant area and leakage overhead (compared to SP-based SRAM). Moreover, to increase the number of ports beyond 2, specific multi-port bitcells are required, further exacerbating this issue. To develop multi-port SRAM compilers without designing specific bitcells or limiting the area and leakage overhead, a circuit design technique, called double pumping, has been adopted by most of the SRAM compiler vendors in advanced technology nodes. This technique consists in artificially duplicating the number of ports by generating a second internal memory clock to start a consecutive memory access without waiting for the end of the first clock. This is made possible by modifiying decoding and IO circuitry in order to add intermediate sequences of dynamic logic, while reducing the operating frequency (up to ~40%) [START_REF] Ditlow | A 4R2W Register File for a 2.3GHz Wire-Speed POWER Processor with Double-Pumped Write Operation[END_REF][START_REF] Wu | A Configurable 2-in-1 SRAM Compiler with Constant-Negative-Level Write Driver for Low Vmin in 16nm Fin-FET CMOS[END_REF][START_REF] Yabuuchi | A 6.05-Mb/mm² 16-nm FinFET Double Pumping 1W1R 2-port SRAM with 313ps Read Access Time[END_REF][START_REF] Ishii | A 5.92-Mb/mm² 28-nm Pseudo 2-Read/Write Dual-port SRAM using Double Pumping Circuitry[END_REF][START_REF]An Ultra High Density Pseudo Dual-Port SRAM in 16nm FINFET Process for Graphics Processors[END_REF]. Nowadays, most of founders and IP vendors propose 2-Port SRAM compilers based on SP bitcells that use this technique to improve both memory density and leakage power consumption at the expense of operating frequency compared to the reference solution (Fig. 5). The more the double pumping technique is optimized, the lower the frequency loss. This technique can also be advantagously used with DP and TP bitcells to achieve up to 4-Port SRAM compilers enabling 4RW, 2R2W or 2R2RW memory accesses [START_REF] Nguyen | A 7nm Double-Pumped 6R6W Register File for Machine Learning Memory[END_REF].

Figure 6 exposes the range of possibilities of SRAM compiler types based on the conventional pushed-rule bitcells (SP, DP and TP) and the use of the double pumping technique for duplicating the number of ports. We will see in the next section the benefits of use 4-Port SRAM compilers, notably those developed from the double pumping technique.

IV. SIMULATION RESULTS

In this section, we quantify the impact of several SRAM types on the energy efficiency of vector processing. Furthermore, we also quantify the additional area and energy cost of the digital wrapper for several sizes of C-SRAM.

Methodology

This study is based on the 22nm FD-SOI (22FDX) design platform developed by Globalfoundries. 8-Track standard cell libraries are used (in order to minimize the area overhead) to design the digital wrapper and a set of representative SRAM compilers (all using pushed-rule bitcells described in the previous section) are selected for the storage part of the C-SRAM macros (1RW, 2RW, 2RW-DP and 4RW-DP). Simulation results are obtained from post P&R gate-level netlists close to typical conditions (TT/0.8V/85°C).

Testcase: Multiply-Accumulate Operation

To make the study relevant to edge AI applications, we have decided to base our testcase on a vectorized (16-element of 8bit) MAC operation. This instruction is executed by pipelining the breakdown operations performed either in the memory (read & write) or in the digital wrapper (multiply & addition). Figure 7 describes the chronogram of the vectorized MAC operation for different memory types of C-SRAM macros. The breakdown operations of this instruction are made of 4 memory accesses (RD A, RD B, RD C, WR Z), 1 instruction decoding, 1 multiplication (A x B) and 1 addition (A x B + C). In this testcase, the operations performed in the digital wrapper are parallelized as much as possible with the memory access clocked at the same frequency (fC-SRAM). The duration of the MAC instruction (latency) can vary between 5 to 6 clock cycles. This variation corresponds to the ability to get A and B operands in the same clock cycle or not, which is the case for 2RW and 1R1RW memories (5 clock cycles) but not for 1RW and 1R1W memories (6 clock cycles). This figure also describes the chronogram corresponding to the use of a pipeline in the digital wrapper, introducing the use of the 4RW memory. This technique enables to execute successive instructions without waiting for the end of the previous one. Therefore, it is possible to reduce the instruction throughput from 5 cycles to only one. Table 1 describes the different configurations (memory size, maximum operating frequency …) of the C-SRAM macros selected for this study. Note that he results obtained in this section for 2RW-DP and 4RW-DP C-SRAM configurations have been extrapolated from simulation of 1RW and 2RW C-SRAM configurations based on [START_REF] Ditlow | A 4R2W Register File for a 2.3GHz Wire-Speed POWER Processor with Double-Pumped Write Operation[END_REF][START_REF] Wu | A Configurable 2-in-1 SRAM Compiler with Constant-Negative-Level Write Driver for Low Vmin in 16nm Fin-FET CMOS[END_REF][START_REF] Yabuuchi | A 6.05-Mb/mm² 16-nm FinFET Double Pumping 1W1R 2-port SRAM with 313ps Read Access Time[END_REF][START_REF] Ishii | A 5.92-Mb/mm² 28-nm Pseudo 2-Read/Write Dual-port SRAM using Double Pumping Circuitry[END_REF][START_REF]An Ultra High Density Pseudo Dual-Port SRAM in 16nm FINFET Process for Graphics Processors[END_REF][START_REF] Nguyen | A 7nm Double-Pumped 6R6W Register File for Machine Learning Memory[END_REF]. The maximum operating frequency varies between 1.5GHz to 2.2GHz for 4RW-DP and 2RW cuts of 128-word (minimum size) of 128bit, respectively. The frequency decrease down to 25% for the maximum size of cut (16k-word for 1RW and 4k-word for 2RW and 4RW). Based on these data, it is possible to get the maximum number of vectorized MAC operations performed by a C-SRAM macro, knowing that 16 MAC operations of 8-bit are performed in parallel (Fig. 8). Duplicating the number of ports of the 1RW memory and pipelining the instruction execution in the digital wrapper increase the performance up to 2.7x. Furthermore, applying the double pumping technique to the 2RW memory to get the 4RW C-SRAM configuration enables to achieve a speedup up to 3.6x. Nevertheless, the use of bigger (128-to 16k-words) and slower macro reduce these performance down to 22%. Figure 9 shows in detail the computing part (included in the digital wrapper) of C-SRAM macros in terms of area and power ratio. As expected, the biggest ratio is achieved for the smallest memories (area: 45% and power: 37%), in particular for those designed with the double pumping technique (area: 51% and power: 37%). The computing part can be reduced by using bigger memory size down to 10% and 20% in area and power ratio, respectively. The physical implementation and the memory/computing partitioning are illustrated in Figure 10 for the 1RW and 2RW C-SRAM macros representing the smallest and the biggest size allowed by the memory compilers. Figure 11 shows the relative area cost of 2RW, 2RW-DP or 4RW-DP memories compared to a 1RW memory. As expected, the 2RW and 4RW using the Table 1 Features of selected C-SRAM macros in terms of size (word length and number) and speed (max. operating frequency and throughput/latency of instructions). double pumping technique enable to limit this cost below 2x (vs. 3x) and 4x, respectively. Figure 12 shows the number of operations (TOPS) for different C-SRAM configurations with a budget of 1W and 1mm² according to the instruction frequency. This illustrates that both the best tradeoff (TOPS/W/mm²) and the highest instruction frequency are achieved by the 4RW-DP C-SRAM configurations. Regarding the configurations executing sequentially the instructions (i.e. w/o pipeline), the 1RW remains the best solution, with a loss of the instruction frequency of 26% for 2.9x of performance w.r.t. 2RW solutions. In all cases, the performance are drastically reduced (down to 33x) by using the biggest memory size (16kwords).

V. CONCLUSION

This paper proposed a design methodology for building energyefficient C-SRAM macros from on-the-shelf memory compilers (based on pushed-rule foundry bitcells) and a configurable RTL IP. The benefit of this approach is to minimize development efforts and production risks using proven design techniques (i.e. double pumping) to take advantage of multi-port SRAM compilers (up to 4-Port). It was demonstrated that this C-SRAM configuration is the most energy efficient to perform vectorized MAC operations (16element of 8-bit) when the instructions sent by the CPU can be pipelined (no data dependencies from one cycle to the next). Otherwise, the 1RW C-SRAM configuration remains the best choice. It was also demonstrated that increasing the memory capacity significantly reduces the energy efficiency of vectorized MAC operation per unit area (up to 33x). The good tradeoff depends on the user's constraints in terms of area and power budgets as well as system performance (instruction frequency…). The next step in this work is to propose an ADMS to build custom C-SRAMs that match each specific request. This design automation methodology paves the way for the integration of energy-efficient vector processing accelerators into energy-constrained systems using scalar processors (microcontrollers ...).

Figure 3

 3 Figure 3 128-bit SRAM cut partitioning constrained by the physical limitations: a) word & bit number < physical limits, b) word number > physical limits, c) bit number > physical limits and d) word & bit number > physical limits.

Figure 1

 1 Figure 1 Illustration of scalar (left), vector (middle) and scalar/vector (right) computing architectures. C-SRAM (right case) enables to drastically reduce the data transfer on the system bus (leading to significant energy savings) based on near-memory computing.

Figure 2

 2 Figure 2 Proposed design methodology to automate the C-SRAM macro generation w.r.t. full custom design solutions. Automation simplifies the design phase but makes it more difficult to implement analog or precomputing functions.

Figure 4

 4 Figure 4 Proposed C-SRAM design flow based on off-the-shelf SRAM compilers and a configurable RTL IP dedicated to decode and execute specific instructions coming from the CPU.

Figure 5

 5 Figure 5 Basic concept of the double pumping technique enabling the port number duplication, while keeping a clock frequency ratio ≥ 0.6 w.r.t. the reference solution.

Figure 7

 7 Figure 7 Chronogram of the multiply-accumulate operation (Z = A * B + C) performed with various C-SRAM configurations. These configurations depend on the selected SRAM type (1RW, 1R1W ...), the frequency ratio between C-SRAM and instruction flow and the way the instructions are performed in the digital wrapper (pipeline enabled or not).

Figure 6

 6 Figure 6 Possible types of SRAM compilers developed from the 3 main pushed-rule foundry bitcells (Single-Port, Dual-Port and Two-Port) and the use of the double pumping technique for duplicating the number of port.

Figure 8

 8 Figure 8 Number of vectorized MAC operations performed by a C-SRAM macro during 1 second (GOPS) for various memory types and sizes (min. and max. macro size).

  maximum C-SRAM frequency (GHz) / throughput (cycle #) * 16-operations

Figure 9 Figure 10

 910 Figure 9 Computing part (area and power ratio) in a C-SRAM macro for various memory types and sizes (min. and max. macro size).

Figure 12 Figure 11

 1211 Figure 12 Energy efficiency of vectorized MAC operations per area unit (TOPS/W/mm²) performed by a C-SRAM macro according to the instruction frequency for various memory types and sizes (128-and 16kwords).