N Lahellec

P Ponte Castañeda

P Suquet
email: suquet@lma.cnrs-mrs.fr

Variational estimates for the effective response and field statistics in thermoelastic composites with intra-phase property fluctuations

Composites, Thermoelasticity, Variational methods. In this work, variational estimates are provided for the macroscopic response, as well as for the first and second moments of the stress and strain fields, in thermoelastic composites with nonuniform distributions of the thermal stress and elastic moduli in the constituent phases. These estimates are obtained in terms of a 'comparison composite' with uniform phase properties depending on the first and second moments of a certain combination of the given intra-phase thermal stresses and modulus field distributions. Under certain hypotheses, these estimates can be shown to lead to upper and lower bounds for the free energy of the composite, which reduce to standard results when the intra-phase fluctuations vanish. An illustrative application is given for rigidly reinforced composites with a non-uniform distribution of the thermal stress in the matrix phase.

Introduction

There are many problems of current interest which can be described as 'thermoelastic' composites, where the 'thermal' stress (also called 'eigenstress' or 'transformation stress' in the literature), or the elasticity tensor, or both, need not be uniform within a given phase. One standard example is provided by residual stresses in polycrystals [START_REF] Noyan | Characterization of Inhomogeneous Elastic Deformation with X-Ray Diffraction[END_REF], where due to the forming process, non-uniform distributions of the stresses may develop within the grains of the polycrystal, even when the applied macroscopic stress vanishes. A second example is provided by elasto-(visco)plastic composites (see Lahellec & Suquet (2007a)) where the (visco)plastic strain field may be heterogeneous in a given phase. As a final example, it is noted that the 'incremental' linearization scheme proposed by [START_REF] Hill | Continuum micro-mechanics of elastoplastic polycrystals[END_REF] for estimating the macroscopic behavior of plastic polycrystals involves the 'tangent' modulus of the grains, which is non-uniform in the grains, due to the inhomogeneity of the stress and strain fields that are produced in the grains of the polycrystal under loading.

Motivated by these examples, we consider a representative volume element V of a composite material made of different linear thermoelastic phases. The stress-strain relation in an individual phase reads as:

σ(x) = L(x) : ε(x) + τ (x), (1)
where both the elastic moduli L(x) and the polarization field τ (x) depend on the position x and are non uniform per phase.

The question addressed in this study is to estimate the macroscopic behavior, as well as the first and the second moments of the strain and stress fields, ε and σ, in a thermoelastic composite with phase constitutive response described by relation (1)-utilizing only partial statistical information about L and τ in each phase (assumed to be statistically uniform in the composite). To the best of the authors' knowledge, the exact response to this question is not known in general. However it is known in the special case where the elastic moduli L and the polarization field τ are uniform in each phase r, i.e. piecewise uniform through the entire volume element V . Quite naturally, our approach to the problem is to replace the initial thermoelastic composite with nonuniform moduli and nonuniform polarization field by a linear (thermoelastic) comparison composite (LCC) with piecewise uniform elastic moduli L (r) 0 and piecewise uniform polarization stress τ (r) 0 .

To this end, a very natural first choice is to replace the nonuniform fields τ (x) and L(x) by their average, or first moment, over each individual phase: r) .

L (r) 0 = L (r) , τ (r) 0 = τ (
(2) However, even in the one-dimensional case where calculations can be carried out explicitly (see appendix), this simple choice does not yield the exact result for the second moment of the stress and strain fields.

The criterion adopted here to determine the improved estimates for the quantities L (r) 0 and τ (r) 0 in the comparison composite is based on a variational approach for a 'comparison composite,' which follows closely the variational approach of Ponte Castañeda (1991) for nonlinear composites (see also [START_REF] Ponte Castañeda | Nonlinear composites[END_REF] and Ponte Castañeda & [START_REF] Ponte Castañeda | Variational second-order estimates for nonlinear composites[END_REF], and references therein).

Effective energies and elementary bounds

Select a volume element V of the nonhomogeneous material under consideration. Although neither the elastic moduli or the polarizations are piecewise uniform, we shall assume that there exists an underlying geometrical decomposition V in the form of N subdomains V (r) | r=1,..,N which is a natural microstructure. An example of such a situation corresponds to a polycrystalline material where the grains are well identified, even though the plastic strain γ can be nonuniform even within each grain. The characteristic function and the volume fraction of V (r) are denoted by χ (r) (x) and c (r) , respectively. The angular brackets,

. and . (r) , will be used to denote the spatial averages over V and V (r) , respectively. Total and partial averages will be denoted by an upper bar:

ε = ε , ε (r) = ε (r)

Effective potentials

The constitutive relation (1) derives from a (convex) thermoelastic energy w which reads at point x as:

w(x, ε) = 1 2 ε : L(x) : ε + τ (x) : ε. (3)
The local stress-strain relation (1) can be inverted and derives from a dual potential denoted by u, conjugate to w in the sense of convex dual functions:

u(x, σ) = sup ε {σ : ε -w(x, ε)} . (4)
A straightforward calculation shows that:

u(x, σ) = 1 2 σ : M (x) : σ + γ(x) : σ + 1 2 γ(x) : L(x) : γ(x), (5)
where

M (x) = L -1 (x), γ(x) = -M (x) : τ (x).
(6) Although the elastic moduli and the polarizations are not piecewise uniform, it will be convenient to use a different notation for the restriction of w and u to V (r) :

w (r) (x, ε) = w(x, ε)χ (r) (x), u (r) (x, σ) = u(x, σ)χ (r) (x).
The effective potential, or energy function, corresponding to w is given by

w(ε) = inf ε ∈K(ε) w(x, ε) , (7)
where K(ε) is the set of fields ε(x) such that there exists a displacement field v(x) satisfying

ε(x) = 1 2 (∇v + T ∇v) and v(x) = ε.x on ∂V. (8)
The complementary effective potential u is characterized as:

u(σ) = inf σ ∈S(σ) u(x, σ) , (9)
where S(σ) is the set of fields σ(x) such that div (σ) = 0, σ = σ.

Note that u and w are dual convex potentials [START_REF] Suquet | Elements of Homogenization for Inelastic Solid Mechanics[END_REF] in the linear case and [START_REF] Willis | The structure of overall constitutive relations for a class of nonlinear composites[END_REF] in general) satisfying in particular:

u(σ) = sup ε {σ : ε -w(ε)} . ( 10
)
Note also that the macroscopic stress σ and macroscopic strain ε are related by the constitutive relations

σ = ∂ w ∂ε (ε) , and ε = ∂ u ∂σ (σ) . ( 11
)

Elementary bounds

The Voigt and Reuss bounds correspond to the particular choice ε = ε in ( 7) and σ = σ in (9) respectively:

w(ε) ≤ w Voigt (ε) = 1 2 ε : L : ε + τ : ε, u(σ) ≤ u Reuss (σ) = 1 2 σ : M : σ -σ : M : τ + 1 2 τ : M : τ .
A straightforward calculation gives the dual of u Reuss as:

w Reuss (ε) = 1 2 ε : M -1 : ε + ε : M -1 : M : τ + 1 2 M : τ : M -1 : M : τ -1 2 τ : M : τ .
These bounds can be given a general form:

w(ε) = 1 2 ε : L : ε + τ : ε + f (τ ), (12)
where for the Voigt bound:

L = L Voigt = L , τ = τ Voigt = τ , f (τ ) = f Voigt (τ ) = 0, (13)
and for the Reuss bound:

L Reuss = M -1 , τ Reuss = M -1 : M : τ , f Reuss (τ ) = 1 2 τ : M : M -1 : M : τ -1 2 τ : M : τ .      (14)

General structure of the exact energy

We now come back to the general problem

σ(x) = L(x) : ε(x) + τ (x), div (σ(x)) = 0, ε ∈ K(ε). (15)
The tensors ε and τ (x) being given, the solution ε of this problem can be obtained by the superposition principle. Consider first the case where τ is identically 0. Problem ( 15) is then a standard elasticity problem and its solution can be expressed by means of the elastic strain-localization tensor A(x) as

ε(x) = A(x) : ε. ( 16
)
Consider next the case where ε = 0 and τ (x) is arbitrary. The solution of (15) can be expressed by means of the nonlocal elastic Green operator Γ(x, x ) of the nonhomogeneous elastic medium as:

ε(x) = -Γ * τ (x), where Γ * τ (x) def = 1 |V | V Γ(x, x ) : τ (x ) dx . ( 17
)
It follows from the superposition principle that the solution of (15) reads as:

ε(x) = A(x) : ε -Γ * τ (x). (18)
Using this form of the strain field in the effective energy

w(ε) = 1 2 ε : L : ε + ε : τ ,
yields an expression which can be simplified by noting that:

τ : Γ * τ = (Γ * τ ) : L : (Γ * τ ) , (19)
and (Γ * τ ) : L :

A : ε = 0. ( 20
)
The relation ( 19) results from the fact that Γ * τ is a compatible strain field with average 0 and that -L : Γ * τ + τ is a divergence-free stress field. Therefore, it follows from Hill's lemma that: 19). Similarly, the relation (20) results from the fact that L : A : ε is divergencefree stress field and that Γ * τ is a compatible strain field with average 0. In addition, since A(x) : ε is a compatible strain field, and A T = I, it follows that

(Γ * τ ) : (-L : Γ * τ + τ ) = 0, hence (
ε : A T : L : A : ε = ε : A T : L : A : ε = ε : L : A : ε. (21)
Then, using the relations ( 19), (20), and (21), the effective energy w can be put in the form (12) with:

L = L : A , τ = A T : τ , f (τ ) = - 1 2 τ : Γ * τ . (22)
It is useful to note that, since Γ * τ is a strain field with average 0, one has, for every divergence-free stress field σ * :

τ : Γ * τ = (τ -σ * ) : Γ * (τ -σ * ) . ( 23
)
This last relation shows that f is in fact due to the part of τ which is not in equilibrium.

Remarks:

(1) It follows from ( 22) and ( 19) that

f (τ ) = - 1 2 (Γ * τ ) : L : (Γ * τ ) ≤ 0. ( 24
)
In other words, f is always negative. In addition f = 0 if and only if Γ * τ = 0, which is the case only when τ is a divergence-free field.

(2) Assume that τ is not a divergence-free field. Taking λτ in the effective energy (7) and letting λ go to +∞ (ε being fixed), it is seen that the exact effective energy goes to -∞ as λ 2 f (τ ).

(3) When ε = 0, the Voigt and Reuss bounds read:

f Reuss (τ ) ≤ f (τ ) ≤ f Voigt (τ ). ( 25
)
The term f Voigt (τ ) is identically 0 and the second inequality in (25) is nothing else than the fact that f (τ ) is always negative (already derived by different means in (24)). On the other hand , the Reuss bound does reproduce the divergence of the f (λτ ) as λ goes to +∞. In this connection, note that the energy term f Reuss (τ ) can be alternatively written

f Reuss (τ ) = - 1 2 τ * : M : τ * , τ * = τ -M -1 : M : τ = τ -τ Reuss .
Thus, the Reuss bound on w can diverge as -λ 2 when λ goes to +∞, except when τ * = 0 which corresponds to a completely uniform polarization τ . (4) It is recalled [START_REF] Laws | On the thermostatics of composite materials[END_REF]; [START_REF] Willis | Variational and related methods for the overall properties of composites[END_REF]) here for future reference that when L(x) and τ (x) are constant per phase such that r) , and τ

L(x) = N r=1 χ (r) (x)L (
(x) = N r=1 χ (r) (x)τ (r) , (26)
where L (r) and τ (r) are constant, the effective properties (22) associated with the effective energy (12) reduce to

L = N r=1 c (r) L (r) : A (r) , τ = N r=1 c (r) A (r) T : τ (r) , f (τ ) = 1 2 N r=1 c (r) a (r) : τ (r) , (27)
where A (r) = A (r) and a (r) = -Γ * τ (r) are strain concentration tensors such that the average of the strain over phase r is given by ε r) . In particular, for a two-phase composite, it is known [START_REF] Levin | Thermal expansion coefficients of heterogeneous materials[END_REF]) that τ and f may be obtained directly in terms of L via

(r) = A (r) ε + a (
τ = τ + ( L -L ) : L (1) -L (2) -1 : τ (1) -τ (2) , (28)
f (τ ) = 1 2 τ (1) -τ (2) : L (1) -L (2) -1 : ( L -L ) : L (1) -L (2) -1 : τ (1) -τ (2) . ( 29
)
3 Bounds on the effective potentials through linear comparison composites

Lower bound on the primal energy

Introduce a linear comparison composite (LCC) with piecewise constant elastic moduli and polarization:

w 0 (x, ε) = 1 2 ε : L 0 (x) : ε + τ 0 (x) : ε, (30)
where

L 0 (x) = N r=1 χ (r) (x)L (r) 0 , τ 0 (x) = N r=1 χ (r) (x)τ (r) 0 , (31)
and let w (r)

0 (ε) = 1 2 ε : L (r) 0 : ε + τ (r) 0 : ε. ( 32
)
Assume that

∆L (r) (x) = L(x) -L (r) 0 ≥ 0, (33)
where the inequality in ( 33) is to be understood in the sense of positive semi-definite quadratic forms and set:

∆τ (r) (x) = τ (x) -τ (r) 0 . (34)
Consider the infimum problem:

inf e∈R 3×3 s w (r) (x, e) -w (r) 0 (e) .
Note that due to assumption (33), the function w (r)w (r) 0 is convex .Therefore, the above infimum problem has a solution and one can define the 'error' function V (r) between w (r) and w (r) 0 as:

V (r) (x, τ (r) 0 , L (r) 0 ) = inf e∈R 3×3 s w (r) (x, e) -w (r) 0 (e) = - 1 2 ∆τ (r) (x) : (∆L (r) (x)) -1 : ∆τ (r) (x). ( 35
)
Then, for every τ (r) 0 and L (r) 0 satisfying (33), one has:

w (r) (x, ε) ≥ w (r) 0 (ε) + V (r) (x, τ (r) 0 , L (r) 0 ). ( 36
)
By definition, w(ε) is such that

w(ε) = inf ε ∈K(ε) w(x, ε) = inf ε ∈K(ε) N r=1 c (r) w (r) (x, ε) (r) .
According to (36), the following inequality holds for every ε ∈ K(ε) and for every τ satisfying (33):

w(x, ε) ≥ N r=1 c (r) w (r) 0 (ε) + V (r) (x, τ (r) 0 , L (r) 0 ) (r) . ( 37
)
Taking the infimum of the two sides of (37) with respect to ε, one has

w(ε) ≥ w 0 (ε) + V(x, τ 0 , L 0 ) , (38)
for every τ

(r) 0 , L (r)
0 satisfying (33), where

w 0 (ε) = inf ε ∈K(ε) w 0 (x, ε) = N r=1 c (r) w (r) 0 (ε) (r) , and
V(x, L 0 , τ 0 ) = N r=1 c (r) V (r) (x, τ (r) 0 , L (r) 0 ) (r) . (39)
Finally, making use of ( 35) in (39), and taking the supremum of ( 38) over all possible τ (r) 0

and L (r) 0 , one obtains the result:

w(ε) ≥ sup τ (r) 0 ,L (r) 0 w 0 (ε) - 1 2 N r=1 c (r) ∆τ (r) (x) : (∆L (r) (x)) -1 : ∆τ (r) (x) (r) . (40)
Stationarity conditions may be written for the supremum problem (40). Thus, stationarity with respect to τ (r) 0 gives:

ε 0 (r) = -∆L (r) -1 : ∆τ (r) , (41)
where ε 0 denotes the strain field in the LCC, while stationarity with respect to L (r) 0 leads to

ε 0 ⊗ ε 0 (r) = ∆L (r) -1 : ∆τ ⊗ ∆L (r) -1 : ∆τ (r) . ( 42
)
The two equations ( 41) and ( 42) are to be solved for τ 0 | r=1,..,N , keeping in mind that ε 0 depends on these unknowns. This system of nonlinear coupled equations may have one solution, or several solutions, or no solution.

Let us assume that there exists a solution for the L (r) 0 's and τ (r) 0 's which satisfies the assumed positivity condition (33). In order to have the best possible lower bound one has to make sure that the stationarity point in ( 40) is a maximum. For this purpose, it is sufficient to prove that the function w 0 (ε) + V(x, τ 0 , L 0 ) is a concave function of the unknowns (τ

(s) 0 , L (s) 0 )| s=1,..,N .
To show this, first note that, for every ε, w (r) 0 (ε) is an affine, and therefore concave, function of (L (r) 0 , τ (r) 0 ). Similarly it follows from (35) that, under the inequality (33), V (r) is a concave function of (τ

(r) 0 , L (r) 0 ). Consequently the function w (r) 0 (ε) + V (r) (x, τ (r) 0 , L (r) 0 ) is, for every ε, a concave function of (L (r) 0 , τ (r) 0 ). Therefore, N r=1 c (r) w (r) 0 (ε) + V (r) (x, τ (r) 0 , L (r) 0 ) (r)
is a concave function of the unknowns (τ

(s) 0 , L (s) 0 )| s=1,.
.,N . Then taking the infimum over ε of these concave functions is again a concave function which proves that w 0 (ε) + V(x, τ 0 , L 0 ) is a concave function of the unknowns (τ

(s) 0 , L (s) 0 )| s=1,.
.,N . Therefore, every stationary point for the right-hand-side of ( 40) is a supremum provided that the solutions to ( 41) and ( 42) satisfy the inequality (33).

Comments:

(1) An upper bound for the energy can be obtained in a similar way when

∆L (r) (x) = L(x) -L (r) 0 ≤ 0. ( 43
)
The modifications to be made to the lower bound procedure go as follows. First the infimum in ( 35) is changed into a supremum. Therefore the inequalities in (36), ( 37) and ( 38) are reversed. The bound (40) becomes an upper bound:

w(ε) ≤ inf τ (r) 0 ,L (r) 0 { w 0 (ε) + V(x, τ 0 , L 0 ) } , (44)
and the stationarity conditions ( 41) and ( 42) remain unchanged. In order to prove that the stationary points of the right-hand-side of ( 44) are in fact infima, it is sufficient to show that w 0 (ε) + V(x, τ 0 , L 0 ) is a convex function of (τ

(s) 0 , L (s) 0 )| s=1,.
.,N . The arguments follow closely those used to prove that it is a concave function of these arguments when the inequality (33) holds.

(2) Define the fluctuations C (r) (e) of a field e in phase r as: r) ⊗ ee (r) (r) .

C (r) (e) = e -e (
(45)

Then, making use of (41), the stationarity condition (42) can be rewritten as:

C (r) (ε 0 ) = C (r) ∆L -1 : ∆τ , (46)
In other words, the two stationarity conditions amount to equating the per phase averages and fluctuations of the fields ε 0 and ∆L -1 : ∆τ . (3) When the elastic moduli are uniform per phase and denoted by L (r) , the stationarity conditions take a simpler form: r) .

∆L (r) : ε 0 (r) = -∆τ (r) , ∆L (r) : ε 0 ⊗ ∆L (r) : ε 0 (r) = ∆τ ⊗ ∆τ (
(47) where ∆L (r) = L (r) -L 0 , it is easy to show that the macroscopic stress-strain relation of the original composite, as determined by ( 11), leads to the approximation

σ = L 0 ε + τ 0 , (48)
where the tensors L 0 and τ 0 have to be evaluated at the optimal values of L (r) 0 and τ (r) 0 resulting from expressions ( 41) and (42).

Dual approach

A completely parallel approach can be developed with the dual potential u defined by (4) and the expression of which is given by (5). The dual potential for the linear comparison composite has a similar expression:

u (r) 0 (σ) = sup ε σ : ε -w (r) 0 (ε) = 1 2 σ : M (r) 0 : σ + γ (r) 0 : σ + 1 2 γ (r) 0 : L (r) 0 : γ (r) 0 , (49)
where

M (r) 0 = L (r) 0 -1 , γ (r) 0 = -M (r) 0 : τ (r) 0 . (50)
Now, assuming that (33) holds, the primal effective potential w admits the lower bound (40), which by means of (10), gives an upper bound on the dual effective potential u.

More specifically, for every τ (r) and L (r) :

u(σ) ≤ sup ε {σ : ε -w 0 (ε) -V(x, τ 0 , L 0 ) } ≤ u 0 (σ) -V(x, τ 0 , L 0 ) ,
and therefore

u(σ) ≤ inf τ (r) 0 ,L (r) 0 { u 0 (σ) -V(x, τ 0 , L 0 ) } . ( 51
)
The error function V (r) between w (r) and w (r) 0 , given by ( 35), is also (up to a sign) the error function between u (r) and u (r) 0 . Indeed, we have that

V (r) (x, τ (r) 0 , L (r) 0 ) = inf ε w (r) (x, ε) -w (r) 0 (ε) = inf ε w (r) (x, ε) -sup σ σ : ε -u (r) 0 (x, σ) = inf σ inf ε w (r) (x, ε) -σ : ε + u (r) 0 (x, σ) = inf σ u (r) 0 (x, σ) -u (r) (x, σ) = -sup σ u (r) (x, σ) -u (r) 0 (x, σ) . ( 52
)
The relations (50) can be inverted to express τ

V (r) (x, γ (r) 0 , M (r) 0 ) = 1 2 ∆γ (r) (x) : (∆M (r) (x)) -1 : ∆γ (r) (x) + 1 2 γ (r) 0 : L (r) 0 : γ (r) 0 , ( 53
)
where

∆M (r) (x) = M (x) -M (r) 0 ≤ 0 and ∆γ (r) (x) = γ(x) -γ (r) 0 .
Therefore, the bound (51) can be alternatively written as:

u(σ) ≤ inf γ (r) 0 ,M (r) 0 u 0 (σ) - 1 2 N r=1 c (r) ∆γ (r) (x) : (∆M (r) (x)) -1 : ∆γ (r) (x) (r) , (54)
where the last term in the definition (49) of u (r) 0 should be dropped in the computation of u 0 in the right hand side of this expression, since it cancels out exactly with the second term in expression (53) for the functions V (r) .

It should be noted that the lower bound (54) for u could have been obtained more directly by following a procedure completely analogous to the one outlined earlier for the upper bound (44) for the primal effective energy w under the hypothesis (43). Also, in complete analogy to what was stated in remark 4 of the previous subsection, it can be shown that the macroscopic constitutive relation for the original composite may be expressed in the form

ε = M 0 σ + γ 0 , (55)
where the tensors M 0 and γ 0 have to be evaluated at the optimal values of M (r) 0 and γ (r) 0

resulting from stationarity conditions associated with the optimization problem defined by the right hand side of expression (54).

No strain fluctuation in one of the phases of the LCC

In certain special cases, it can happen that the strain field ε 0 in one of the phases of the LCC is uniform, which implies that its fluctuations C (r) (ε 0 ) vanish. This is the case in all phases of a laminate. This is also the case in one of the phases of composites for which one of the Hashin-Shtrikman bounds is an accurate approximation. When this happens, the stationarity equation ( 46) requires that the LCC be chosen in such a way that the fluctuations of ∆L -1 : τ vanish in the phase where ε 0 has no fluctuation. In other words, one would have to choose L (r) 0

in such a way that ∆L -1 : τ is uniform in this phase which is impossible in general. Consider the case where the elastic moduli are piecewise uniform. Then the condition (46) can be rewritten as

C (r) (∆L : ε 0 ) = C (r) (∆τ ) . ( 56
)
When the fluctuations of ε 0 vanish in one phase, the left-hand-side of (56) vanishes, whereas

C (r) (∆τ ) = C (r) (τ ) + τ (r) -τ (r) 0 ⊗ τ (r) -τ (r) 0
, which cannot vanish except when τ is uniform in this phase. The stationarity equation ( 46) has no solution in this particular phase. This means that the extremum in ( 40) or ( 44) is attained on the boundary of the admissible domain for L (r) 0 , which in turn suggests identifying the elastic moduli in this phase with the original piecewise constant elastic moduli. As will be shown next, this choice combined with other choices in the other phases, leads to an upper bound for the effective energy of the composite.

For definiteness, assume that the strain field in the LCC has fluctuations in all phases except in phase 1. Set: L

(1)

0 = L (1) , τ (1)
0 = τ (1) . (57)
The other L (r) 0 (r = 1) are searched such that ∆L (r) ≤ 0. With the choice (57), one has in phase 1, w

0 (ε) = w (1) (x, ε) + (τ (1)
0τ (x)) : ε, whereas in the other phases (r = 1)

w (r) 0 (x, ε) ≥ w (r) (x, ε) -V (r) (x, τ (r) 0 , L (r) 0 ).
Therefore, denoting by ε 0 the strain field in the LCC:

w 0 (ε) = w 0 (ε 0 ) = N r=1 c (r) w (r) 0 (ε 0 ) (r) ≥ c (1) w (1) (x, ε 0 ) (1) -(τ (1) 0 -τ ) : ε 0 (1) + N r=2 c (r) w (r) (ε 0 ) (r) -V (r) (x, τ (r) 0 , L (r) 0 ) (r)
But since ε 0 has no fluctuation in phase 1 and since τ

(1) 0 is the average of τ over phase 1, the above result reduces to

w 0 (ε) ≥ w(ε 0 ) - N r=2 c (r) V (r) (x, τ (r) 0 , L (r) 0 ) (r) .
Finally, we obtain

w(ε) ≤ w(ε 0 ) ≤ inf τ (r) 0 ,L (r) 0 w 0 (ε) + N r=2 c (r) V (r) (x, τ (r) 0 , L (r) 0 ) (r) ,
which is an upper bound for the effective potential w under the hypothesis that ∆L (r) ≤ 0 (r = 1).

Statistics of the fields

The aim of this section is to show that the first and second moments of the strain and stress fields in the actual non-piecewise uniform composite can be approximated at a good degree of accuracy by the first and second moments of the fields in the LCC.

Following [START_REF] Idiart | Field statistics in nonlinear composites[END_REF], the first and second moments of the strain and stress field in the actual non piecewise-uniform composite can be obtained by perturbing the initial potentials by terms which are piecewise-uniform in each phase. More specifically, consider the following perturbed potentials:

w t (x, ε) = N r=1 χ (r) (x)w (r) t (x, ε), w (r) t (x, ε) = w (r) (x, ε) + t (r) : ε, (58)
and

w λ (x, ε) = N r=1 χ (r) (x)w (r) λ (x, ε), w (r) λ (x, ε) = w (r) (x, ε) + 1 2 ε : λ (r) : ε, (59)
where the t (r) 's and the λ (r) 's are constant second-order and fourth-order tensors respectively. Then, following [START_REF] Idiart | Field statistics in nonlinear composites[END_REF]:

ε (r) = 1 c (r) ∂ w t (ε) ∂t (r) | t=0 , ε ⊗ ε (r) = 2 c (r) ∂ ∂λ (r) w λ (ε)| λ=0 . ( 60
)
Similar relations apply to the first and second moments of the strain in the linear comparison composite.

Using the variational estimate (40) for the effective potential w amounts to making the following approximation:

w(ε) sup τ (r) 0 ,L (r) 0 { w 0 (ε) + V(x, τ 0 , L 0 ) } .
Note that the perturbation by a piecewise uniform polarization t (r) or piecewise uniform moduli λ (r) , does not affect the error function V (r) . In other words, from the inequality ( 36)

w (r) (x, ε) ≥ w (r) 0 (ε) + V (r) (x, τ (r) 0 , L (r) 0 ), one obtains w (r) t (x, ε) ≥ w (r) 0,t (ε) + V (r) (x, τ (r) 0 , L (r)
0 ), with the same function V (r) as above. Consequently, the lower bound (40) still applies to the perturbed potential and can be taken as an approximation for the effective potential:

w t (ε) sup τ (r) 0 ,L (r) 0 { w 0,t (ε) + V(x, τ 0 , L 0 ) } . ( 61
)
Taking the derivative with respect to t (r) of both sides of this equality, and noting that V does not depend on t, we obtain that

ε (r) ε 0 (r) . ( 62
)
The same procedure applied to w λ and w 0,λ shows that

ε ⊗ ε (r) ε 0 ⊗ ε 0 (r) . ( 63
)
As for the moments of the stress field, they are obtained by suitable perturbations of the dual potential u. These perturbations read as:

u η (x, σ) = N r=1 χ (r) (x)u (r) η (x, σ), u (r) η (x, σ) = u (r) (x, σ) + η (r) : σ, (64)
and

u µ (x, σ) = N r=1 χ (r) (x)u (r) µ (x, σ), u (r) µ (x, σ) = u (r) (x, σ) + 1 2 σ : µ (r) : σ. ( 65
)
Then,

σ (r) = 1 c (r) ∂ u η ∂η (r) (σ)| η=0 , σ ⊗ σ (r) = 2 c (r) ∂ u µ ∂µ (r) (σ)| µ=0 . ( 66
)
It is often convenient to compute the above moments of the stress for a given macroscopic strain ε, by means of the primal potential w. Following [START_REF] Idiart | Field statistics in nonlinear composites[END_REF], one can introduce the Legendre transforms w η and w µ of u η and u µ and note (as in proposition 3.4 and corollary 3.7 of [START_REF] Idiart | Field statistics in nonlinear composites[END_REF]) that:

σ (r) = - 1 c (r) ∂ w η ∂η (r) (ε)| η=0 , σ ⊗ σ (r) = - 2 c (r) ∂ w µ ∂µ (r) (σ)| µ=0 . ( 67
)
Similar relations hold for the first and second moment of the stress field in the comparison composite. It should be noted that the error function V (r) between w (r) η and w (r) 0,η is, up to a sign, the same as the error function between u (r) η and u (r) 0,η (see ( 52)) and, since these potentials are obtained from u (r) and u (r) 0 by the same perturbation term, this error function is the original function V (r) between w (r) and w (r) 0 . This leads to variational estimates for the perturbed potentials:

w η (ε) sup τ (r) 0 ,L (r) 0 { w 0,η (ε) + V(x, τ 0 , L 0 ) } , ( 68
)
with a similar relation for w µ . Taking the derivative of the above relation with respect to η (r) (which does not appear in V ) and of the similar relation for w µ with respect to µ, gives:

σ (r) σ (r) 0 , σ ⊗ σ (r) σ 0 ⊗ σ 0 (r) . ( 69
)
In conclusion, the variational estimates for the first and second moments of the local strain and stress fields coincide with those in the LCC.

Application to rigidly reinforced composites

With the objective of illustrating the general results of the previous sections, we provide in this section an application to rigidly reinforced composites with a prescribed distribution of non-uniform polarization field in the matrix phase. For simplicity, the microstructure will be taken to consist of aligned, cylindrical fibers of the rigid phase distributed randomly in a uniform elastic matrix material. Assuming that the modulus tensor of the matrix phase is invariant under a reflection about the fiber direction, the three-dimensional elasticity problem is known [START_REF] Hill | Theory of mechanical properties of fiber-strengthened materials : I. Elastic behavior[END_REF]) to reduce to uncoupled anti-plane strain and generalized plane strain problems. Motivated by this result, we restrict our attention here to the two-dimensional problem defined by plane strain loading of the fiber-reinforced composite in the transverse plane. The matrix elasticity tensor will be further taken to be incompressible and to have the special orthogonal anisotropy (see Ponte Castañeda (2002b)) defined by

L (1) = 2λE + 2µF , (70)
where E and F are orthogonal projection operators, such that E :

F = F : E = 0, E : E = E, F : F = F , E + F = K,
and K is the identity in the space of incompressible, fourth-order tensors in two dimensions. Note that a general representation for E may be given in terms of deviatoric, second-order tensors ˆ with ˆ : ˆ = 1, such that E = ˆ ⊗ ˆ . In addition, it is assumed that the matrix is subjected to a distribution of polarization τ (1) (x) with prescribed phase average τ (1) , and fluctuations C (1) (τ ).

The computation of the various variational estimates requires a 'comparison' composite with the same microstructure as the above-described composite, where the matrix phase now has a modulus tensor L

(1)

0 = 2λ 0 E + 2µ 0 F , (71)
and a constant polarization τ

(1) 0 , while the fiber phase is still rigid. As discussed in remark 4 of section 2, the effective energy of a two-phase composite with constant-per-phase moduli and polarization is completely determined (via expressions ( 28) and ( 29)) by the effective modulus tensor L 0 of the corresponding two-phase elastic composite. However, for the case of rigid fibers, the results simplify dramatically, and we obtain

τ 0 = τ (1) 0 , and f 0 = 0. ( 72
)
Further assuming that the fibers have circular cross-section and are distributed isotropically in the transverse plane, the following estimate of Hashin-Shtrikman type is available (see Ponte Castañeda & Willis (1995)) for the elasticity tensor of the comparison composite

L 0 = L (1) 0 + c 1 -c P (1) 0 -1 , ( 73
)
where c is the volume fraction of the fibers, and the Eshelby-type, microstructural tensor P

(1) 0 is determined by

P (1) 0 = 1 4π |ξ| = 1 H (1) 0 (ξ) dS(ξ), (74)
where H

(1)

0 ijkh (ξ) = N (1) 0 ik ξ j ξ h | (ij)(kh) , with N (1) 0 = K (1) 0
-1 , and K

(1)

0 ik = L (1)
0 ijkh ξ j ξ h . The expression ( 73) is in general a lower bound for the effective stiffness of all isotropic rigidlyreinforced composites. In addition it is known to be attained for specific microstructures and it will assumed in the sequel that the composite under consideration has such a microstructure. Then, making use of expression (71), it follows (Ponte Castañeda (2002b)) that

P (1) 0 = 1 2(λ 0 + √ µ 0 λ 0 ) E + 1 2(µ 0 + √ µ 0 λ 0 ) F , (75)
and also that

L 0 = 2λ 0 1 -c 1 + c µ 0 λ 0 E + 2µ 0 1 -c 1 + c λ 0 µ 0 F . (76)
In conclusion, assuming that the applied strain ε is co-axial with the material anisotropy tensor ˆ , and defining the equivalent strain via ε 2 e = (2/3)ε : ε, the energy function of the comparison composite is given by

w 0 (ε) = 3 2 λ 0 1 -c 1 + c µ 0 λ 0 ε 2 e + τ (1) 0 : ε, (77)
which is seen to be a function of λ 0 , µ 0 and τ

(1) 0 . In addition, the phase average and fluctuations of the strain may be easily computed (see [START_REF] Idiart | Macroscopic behavior and field fluctuations in viscoplastic composites: second-order estimates vs full-field simulations[END_REF]), and are given by

ε 0 (1) = 1 1 -c ε 0 , and C (1) (ε) = 3 4 c (1 -c) 2 µ 0 λ 0 E + λ 0 µ 0 F ε 2 e . (78)
The result (77) may be interpreted as a lower bound for fiber-reinforced composites with the special classes of microstructures described earlier. Alternatively, it can be viewed as an exact estimate for those members of the above class that attain this bound, since it is known that the bound is attainable [START_REF] Milton | The Theory of Composites[END_REF]). In the discussion below, we take this second point of view and make use of the estimate (77), together with the theory developed in section 3, to obtain 'upper' and 'lower' bounds for the corresponding effective energy of the fiber-reinforced composite with a non-uniform distribution of polarization in the matrix phase.

Thus, assuming that λ 0 ≤ λ and µ 0 ≤ µ (so that L

(1) 0 ≤ L (1) ), it follows from expression (40) that

w(ε) ≥ sup τ (1) 0 ,L (1) 0 w 0 (ε) - (1 -c) 2 τ (1) -τ (1) 0 : L (1) -L (1) 0 -1 : τ (1) -τ (1) 0 (1)
.

(79)

Then, using expressions (70) and (71) for L (1) and L

(1) 0 , respectively, and the pertinent stationarity conditions (41) and (42) (or ( 46)), we arrive at the following equations for τ

(1) 0 , λ 0 and µ 0 , namely, τ

(1)

0 = τ (1) + 2(λ -λ 0 ) 1 -c ε, (80)
which is a 'lower' bound for the effective energy when the branch of (85) with λ 0 ≤ µ is used. On the other hand, an 'upper' bound is provided by the same expression when the λ 0 ≥ µ branch is used instead. In addition, the expression for the macroscopic stress (83) reduces to

σ e = 3µ 1 -c   1 + c λ 0 µ 1/2   ε e , (87)
where again the λ 0 ≤ µ and λ 0 ≥ µ branches of (85) must be used for the stress associated with the 'lower' and 'upper' bounds on the energy (not necessarily bound on the stress). Finally, it is pointed out that estimates for the average and fluctuations of the strain in the matrix phase may be obtained via expressions (78) for the comparison composite with µ 0 = µ and the appropriate choice of λ 0 .

For comparison purposes, we also include here the corresponding specialization for the Reuss lower bound for this special case of a rigidly reinforced two-dimensional composite, as given by

w Reuss (ε) = 3 2 µ 1 -c ε 2 e - 1 -c 4µ τ 2 . ( 88
)
Note that this bound grows quadratically with the applied strain ε e -and decreases quadratically with the polarization fluctuations τ . This behavior should be contrasted with the corresponding behavior of the 'lower' and 'upper' bounds, which as discussed in more detail in the next paragraph, is non-quadratic in both the applied strain and polarization fluctuations. In addition, we provide the corresponding specialization of the effective energy making use of the 'first moment' approximation ( 2)

w(ε) = 3µ 2 1 + c 1 -c ε 2 e , (89)
which is seen to be independent of the polarization fluctuations in view of the fact that τ (1) = 0 in this particular example. Note that this estimate lies between the bounds, and agrees exactly with them when the polarization fluctuations vanish.

The results for the above-described composite are shown in Fig. 1 for a fiber volume fraction of 25%. Part (a) shows the value of λ 0 , normalized by µ, for the 'lower' and 'upper' bounds, as functions of the dimensionless parameter τ / (µε e ). It can be seen that λ 0 → µ not only when the polarization fluctuations τ vanish, but also when the applied strain ε e increases. On the other hand, λ 0 tends to its extreme values of zero and infinity, when either the fluctuations increase or the applied strain decreases. Parts (b) and (c) show the new 'lower' and 'upper' bounds, as determined by ( 86), as well as the Reuss lower bound (88) and 'first moment' estimate (89) for the effective energy, normalized with respect to 3τ 2 /(2µ) and (3/2)µε 2 e , as functions of µε e /τ and its inverse τ / (µε e ), respectively. The first of these plots illustrates the behavior of the effective energy as a function of the applied strain for fixed polarization fluctuations, while the second shows the corresponding behavior of the effective energy as a function of the polarization fluctuations, for fixed applied strain. Thus, it can be seen that while the new 'lower' bound improves on the Reuss lower bound only for sufficiently large strains, the opposite behavior is observed as a function of the polarization fluctuations. In fact, it can be shown that the new 'lower' bound tends to the Reuss lower bound as µε e /τ → 0, so that both lower bounds on the effective energy are seen to tend to -∞ quadratically with the polarization fluctuations, in agreement with the comments provided in remark 2 of section 2. On the other hand, for the microstructures for which the HS estimate ( 73) is attained, the new 'upper' bound is significantly tighter than the classical Voigt bound (which corresponds to a much larger class of microstructures and happens to be infinite in this case). Thus, for fixed polarizations, the range of possible behaviors-as determined by the bounds on the effective energy of the composite-remains bounded as the strain increases, although the differences continue to increase in magnitude as the polarization fluctuations increase. We also observe that the 'first moment' estimate satisfies the bounds, but as already mentioned, it is completely insensitive to the polarization fluctuations. Finally, part (d) show plots of the macroscopic stress, normalized by the applied strain, as a function of τ / (µε e ). Although the results for the 'upper' and 'lower' estimates are not necessarily bounds, it can be seen that the stress (unlike the effective energy) associated with the 'lower' bound does not decrease without a bound with increasing polarization fluctuations, tending to the Reuss lower estimate as the polarization fluctuations increase. This makes it clear that the quadratic drop in the effective energy with the polarization fluctuations arises from the polarization terms in the energy and not from the strain terms.

Concluding remarks

In this work, variational approximations have been developed to estimate the macroscopic response, as well as the phase averages and covariance of the fluctuations of the stress and strain fields in elastic composite materials with non-uniform distributions of polarization and/or moduli. The approximation consists in the replacement of the actual composite by a 'comparison' composite with uniform, per-phase polarization and moduli that are determined by a suitably designed variational statement. When the moduli of the actual composite are uniform per phase, the resulting estimates for the macroscopic behavior are found to depend on the first and second moments of the polarization distributions in each phase. On the other hand, it has been shown that the first and second moments of the stress and strain fields in the actual composite may be estimated consistently from the corresponding first and second moments of the stress and strain fields in the comparison composite. Furthermore, when the modulus tensor of the comparison composite is chosen to be greater (lower) than the modulus of the actual composite (in the sense of quadratic forms), the resulting estimate for the effective free energy of the composite is shown to be an upper (lower) bound. Although this was not explored in this work, it is also possible to make other choices for the phase moduli of the comparison composite (see Ponte Castañeda (2002a) and Lahellec & Suquet (2007b)), also leading to stationary variational estimates, but which may not necessarily be bounds. Such possibilities are application-dependent and will continue to be explored in future work.

Appendix: Unidimensional case

In dimension 1, the constitutive relations read as σ(x) = (x) : ε(x) + τ (x) where both and τ may have fluctuations within a single phase. The equilibrium equations, the compatibility conditions and the boundary conditions take the form dσ dx = 0, ε(x) = du dx (x), u(0) = 0, u(L) = εL.

It follows from the equilibrium equation that the Reuss model is exact σ(x) = σ. Therefore:

w(ε) = 1 2 1 m ε 2 + mτ m ε + 1 2 mτ 2 m - 1 2 mτ 2 and σ = 1 m ε + mτ m .
In the LCC with piecewise uniform moduli and polarization fields, there is no strain fluctuation in the phases. Therefore the variational procedure is not applicable in this case and a direct approach is preferable. A piecewise uniform LCC will give the exact energy and the exact overall stress iff This system has two solutions:

τ (1) 0 = mτ m ± c (2) m (2) c (1) m (1) mτ 2 m - mτ 2 m 2 , τ (2)
0 = mτ m -± c (1) m (1) c (2) m (2) mτ 2 m - mτ 2 m 2 ,
where the ± signs are the same in both expressions. These solutions differ from τ (1) and τ (2) respectively.

  The right-hand side of (38) can be considered as an approximation of the actual potential w. Taking advantage of the stationarity of the variational lower bound (40) (or the corresponding upper bound (44)) with respect to the variables L

  (r) can be considered as a function of the strain polarization γ fact, it is easy to show that

  Fig. 1. (a) Solutions for λ 0 of 'comparison' composite for 'lower' and 'upper' bounds, as function of applied strain ε e and parallel polarization fluctuation τ . (b) Comparison of the new 'lower' and 'upper' bound for the effective energy with the Reuss lower bound and 'first moment' estimate, plotted as functions of the applied strain (for fixed polarization fluctuation). (c) Comparison of the new 'lower' and 'upper' bound for the effective energy with the Reuss lower bound and 'first moment' estimate, plotted as functions of the polarization fluctuation (for fixed applied strain). (d) Comparison of the new 'lower' and 'upper' estimates with the Reuss and 'first moment' estimates for the normalized stress, as functions of the polarization fluctuation (for fixed applied strain).

  us give one of them for which the elastic moduli in the LCC coincide with those in the actual composite (where the elastic moduli are assumed to be piecewise uniform). The LCC is searched such that m (r) 0 = m (r) . (90) gives: c (1) m (1) τ

Acknowledgements

This work was partially supported by the Centre National de la Recherche Scientifique

3c

(1c) 2 µ 0 λ 0 (λλ 0 ) 2 ε 2 e = τ (1)τ (1) : E : τ (1)τ (1) (1) , (81) 3c (1c) 2 λ 0 µ 0 (µµ 0 ) 2 ε 2 e = τ (1)τ (1) : F : τ (1)τ (1) (1) .

(82)

When solutions are found to the system of the two equations ( 81) and ( 82) for λ 0 ≤ λ and µ 0 ≤ µ (note that τ

(1) 0 is determined in terms of λ 0 from (80)), and the results are used in the right-hand-side of (79), a 'lower' bound is obtained for the effective energy w of the fiber-reinforced composite with a non-uniform distribution of polarization in the matrix phase. On the other hand, according to what was stated in remark 1 of section 3, when solutions of the system (81) and ( 82) with λ 0 ≥ λ and µ 0 ≥ µ are obtained, the expression in the right-hand-side of ( 79) provides an 'upper' bound for the effective energy w of the fiber-reinforced composite with a non-uniform distribution of polarization in the matrix phase. In addition, the corresponding estimates for the macroscopic stress may be obtained from expression (48) such that

where we have already made use of expression (80) for τ

(1) 0 to simplify the result. Thus, this expression provides a 'lower' estimate (not necessarily a bound) for the stress when solutions such that λ 0 ≤ λ and µ 0 ≤ µ are obtained from expressions ( 81) and ( 82), as well as an 'upper' estimate (again not necessarily a bound) for the stress when solutions such that λ 0 ≥ λ and µ 0 ≥ µ are obtained from expressions ( 81) and ( 82).

Given that our interest here is to explore the effect of fluctuations in the distribution of the polarization field, we will assume in the illustrative example below that the matrix phase is isotropic, so that λ = µ, and that the polarization field has zero average, so that τ (1) = 0, and is co-axial with the direction defined by the applied field ε, so that τ (1) (x) = τ (1) (x)ˆ . Then, the 'orthogonal' component (F projection) of the polarization fluctuations vanishes, and defining the 'parallel' component of the polarization via τ = τ : ˆ , the parallel fluctuations are described by

Under these further assumptions, it follows from (82) that µ 0 = µ, and from (81) that the value of λ 0 ≤ µ for the 'lower' bound is determined from the solution of the equation

The corresponding value of λ 0 ≥ µ for the 'upper' bound is determined by replacing the factor 1λ 0 /µ by λ 0 /µ -1 in the same equation. Then, the expression (79) for the lower bound on the effective energy of the composite with the non-uniform distribution of polarization specializes to


Michel Planat

Fabio Anselmi

Patrick Solé

Pauli graphs, Riemann hypothesis, Goldbach pairs

Keywords: numbers: 02.10.De, 02.10.0x, 03.65.Fd, 03.67.Lx Mathematics Subject Classification: 11M26, 11P32, 05C10, 81P68, 11A25

. In the extreme high dimensional case, at primorial numbers N r , it is known that the inequality ψ(Nr) Nr log log Nr e γ ζ(2) (for every r > 2) is equivalent to Riemann hypothesis. Introducing the Hardy-Littlewood function R(q) = 2C 2 p|n p-1 p-2 (with C 2 ∼ 0.660 the twin prime constant), that is used for estimating the number g(q) ∼ R(q) q ln 2 q of Goldbach pairs, one shows that the new inequality R(Nr) log log Nr e γ is also equivalent to Riemann hypothesis. In this paper, these number theoretical properties are discusssed in the context of the qudit commutation structure.

Introduction

We propose new connections between the Pauli graphs [START_REF] Planat | On the Pauli graphs on N -qudits[END_REF][START_REF] Planat | Pauli graphs when the Hilbert space dimension contains a square: why the Dedekind psi function?[END_REF], that encode the commutation relations of qudit observables, and prime number theory. We already emphasized that the Dedekind psi function ψ(q) = q p|q (1 + 1 p ) (with p a prime) is used to count the number of maximal commuting sets of the qudits [START_REF] Planat | Pauli graphs when the Hilbert space dimension contains a square: why the Dedekind psi function?[END_REF] and meets the Riemann hypothesis (RH) at primorial numbers q ≡ N r = 2 • • • p r [START_REF] Solé | Extreme values of the Dedekind ψ function[END_REF]. Similarly, there exist striking connections between ψ(q) and the Hardy-Littlewood function g(q) = R(q) q log 2 (q)

for the Goldbach distribution of prime pairs [see Sec. 3 for the definition of R(q)]. In particular, we observe that ψ(q) meets the so-called totally Goldbach numbers at small q's and that R(q) < ζ(2) ψ(q) q also meets RH at primorial numbers.

The Euler constant γ = lim n→∞ ( n k=1 1 klog n) ∼ 0.57721, through the Mertens formula e γ = lim n→∞

1 log n p≤n (1-1 p ) -1
, is an important ingredient of all the inequalities involved in this correspondence.

In the rest of this section, we report on the number theoretical coincidence between ψ(q) and the totally Goldbach numbers at small q's, as well on the already known theorem connecting ψ(q) and RH at primorial numbers N r . In Section 2, we explore in detail this coincidence by refering to the qudit Pauli graphs. Then, in Section 3, we establish the connection between R(q) and RH at primorial numbers. In the discussion, the concept of a Goldbach defect for encompassing the statements at low and high q's is proposed.

A number theoretical coincidence

Let us start the exposition of our ideas with a few theorems and definitions.

Goldbach's conjecture, formulated in 1742, is that every even integer greater than 2 is the sum of two primes. To date it has been checked for q up to 2.10 18 [START_REF] Oliveira E Silva | Goldbach conjecture verification[END_REF]. A pair (p 1 , p 2 ) of primes such that the even integer n = p 1 + p 2 is called a Golbach partition.

Definition 1: A positive integer n is totally Goldbach if for all primes p < n -1, with p not dividing n (except when p = np) we have that np is prime [START_REF] Brouwers | Totally Goldbach numbers and related conjectures[END_REF].

Theorem 1: The inequality ψ(q) q > e γ log log q is only satisfied at a totally Goldbach number q ∈ A = {2, 3, 4, 5, 6, 8, 10, 12, 18, 30}. The only totally Goldbach numbers not satisfying the inequality are q = 1 and q = 24.

Proof: The proof of theorem 1 follows from a combination of results in [START_REF] Brouwers | Totally Goldbach numbers and related conjectures[END_REF] and [START_REF] Solé | Extreme values of the Dedekind ψ function[END_REF] (for a more general setting, see [START_REF] Planat | Riemann hypothesis and quantum mechanics[END_REF]).

Definition 2: A positive integer n is almost totally Goldbach of index r if for all primes p < n -1, with p not dividing n (except when p = np) we have that np is prime with r exceptions. Let g(n) be the number of ways of representing the integer n as the sum of two primes. The maximum value of g(n) is indeed less or equal than the number of primes

n 2 ≤ p ≤ n -1.
Values of n such that g(n) reaches its maximum are in the set B = A 0 ∪ {7, 14, 16, 36, 42, 48, 60, 90, 210} [START_REF] Deshouillers | An upper bound in Golbach's problem[END_REF], where A 0 is the set of totally Goldbach numbers. It is not surprising that numbers in B that are not totally Goldbach are almost totally Goldbach with a small index r (as shown in Table 1). The first five and the integer 60 have index 1, while 48 and 90 have index 2 and 210 has index 4. For a prime number p > 3, the index r(p) is the Sloane's sequence A062302. This number theoretical coincidence is further explored in Sec. 2 in the context of the qudit commutation structure.

Number theoretical inequalities at primorial numbers

Let us start with

Theorem 2: Let N r = 2 • • • p r be the primorial of order r. The statement

ψ(Nr )
Nr log log Nr > e γ ζ(2) for every r > 2 is equivalent to the Riemann hypothesis. Proof: Theorem 2 is proved in [START_REF] Solé | Extreme values of the Dedekind ψ function[END_REF] and in a more general form in [START_REF] Planat | Riemann hypothesis and quantum mechanics[END_REF]. In Sec. 3, we establish further statements relating ψ(q) and the Hardy-Littlewood function R(q) about the number of Goldbach pairs. An important result is conjecture 1 that is found to be equivalent to RH by theorem 5.

The number theoretical coincidence for qudits

In this section, we recall the definition of a q-level system, or qudit, and the number theoretical structure, involving the Dedekind psi function ψ(q), of maximal sets of pairwise commuting operators. Calculations about the inequality of theorem 1, at the number q, and the symmetries of the qudit system in the corresponding dimension q are given in two separate tables 2 and 3, in order to display the distinctive features of both types of systems. In table 2, one selects q ∈ A 0 , a totally Goldbach number, and in table 3 one selects q / ∈ A 0 , with q ≤ 36.

The single qudit commutation structure Definition 3: In a loose sense, a qudit is defined as the group P q = X, Z with unitary quantum generators

X =        0 0 . . . 0 1 1 0 . . . 0 0 . . . . . . . . . . . . . . 0 0 . . . 1 0        , Z = diag(1, ω, ω 2 , . . . , ω q-1 ), (1)
with ω = exp( 2iπ q ). The Weyl pair (X, Z) satisfies ZX -ωXZ = 0, and each element of P q may be uniquely written as ω a X b Y c , with a, b, c in the modular ring Z q . One has |P q | = q 3 , P ′ q = Z(P q ) (the derived subgroup equals the center of the group, and is of order q). One is interested in the maximal mutually commuting sets within the central quotient Pq = P q /Z(P q ), with the identity I q removed. Theorem 3: Let S be the set of maximal mutually commuting sets M in Pq . The cardinality of M is |M| = q -1 and that of S is |S| = σ(q), where σ(q) is the sum of divisor function.

Proof: The proof of theorem 3 is given in [START_REF] Albouy | The isotropic lines of Z 2 d[END_REF] (for the relation to σ(q), see [START_REF] Planat | Pauli graphs when the Hilbert space dimension contains a square: why the Dedekind psi function?[END_REF]). Hint: Using the Weyl pair property one writes the group theoretical commutator as [X, Z] = XZX -1 Z -1 = ω -1 I q so that one gets the expression

ω a X b Z c , ω a ′ X b ′ Z c ′ = ω cb ′ -c ′ b I q , meaning that two elements of P q commute if only if the determinant ∆ = det b ′ b c ′ c vanishes. Two vectors such that their symplectic inner product [(b ′ , c ′ ).(b, c)] = ∆ = b ′ c -bc ′ vanishes
are called perpendicular. Thus, one can transfer the study of commutation relations within the group P q to the study of perpendicularity of vectors in the ring Z 2 q [9] with the isomorphism

(P q /Z(P q ), ×) ∼ = (Z 2 q , +) (2)
between the central quotient of P q and the algebra of vectors in the Z q -modular lattice Z 2 q , endowed with the symplectic inner product ".". Then, one defines a isotropic line as a set of q points on the lattice Z 2 q such that the symplectic product of any two of them is 0(mod q). To such an isotropic line corresponds a maximal commuting set in P q /Z(P q ).

Taking the prime power decomposition of the Hilbert space dimension as q = i p s i i , it is shown in (18) of [START_REF] Albouy | The isotropic lines of Z 2 d[END_REF] that the number of isotropic lines of the lattice

Z 2 q reads i p s i +1 i -1 p i -1 ≡ σ(q). (3)
Theorem 4: Let S be the set of maximal mutually commuting sets M in Pq and D ⊂ S the subset of S such that for any d ∈ D and s ∈ S the intersection of d and s is non empty, i.e. d ∩ s = ∅. The cardinality of S -D is the Dedekind psi function ψ(q).

Proof: As shown in [START_REF] Albouy | The isotropic lines of Z 2 d[END_REF], a isotropic line of Z 2 q corresponds to a Lagrangian submodule, i.e. a maximal module such that the perpendicular module M ⊥ = M. Let us now specialize to Lagrangian submodules that are free cyclic submodules

Z q (b, c) = {(ub, uc)|u ∈ Z q } ,
for which the application u → (ub, uc) is injective. Not all Lagrangian submodules are free cyclic submodules. A point x = (b, c) such that Z q (b, c) is free is called an admissible vector, and the set of free cyclic submodules defined by the admissible vectors is called the projective line

P 1 (Z q ) = {Z q (b, c)|(b, c) is admissible} .
Following theorem 5 in [START_REF] Havlicek | Projective line of an arbitrary single qudit[END_REF], the number of points of the projective line is

|P 1 (Z q )| = i (p s i i + p s i -1 i ) ≡ ψ(q), (4)
where ψ(q) is the Dedekind psi function (see also [START_REF] Vourdas | Symplectic transformations and quantum tomography in finite quantum systems[END_REF]). Then, using the isomorphism relation (2) between the isotropic lines of Z 2 q and the maximal commuting sets of S, it is clear that a free cyclic submodule is isomorphic to an element of the set S -D, and that the projective line is in bijective correspondance with S -D. This completes the proof.

Pauli graphs over a single qudit system Definition 4: Let G q be the Pauli graph constructed by taking the observables as vertices and an edge joining two mutually commuting observables (with the identity I q removed). Indeed, |G q | = q 2 -1. One further defines a point/line incidence geometry with observables as points and the maximum cliques as lines. One characterizes this geometry by creating a dual graph G ⋆ q such that the vertices are the cliques (there are σ(q) of them) and an edge joins two non-intersecting maximum cliques. The connected component of G ⋆ q corresponds to the graph of the projective line Table 2. Data about the structure maximal mutually commuting sets for qudits when q ∈ A 0 (the set of totally Goldbach numbers). The symbols #tg and #mc denote the number of totally Goldbach pairs and the number of maximal mutually commuting sets, respectively. In column 4, one computes the difference ǫ(q) = ψ(q) qe γ log log q and in column 5 the automorphism group of the corresponding projective line. The graphs or the non trivial subgraphs that may identified at the single or multiple point intersection set are displayed in colum 6. The notations K n , K n1,•••,n l , C n , Gr p,q and Cu are for the complete graph with n vertices, a l-partite graph, the polygon graph with n vertices, the grid graph p × q and the cube graph respectively, respectively. The notation BP 6 means a bipartite graph of spectrum {0 8 , (-6) 1 , 6 1 , (-2) 3 , 2 3 }. Here, the symbol × corresponds to non-intersecting copies of the relevant subgraph. 3. Data about the structure of maximal mutually commuting sets for qudits when q / ∈ A 0 . The two symbols r and #atg denotes the index of a almost totally Goldbach pair (see definition 2) and the number of them. The symbol #mc denotes the number of maximal mutually commuting sets. In column 4, one computes the difference ǫ(q) = ψ(q)

P 1 (Z q ) [of size ψ(q)]. q #tg #mc ǫ(q) aut(G * q ) graph & subgraphs 2 0 3 2.15 S 3 K 3 3 0 4 1.16 S 4 K 4 4 1 6 + 1 0.92 G 48 = Z 2 × S 4 K 2,
qe γ log log q and in column 5 the automorphism group of the corresponding projective line. The graphs or the non trivial subgraphs that may identified at the single or multiple point intersection set are displayed in colum 6. Notations are the same than for table 2. The notation BP r means the bipartite graph of prime index r and spectrum {-r 1 , r 1 , (-1) r , 1 r }; the graph BP 3 is the cube graph Cu. The notation BP 9 means a bipartite graph of spectrum {0 16 , (-9) 1 , 9 1 , (-3) 3 , 3 3 }. The spectrum of the subgraph SG found for q = 36 is {8 1 , 0 9 , (4) 4 , 2 4 }.

Definition 5: Let G be a graph. The set of graph eigenvalues λ i (with multiplicities n i ), also called the graph spectrum of G, is denoted spec(G) = {λ n 1 1 , λ n 2 2 , . . .}. A graph whose spectrum consists entirely of integers is called an integral graph.

Results

Main properties of graphs G q and G * q are displayed in tables 2 and 3. A few noticeable properties are summarized below.

Result 1: Let G ⋆ q the dual Pauli graph of the definition 4, i.e. the graph attached to the projective line P 1 (Z q ), integral graphs of the following type are found type 1: if q = r, with r a prime, spec(G ⋆ q ) = {(r + 1) 1 , (-1) r+1 }. type 2: if q = rs, with r, s two distinct primes, spec(G ⋆ q ) = {q 1 , 1 q , (-r) s , (-s) r }. type 3:

if q = rst, with r, s, t three distinct primes, spec(G ⋆ q ) = {q 1 , t rs , s rt , r st , (-1) q , (-rs) t , (-rt) s , (-st) r }.

type 4: if q = 4r, with r a odd prime, spec(G ⋆ q ) = {q 1 , 2 2r , 0 ψ(q)/2 , (-4) r , (-2r) 2 }. type 5: if q = p l is the power of a prime p, with l > 1, spec(G ⋆ q ) = {q 1 , 0 ψ(q)-(p+1) , (-q/p) p }. type 6: in dimension q ≤ 36, cases not of the above type occur for q = 18 with spec(G ⋆ q ) = {18 1 , 3 6 , 0 24 , (-6) 3 , (-9) 2 }, for q = 24 with spec(G ⋆ q ) = {24 1 , 4 6 , 0 36 , (-8) 3 , (-12) 2 } and for q = 36 with spec(G ⋆ q ) = {36 1 , 6 6 , 0 60 , (-12) 3 , (-18) 2 }.

Result 2: Let G ⋆ q the dual Pauli graph of the definition 4. Graphs and subgraphs of the following type are found in colum 6 of tables 2 and 3.

type 1: if q = r, with r a prime, G ⋆ r = K r+1 . type 2: if q = rs, with r, s two distinct primes, copies of complete subgraphs of the form p j × K p i +1 , i, j ∈ {1, 2}, are found. types 4, 5 and 6: For a dimension q containing a square, the projective line P 1 (Z q ) contains bipartite or multipartite subgraphs. For q ∈ A 0 , they are of the type K 2,2 , K 2,2,2 , Cu or BP 6 (see column 6 in Table 2). For q / ∈ A 0 , bipartite or multipartite subgraphs of the type K 4,4 , K 4,4,4 and BP r , r = 6 are found (see column 6 in Table 3).

for type 3 and higher dimensions of the primorial type one gets result 3. Result 3: For primorial dimensions q = N r = r i=1 p i , r > 3, the automorphism group of the graph G * q is a product of the symmetric groups S p i +1 arising from the factors of q. In addition, according to theorem 2, if Riemann hypothesis is satisfied, the size of the projective line follows as |P 1 (Zq)| q log log q > e γ ζ(2) ∼ 1.08.

Comments

Let us discuss shortly the results collected in tables 2 and 3. At first sight, the symmetries of a qudit in a dimension q ∈ A 0 corresponding to a totally Goldbach number (in table 2) do not contrast with those at small dimensions q / ∈ A 0 , q ≤ 36 (in table 2). Let us restrict our reading of the tables to the cases where q contains a square so that the total number σ(q) of isotropic lines/maximal commuting sets, in (3), is strictly bigger than the size ψ(q) of the projective line, in (4). Hence, our decomposition of the cardinality of the set of maximal commuting sets is as in column 3: σ(q) = ψ(q) + (.) [START_REF] Planat | Pauli graphs when the Hilbert space dimension contains a square: why the Dedekind psi function?[END_REF].

For these non square free cases, in table 2 (column 5) the two groups G 48 of the quartit and G 144 of the sextit are building blocks of the automorphism groups. Looking at the dual Pauli graph G * q as a simplicial complex, both groups correspond to the topology of the sphere and of the torus, respectively. Non trivial grid graphs or bipartite graphs are found as constituents of G * q , as shown in column 6. In table 3, the graph G * q for non square free dimensions shows symmetries of a different topology involving higher order symmetric groups and simpler bipartite subgraphs such as BP r (see the legend of the table for a definition).

Dimension 24 deserves a special attention because it is the only totally Goldbach number that violates the inequality of theorem 1. Note that dimension 24 is also very special for lattices (such as the Leech lattice) and marks a kind of phase transition for symmetries [START_REF] Conway | Sphere packings, lattices and groups third edition[END_REF], of a different type than the one described in [START_REF] Planat | Riemann hypothesis and quantum mechanics[END_REF].

On Goldbach pairs at primorial dimensions and RH

According to Hardy and Littlewood [START_REF] Hardy | Some problems of partition numerorum, III: on the expression of a number as a sum of primes[END_REF], the number of Goldbach pairs (p 1 , p 2 ) of primes whose sum is the even integer q is asymtotically given by the Hardy-Littlewood function g(q) = R(q) q log 2 q , with R(q) = 2C

C 2 = r i=2 (1 - 1 (p i -1) 2 ) ∞ i=r+1 (1 - 1 (p i -1) 2 ),
and one uses the equality p i -1

p i -2 = 1/(1 -1 p i -1 ) to obtain R(N r ) = 2 ∞ i=r+1 (1 - 1 (p i -1) 2 ) r i=2 (1 + 1 p i -1
).

(5)

Then, by the convergence of the infinite product for

C 2 R(N r ) ∼ 2 r i=2 (1 + 1 p i -1 ) = 2 r i=2 (1 - 1 p i ) -1 = r i=1 (1 - 1 p i ) -1 .
Finally, proposition 1 follows by using the same line of reasoning than in proposition 3 of [START_REF] Solé | Extreme values of the Dedekind ψ function[END_REF].

Conjecture 1: The Hardy-Littlewood inequality R(Nr ) log log Nr > e γ holds for every r > 1.

This proposed conjecture follows from numerical evidence as shown on table 4. In theorem 5, conjecture 1 is found to be equivalent to Riemann hypothesis.

Proposition 2: For an even integer q, the following inequality holds

R(q) < ζ(2) ψ(q) q . ( 6
)
Proof: One introduces the ratio x(q) = 1

ζ(2)
R(q) ψ(q)/q that reads.

x(q) = 2 p>2 (1 -1 (p-1) 2 )

4 3 p>2 1 1-1 p 2 p>2,p|q 1 3 2 (1 + 1 p )(1 -1 p-1 )
. 9.1 × 10 45336 0.0010 10 5 1.9 × 10 563920 0.00023 Table 4. Numerical evidence for conjecture 1 (in column 3).

By splitting the product as p∤q () × p|q () the ratio x(q) is rewritten as

x(q) = p>2,p∤q 1 -1 (p-1) 2 1 -1 p 2 p>2,p|q (1 + 1 p -1 )(1 - 1 p ).
The right hand side product p|q () equals 1 and the first one p>2,p∤q is lower than 1 so that Proposition 2 is satisfied.

Corollary 1: The primorial numbers N r are champion numbers, i.e. left to right maxima of the function q → x(q). One gets lim r→∞ x(N r ) = 1.

Proof: According to [START_REF] Montgomery | Multiplicative Number Theory I: Classical Theory[END_REF], p. 90, or [START_REF] Sándor | Handbuch of Number Theory I[END_REF], p. 238, the number of Goldbach pairs is such that there exists a constant c 1 with R(q) < c 1 ψ(q) q . Using the same method than the one used in [START_REF] Solé | Extreme values of the Dedekind ψ function[END_REF] for the champions of ψ(q)/q one gets c 1 ≡ ζ(2).

Theorem 5: Conjecture 1 is equivalent to RH. Proof: If Riemann hypothesis is not true, then according to Theorem 4.2 of [3] and the inequality (6) there are infinitely many r such R(Nr ) log log Nr < e γ , corresponding to as many violations of the Hardy-Littlewood inequality in Conjecture 1. Now we prove that Conjecture 1 follows from Riemann hypothesis. It is known from theorem 3 (b) of [START_REF] Nicolas | Petites valeurs de la fonction d'Euler[END_REF] that if RH is true then, for any x ≥ 2,

f (x) = e γ log θ(x) p≤x (1 - 1 p ) < 1,
with θ(x) = p≤x log p the first Tchebycheff function.

We introduce the new function

g(x) = e γ 2C 2 log θ(x) p≤x p -2 p -1
and show below that, under RH, we have

g(x) < 1. (7)
Let us explicit the ratio g

(x) f (x) = 1 C 2 2<p≤x p(p -2) (p -1) 2 ,
where

C 2 = 2<p≤x (1 -1 (p-1) 2 ) p>x (1 -1 (p-1) 2 ). One immediately gets log g(x) = log f (x) - p>x log(1 - 1 (p -1) 2 ). (8)
Next, using the change of variables y = x -1, one is interested with the range of values such that p -1 > y so that log(1 -1 (p-1) 2 ) < log(1 -1 y 2 ). Replacing the sum by an integral at the right hand side of ( 8), one obtains

- p>x log(1- 1 (p -1) 2 ) < - ∞ y log(1- 1 u 2 )du = -u log(1 - 1 u 2 ) + log( u + 1 u -1 ) ∞ y ∼ 1 y ∼ 1 x ,
at large y. By [START_REF] Nicolas | Petites valeurs de la fonction d'Euler[END_REF]Theorem 3 (b)] this is dominated by log f. Hence, by (8) log g < 0 and the inequality (7) follows. Then, by using the same reasoning that in [START_REF] Nicolas | Petites valeurs de la fonction d'Euler[END_REF], Conjecture 1 follows from the inequality g(x) < 1 by observing that g(p r ) = e γ log log N r R(N r ) .

This completes the proof of Theorem 5.

Discussion

It is useful to encapsulate some of the results obtained for Goldbach pairs at small even dimensions and at primorial dimensions by introducing a measure that we denote the Goldbach defect. Definition 6: For a positive even integer q, one defines a Goldbach defect gd q from the relation gd q = ζ(2)qR(q)ψ(q).

Result 4: The champions of the function gd q (the numbers at which gd q reaches a new record) form the sequence C = {2, [START_REF] Oliveira E Silva | Goldbach conjecture verification[END_REF][START_REF] Planat | Riemann hypothesis and quantum mechanics[END_REF][START_REF] Hardy | Some problems of partition numerorum, III: on the expression of a number as a sum of primes[END_REF]18,24,30,42,54,60,84,90,120,150,180,210,270,300,330,390,420,510,570,630,780,840,990,1050,1260,1470,1650,1680,1890,2100,2310,2730,3150,3360,3570,3990,4290,4410,4620,5250,5460,6090,6510,6930,7770,7980,8190,9030,9240 . . .}

The sequence C is found to be a subset of the champions of the cototient function (A051953 in Sloane's encyclopedia). By its definition, sequence C encompasses the even totally Goldbach numbers, except for the number 10 (see definition 1) and the primorial numbers (see Conjecture 1 and Corollary 1). None of the (small) even dimensions in table 3 belongs to the sequence. INumbers in C satisfy ζ(2)ψ(i) i log log i > R(i) log log i > e γ for all i ∈ C, i = 2.

In future work, the Goldbach defect may turn to be useful for discussing the change in the symmetries of the Pauli graphs. The first not square free and not totally Golbach dimension in the sequence C is q = 60. The cardinality of maximal commuting sets for this case decomposes as σ(60) = ψ(60) + (.), i.e. 168 = 144 + 24.

Table 1 .

 1 Almost

totally Goldbach numbers of index r ≤ 10.

  Let N r = 2 • • • p r be the primorial number of order r, we have lim r→∞ R(Nr ) log log Nr = e γ . Proof: One splits the infinite product for C 2 as

	2 p|q	p-1 p-2 , and C 2 = p>2 (1 -1 (p-1) 2 ) ∼ 0.660 is the
	twin prime constant.
	Proposition 1: