B E Peigney

M Peigney

Bounds for nonlinear composite conductors via the translation method

Keywords: composite materials, nonlinear homogenization, conductivity, bounds, translation method

Hashin-Shtrikman type bounds are proposed for nonlinear isotropic composite conductors in two dimensions. Those bounds are obtained by combining the translation method with the idea of embedding the original two-dimensional problem in an extended problem of dimension 6. Invariance properties allow the evaluation of the bounds to be dramatically simplified. Explicit results are obtained for the problem of dielectric breakdown. Numerical results are given for two-phase composites governed by power-law energy functions. The obtained bounds are shown to improve on the linear comparison bounds of the Hashin-Shtrikman type that are delivered by the Talbot-Willis (1985) approach and the Ponte Castañeda (1991) variational method.

Introduction

This paper is concerned with bounding the effective energy of a nonlinear composite conductor in terms of the properties of its constitutive phases. A lot of work has been devoted to that topic, especially for linear composites (see e.g. the monograph by [START_REF] Milton | The Theory of Composites[END_REF]). [START_REF] Wiener | Die theorie des mischkörpers für das feld der stationären strömung. 1. abhandlung: Die mittelwertsätze für kraft, polarisation und energie[END_REF] proposed general bounds on the effective conductivity that hold for any microstructure. For statically isotropic microstructures, tighter bounds have been obtained by [START_REF] Hashin | A variational approach to the theory of the effective magnetic permeability of multiphase materials[END_REF] using a variational approach that employs a homogeneous 'linear comparison' medium. The Hashin-Shtrikman bounds are known to be optimal, i.e. they are achieved by certain isotropic microstructures. The variational approach of [START_REF] Hashin | A variational approach to the theory of the effective magnetic permeability of multiphase materials[END_REF] has later been extended by [START_REF] Willis | Bounds and self-consistent estimates for the overall properties of anisotropic composites[END_REF] and further generalized to nonlinear composites by [START_REF] Talbot | Variational principles for inhomogeneous non-linear media[END_REF]. A related approach for bounding the effective energy of nonlinear composites has been advanced by [START_REF] Ponte Castañeda | The effective mechanical properties of nonlinear isotropic solids[END_REF]. The method [START_REF] Ponte Castañeda | The effective mechanical properties of nonlinear isotropic solids[END_REF] makes use of a linear comparison composite with the same microstructure as the given nonlinear composite. Most crucially, any bound on the effective energy of the linear comparison composite can be used to generate a corresponding bound on the energy of the nonlinear composite. In particular, nonlinear bounds that take 3-point statistics into account can be obtained [START_REF] Ponte Castañeda | Bounds and estimates for the properties of nonlinear heterogeneous systems[END_REF]. When used in conjunction with the linear Hashin-Shtrikman bounds, the approach of [START_REF] Ponte Castañeda | The effective mechanical properties of nonlinear isotropic solids[END_REF][START_REF] Ponte Castañeda | Bounds and estimates for the properties of nonlinear heterogeneous systems[END_REF]) recovers1 the bounds of [START_REF] Talbot | Variational principles for inhomogeneous non-linear media[END_REF]. Those bounds will be referred to as the linear comparison bounds in the following. In contrast with the linear case, the linear comparison bounds are generally not optimal and can a priori be improved. For viscoplastic composites, Ponte Castañeda (2012) has shown that improved bounds could be obtained via incremental homogenization. Strictly speaking, the bounds of Ponte Castañeda (2012) hold only for a special class of isotropic microstructures that are compatible with the incremental homogenization procedure (see [START_REF] Ponte Castañeda | Bounds for nonlinear composites via iterated homogenization[END_REF] for details). Aside from rigorous bounds, the issue of finding estimates of the effective behavior has also received a lot of attention [START_REF] Maxwell-Garnett | Colours in metal glasses and in metal films[END_REF]Bruggeman, 1935;[START_REF] Norris | A generalized differential effective medium theory[END_REF][START_REF] Ponte Castañeda | Nonlinear electrical conductivity in heterogeneous media[END_REF][START_REF] Ponte Castañeda | Second-order homogenization estimates for nonlinear composites incorporating field fluctuations: I-theory[END_REF][START_REF] Idiart | Estimates for two-phase nonlinear conductors via iterated homogenization[END_REF]. As pointed out notably by [START_REF] Idiart | Estimates for two-phase nonlinear conductors via iterated homogenization[END_REF], bounds and estimates are complementary. In particular, bounds allow for the validity of estimates to be gauged.

In this paper, we use the 'translation' method (Milton, 1990) to derive rigorous bounds that hold for the whole class of isotropic two-dimensional composites and improve on the linear comparison bounds. Originally introduced by [START_REF] Lurie | Exact estimates of conductivity of composites formed by two isotropically conducting media taken in prescribed proportion[END_REF] and independently by [START_REF] Tartar | Estimations fines des coefficients homogénéisés[END_REF]; [START_REF] Murat | Calcul des variations et homogénéisation[END_REF], the translation method has proved to be a powerful tool in various nonlinear homogenization problems, such as bounding the stress-strain curve of composites [START_REF] Milton | Bounding the current in nonlinear conducting compos[END_REF][START_REF] Talbot | Bounds for the effective constitutive relation of a composite[END_REF][START_REF] Peigney | A pattern-based method for bounding the effective response of a nonlinear composite[END_REF] or bounding the recoverable strains in shape memory alloys [START_REF] Bhattacharya | Energy minimization and the recoverable strains in polycrystalline shape memory alloys[END_REF][START_REF] Peigney | Recoverable strains in composite shape-memory alloys[END_REF][START_REF] Peigney | Improved bounds on the energy-minimizing strains in martensitic poly-Continuum Mech[END_REF]. For linear composites, the Hashin-Shtrikman bounds can be elegantly recovered from the translation method by embedding the original problem in a problem of higher dimension (namely, d 2 where d is the dimension of the physical space). That embedding is achieved by considering simultaneously d loading orientations. For nonlinear composites, combining the translation method with the idea of embedding has proved to be a successful approach notably for bounding the breakdown limit of dielectric composites or the yield limit of ideally-plastic polycrystals [START_REF] Kohn | Some problems of polycrystal plasticity with deficient basic crystals[END_REF][START_REF] Nesi | Translation and related bounds for the response of a nonlinear composite conductor[END_REF][START_REF] Nesi | Improved bounds for the yield stress of a polycrystalline material[END_REF]Garroni and Kohn, 2003). [START_REF] Nesi | Translation and related bounds for the response of a nonlinear composite conductor[END_REF] proved that the linear comparison bound cannot be better than the bounds obtained from the translation method. Those authors compared the two bounds on some examples of twodimensional composites and observed that they actually coincide. In those comparisons, [START_REF] Nesi | Translation and related bounds for the response of a nonlinear composite conductor[END_REF] used the translation bound obtained by embedding the original problem in R 4 , as considered originally by [START_REF] Lurie | Exact estimates of conductivity of composites formed by two isotropically conducting media taken in prescribed proportion[END_REF] and [START_REF] Tartar | Estimations fines des coefficients homogénéisés[END_REF]; [START_REF] Murat | Calcul des variations et homogénéisation[END_REF] in the linear setting. In this paper, we consider an embedding of the original two-dimensional problem in R 6 , instead of R 4 . This is achieved by considering 3 loading orientations, instead of 2. The corresponding bounds are expected to be stronger than the earlier translation bounds (using 2 loading orientations) because they actually apply to a smaller class of composites (namely the class of composites having the same response along 3 loading orientations, rather than 2). Although the calculations become more difficult, the bounds obtained are shown to often improve on the linear comparison bounds.

The general methodology is introduced in Sect. 2, with special emphasis on the problem of deriving lower bounds on the effective energy. The relevance of the proposed approach is demonstrated in Sect. 3 by considering some examples for which explicit results can be obtained. Embedding in R 6 has the critical downside of dramatically increasing the complexity of the calculations. In Sect. 4 we use invariance properties to show that major simplifications occur when the individual phases are isotropic. The developments of Sect. 4 are essential to make the approach tractable and to be able to evaluate the bounds in practice. The proposed approach is applied in Sect. 5 to power law composites. In Sect. 6, the methodology is adapted to bound the effective energy from above. where w is the convex energy-density function, depending both on e and on the location x. Denoting by ē (resp. ) the spatial average of e (resp. j), it can be shown [START_REF] Hill | Elastic properties of reinforced solids: some theoretical principles[END_REF][START_REF] Willis | Variational estimates for the overall response of an inhomogeneous nonlinear dielectric[END_REF] that ē and  are related by the effective constitutive law

A general lower bound on the effective energy

 = dw eff dē (ē) (2.2)
where w eff is the effective energy function of the composite material, defined by

w eff (ē) = inf e∈K(ē) Ω w(e, x)dω (2.3) with K(ē) = {e : Ω → R 2 |e = ∇V for some V : Ω → R verifying V (x) = ē.x on ∂Ω}.
Determining w eff (ē) formally amounts to solve the set of local equations div j = 0 in Ω , j = ∂w ∂e (e, x), e ∈ K(ē).

(2.4)

The problem (2.1)-(2.4) arises in many physical contexts, such as thermal conductivity, magnetic permeability or diffusion. Alternatively, the local constitutive law (2.1) can be rewritten as e = ∂w * ∂j (j, x)

where the complementary energy-density function w * is the Legendre dual of w. Similarly, the effective constitutive law (2.2) can be expressed as

ē = dw * eff d ()
where the effective complementary energy function w * eff is the Legendre dual of w eff and satisfies

w * eff () = inf j∈S() Ω w * (j, x)dω (2.5) with S() = {j : Ω → R 2 |div j = 0 in Ω; j(x) • n =  • n on ∂Ω}.
This paper is concerned with deriving bounds on w eff and w * eff , with special emphasis on the case where the overall behavior is isotropic. Since w eff and w * eff are Legendre duals of each other, a lower (resp. upper) bound on w eff translates as an upper (resp. lower) bound on w * eff , and vice versa. In order to simplify the presentation, we mainly focus on lower bounds on w eff . The issue of deriving upper bounds on w eff is addressed in Sect. 6.

Extended effective energy

Consider n distinct values ē1 , • • • , ēn and let e i (x) be the solution of (2.4) for ē = ēi . It can be seen from ( 2.3) that the extended field

E = (e 1 , • • • , e n ) solves the minimization problem inf E∈K( Ē) Ω n i=1 w(e i , x) dω (2.6) where Ē = (ē 1 , • • • , ēn ) and K( Ē) = {(e 1 , • • • , e n ) : e i = ∇V i for some V i : Ω → R verifying V i (x) = ēi .x on ∂Ω}.
In view of (2.6), it is convenient to introduce the extended energy W defined as

W (E, x) = n i=1
w(e i , x).

Similarly, we define the extended effective energy W eff by

W eff ( Ē) = n i=1
w eff (ē i ).

(2.7)

Combining (2.3) and (2.6) leads to

W eff ( Ē) = inf E∈K( Ē) Ω W (E, x) dω.
(2.8)

The variational structure (2.8) is similar to that of w eff in (2.3), except it is defined in terms of higher dimensional fields.

Translated bounds on the extended energy

We now use the translation method to generate bounds on W eff . Let U (E) be a given function (referred to in the following as a comparison potential) and consider the Legendre transform

(W -U ) * (T , x) = sup E E • T -W (E, x) + U (E) (2.9)
for any T in R 2n . The dot • in (2.9) denotes the Euclidean scalar product in R 2n . It follows from (2.9) that

W (E, x) ≥ E • T + U (E) -(W -U ) * (T , x) (2.10)
for any E and T in R 2n . Taking T independent on x and integrating (2.10) over the domain Ω, we find

Ω W (E, x)dω ≥ Ē • T + Ω U (E)dω - Ω (W -U ) * (T , x)dω (2.11)
for any E ∈ K( Ē). Assume now that the comparison potential U is quasiconvex, i.e. satisfies

Ω U (E)dω ≥ U ( Ē) for any E ∈ K( Ē).
(2.12)

Replacing in (2.11) yields

Ω W (E, x)dω ≥ Ē • T + U ( Ē) - Ω (W -U ) * (T , x)dω (2.13)
for any E ∈ K( Ē). Since the right-hand side of (2.13) is independent on E, taking the infimum over E ∈ K( Ē) gives

W eff ( Ē) ≥ Ē • T + U ( Ē) - Ω (W -U ) * (T , x)dω.
(2.14)

The right-hand side of (2.14) is thus a lower bound on W eff . Each choice of comparison potential U generates a corresponding bound. However, U should be carefully chosen for the bound (2.14) to be meaningful. In that regard, a first requirement is that (W -U ) * (T , x) < +∞, which sets restrictions on the behavior of U at infinity.

The isotropic case

The bound (2.14) applies to any microstructure. We assume from now on that the given composite is isotropic, in which case w eff only depends on the electric field ē through its norm ē. In other words, there exists a function φ eff such that w eff (ē) = φ eff (ē).

(2.15)

Take Ē = (ē 1 , • • • , ēn ) such that ē1 ,• • • ,ē n have the same norm (denoted by ē), i.e. ē1 = • • • = ēn = ē.
(2.16)

We have

w eff (ē 1 ) = • • • = w eff (ē n ) = φ eff (ē),
(2.17) hence Eq. (2.14) implies the following lower bound on φ eff (ē):

φ eff (ē) ≥ 1 n Ē • T + U ( Ē) - Ω (W -U ) * (T , x)dω .
(2.18)

In the derivation exposed, the isotropy is captured only by the relation (2.17), i.e. by stating that the effective energy w eff takes the same value along n prescribed directions having the same norm. For linear materials, the condition (2.17) fully characterizes isotropy as soon as n ≥ 2 (at least for well-chosen directions e i ). For nonlinear behaviour, however, the relation (2.17) is a necessary but not sufficient condition for isotropy to occur. In such case, the bound delivered by (2.18) is expected to become tighter as n increases. [START_REF] Nesi | Translation and related bounds for the response of a nonlinear composite conductor[END_REF] studied bounds of the form (2.18) with n = 2. The comparison potential U was chosen as the determinant in R 2 . The resulting translated bound was shown by [START_REF] Nesi | Translation and related bounds for the response of a nonlinear composite conductor[END_REF] to often coincide with the linear comparison bound. In this paper, we take n = 3 i.e. we consider 3 loading orientations instead of 2. The resulting bounds are expected to be sharper than the translation bounds of [START_REF] Nesi | Translation and related bounds for the response of a nonlinear composite conductor[END_REF] because they apply to a smaller class of composites, namely the class of composites having the same response along 3 different orientations rather than 2. In more detail, we consider comparison potentials of the form

U (E) = αV (E) (2.19)
where α is a parameter and

V (E) = 3 i, j = 1 i < j (u i v j -u j v i ).
(2.20)

In (2.20), (u i , v i ) denotes the components of the vector e i in R 2 , so that

E = (u 1 , v 1 , u 2 , v 2 , u 3 , v 3 ).
For any fixed i and j, the determinant u i v j -u j v i is known to be a null Lagrangian, i.e. it satisfies (2.12) as an equality. Consequently, the function (2.19) is also a null Lagrangian for any value of α. The potential U in (2.19) is quadratic in E and can be written in the form

U (E) = 1 2 αE • M • E
where M is the symmetric matrix given by

M =         0 0 0 1 0 1 0 0 -1 0 -1 0 0 -1 0 0 0 1 1 0 0 0 -1 0 0 -1 0 -1 0 0 1 0 1 0 0 0         . (2.21)
For later reference, we note that the eigenvalues of M are ± √ 3 and 0.

We will use (2.18) with

Ē = ēN (2.22) where N = 1 2 (1, - √ 3, -1, - √ 3, -2, 0).
(2.23)

The vector Ē = ēN satisfies (2.16) and the three corresponding directions ēi make a π/3 angle between each other (Fig. 1). The main motivation for the choice (2.22) is that it allows the linear Hashin-Shtrikman bound to be recovered, as will be seen in Sect. 3.1. It can easily be verified that N is an eigenvector of M for the eigenvalue -√ 3. Therefore, using (2.18) with Ē of the form (2.22), we obtain

φ eff (ē) ≥ ē 3 N • T - √ 3 2 αē 2 - 1 3 Ω (W -αV ) * (T , x) dω.
(2.24)

The bound (2.24) holds for any α ∈ R and T ∈ R 6 , hence

φ eff (ē) ≥ sup α∈R,T ∈R 6 ē 3 N • T - √ 3 2 αē 2 - 1 3 Ω (W -αV ) * (T , x) dω .
(2.25)

Calculations can be made simpler if the vector T in (2.25) is restricted to be parallel to N , i.e. of the form T = τ N . This generates the following lower bound: For an isotropic composite, the effective energies w eff (ē1), w eff (ē2), w eff (ē3) take the same value.

φ eff (ē) ≥ sup α∈R,τ ∈R ēτ - √ 3 2 αē 2 - 1 3 Ω (W -αV ) * (τ N , x) dω . (2.26) -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 π/3 π/3 ē1 ē2 ē3
In general, the bound (2.26) is weaker than the bound (2.25). We will show in Sect.4 that the two bounds actually coincide if all the constitutive phases have an isotropic behavior. It should also be noted that the bounds (2.25) and (2.26) are meaningful only if w(e, x) has faster than quadratic growth in e (so as to ensure that (W -αV ) * (τ N , x) remains finite).

In the next section, the relevance of the bounds (2.25) and (2.26) is demonstrated on some examples for which exact expressions can be obtained.

Exact results

The linear case

It is insightful to first consider the case where w(e, x) = 1 2 e • σ(x) • e and σ(x) is symmetric positive definite. The corresponding extended energy

W (E, x) is equal to 1 2 E • Σ(x) • E where Σ(x) has the block structure Σ(x) =   σ(x) 0 0 0 σ(x) 0 0 0 σ(x)   .
For an isotropic composite, the effective energy w eff (ē) can be written as 1 2 σ eff ē2 where σ eff is the effective conductivity. A lower bound on σ eff can be obtained by using the prescription (2.25) with ē = 1. Provided that Σ(x) -αM is positive definite, (W -αV ) * (T , x) is finite and given by

(W -αV ) * (T , x) = 1 2 T • (Σ(x) -αM) -1 • T .
It follows from (2.24) that

σ eff ≥ 2 3 N • T - √ 3α - 1 3 T • Ω (Σ(x) -αM) -1 dω • T . (3.1)
Maximizing the right hand side of (3.1) with respect to T , we obtain

σ eff ≥ 1 3 N • Ω (Σ(x) -αM) -1 dω -1 • N - √ 3α. (3.2)
We note that the value of T realizing the maximum in (3.1) is

T = Ω (Σ(x) -αM) -1 dω -1 • N . (3.3)
Let us specialize the bound (3.2) in the case where σ(x) is isotropic, i.e. of the form σ(x)1 where σ(x) > 0 and

1 is the identity matrix in R 2 . Since (Σ(x) -αM) • N = (σ(x) + √ 3α)N , a direct calculation from (3.2) leads to σ eff ≥ Ω 1 σ(x) + √ 3α dω -1 - √ 3α. (3.4)
The relation (3.4) holds provided that Σ(x) -αM is positive definite. As noted in Sect. 2.3, the eigenvalues of M are 0 and ± √ 3. Hence Σ(x) -αM is positive definite if |α| < σ(x)/ √ 3. Letting α tend towards inf x σ(x)/ √ 3 from below in (3.4), we obtain

σ eff ≥ Ω 1 σ(x) + inf x σ(x) dω -1 -inf x σ(x).
(3.5)

The formula (3.5) corresponds to the celebrated bound of [START_REF] Hashin | A variational approach to the theory of the effective magnetic permeability of multiphase materials[END_REF]. For instance, in the case of a two-phase composite, (3.5) specializes as

σ eff ≥ c 1 σ 1 + c 2 σ 2 - c 1 c 2 (σ 1 -σ 2 ) 2 σ 1 + c 2 σ 1 + c 1 σ 2 (3.6)
where σ i and c i are respectively the conductivity and the volume fraction of phase i. In (3.6), the phases are labelled in such fashion that σ 1 < σ 2 . The fact that the linear Hashin-Shtrikman bound is recovered allows us to check that the choice (2.22) of Ē is consistent with the potentials U considered. Others choice would do, but Ē cannot be taken in a completely arbitrary fashion (the important point is that

M • Ē = - √ 3 Ē).
In the case of linear isotropic constitutive materials, observe from (3.3) that the optimal value of T is parallel to N , hence nothing is lost by considering (2.26) instead of (2.25). As will be shown in Sect. 4, that property remains true in the nonlinear case.

Dielectric breakdown in a two-phase composite

Consider a two-phase nonlinear composite in which the energy w i of phase i (i = 1, 2) is given by

w 1 (e) = 0 if e ≤ 1, +∞ if e > 1.
(3.7)

w 2 (e) = 0 if e = 0, +∞ otherwise (3.8)
Phase 1 behaves as a perfect dielectric so long as e ≤ 1, and breaks down when e reaches 1. Phase 2 is a perfect conductor. Assuming overall isotropy, the effective energy w eff (ē) takes the form

w eff (ē) = 0 if ē ≤ k eff , +∞ if ē > k eff ;
(3.9)

i.e. the composite material behaves as a perfect dielectric as long as ē < k eff and breaks down when the magnitude of the electric field ē reaches the critical value k eff . For later reference, we note that (3.9) corresponds to an effective complementary energy w * eff of the form

w * eff () = k eff . (3.10)
where  is the magnitude of . The parameter k eff in (3.9) and (3.10) is referred to as the effective strength of the composite dielectric. The formula (2.26) is now use to bound k eff .

Denoting the volume fraction of phase i by c i , Eqs. (2.26) and (3.9) imply that any k in [0,

k eff ] satisfies 0 ≥ kτ - √ 3 2 αk 2 - 2 i=1 c i 3 (W i -αV ) * (τ N ) ∀α > 0, τ ≥ 0; (3.11) with (W i -αV ) * (τ N ) = sup E=(e 1 ,e 2 ,e 3 )    τ N • E + 1 2 αE • M • E - 3 j=1 w i (e j )    .
(3.12)

From (3.7) and (3.8) we have

(W 1 -αV ) * (τ N ) = sup E = (e 1 , e 2 , e 3 ) e i ≤ 1 τ N • E + 1 2 αE • M • E (3.13)
and (W 2 -αV ) * (τ N ) = 0. We can easily get rid of the parameter α in (3.11): For any non-negative α, observe indeed from (3.13) that (W 1 -αV ) * (τ N ) = α g * (τ /α) where g * is defined for τ > 0 by

g * (τ ) = sup E = (e 1 , e 2 , e 3 ) e i ≤ 1 τ N • E + 1 2 E • M • E . (3.14)
Further setting

P (k, τ ) = kτ - √ 3 2 k 2 - c 1 3 g * (τ ), we have sup α>0,τ ≥0 kτ - √ 3 2 αk 2 - 2 i=1 c i 3 (W i -αV ) * (τ N ) = sup α>0,τ ≥0 αP (k, τ α ) = sup α>0,τ ≥0 αP (k, τ ) =    +∞ if 0 < sup τ ≥0 P (k, τ ) 0 if 0 ≥ sup τ ≥0 P (k, τ ).
Hence the condition (3.11) is equivalent to

0 ≥ sup τ ≥0 P (k, τ ).
Any k verifying 0 < sup τ ≥0 P (k, τ ) is necessary above k eff . Therefore, an upper bound k on k eff is given by

k = inf k : 0 < sup τ ≥0 P (k, τ ) . (3.15)
Evaluating k requires some details on the function g * introduced in (3.14). Note that g * is convex because it is defined as the pointwise supremum of a family of linear functions. The constrained quadratic optimization problem (3.14) can be solved in closed form as detailed in Appendix A. Setting

τ 0 = 1 + √ 3, it is shown in Appendix A that: g * (τ ) = τ (cos θ - √ 3 sin θ + 1) + 2(cos θ + 1) sin θ if 0 ≤ τ ≤ τ 0 , 3(τ - √ 3 2 ) if τ ≥ τ 0 , (3.16)
where θ is characterized by the conditions

sin θ =      √ 3 2 if τ = 0, 2 τ (2 cos 2 θ -1 + cos θ) - √ 3 cos θ if τ > 0, cos θ = z - 1 2 + τ 2 √ 3 , (3.17)
and z is the largest root of the cubic equation At τ = τ 0 , g * is subdifferentiable [START_REF] Rockafellar | Convex Analysis[END_REF] and we have

8z 3 -6z + √ 3τ 3 9 -2 = 0. (3.18) 0 1 2 3 4 2 4 6 8 10 τ g * 0 1 2 3 4 0 1 2 3 4 τ 0 τ 0 τ ∂g *
∂g * (τ 0 ) = lim τ -→τ - 0 dg * dt (τ ), lim τ -→τ + 0 dg * dt (τ ) = [τ 0 , 3]. (3.20)
The subdifferential ∂g * is represented in Fig. 2(right). The nondifferentiability of g * at τ = τ 0 corresponds to a (hardly noticeable) discontinuity in the slope of the curve τ → g * (τ ) shown in Fig. 2(left).

We are now in a position to evaluate the bound k in (3.15). By (3.16) we have sup

τ ≥τ 0 P (k, τ ) = +∞ if k > c 1 , P (k, τ 0 ) if k ≤ c 1 .
Hence (3.15) can be rewritten as

k = inf k ≤ c 1 : 0 < sup 0≤τ ≤τ 0 P (k, τ ) . (3.
     τ * = τ 0 if k ≥ c 1 3 τ 0 , τ * < τ 0 , 3 k c 1 = dg * dτ (τ * ) if k < c 1 3 τ 0 .
(3.24)

We first consider the case k

≥ c 1 3 τ 0 . In such case, the requirement (3.22) is equivalent to P ( k, τ 0 ) = 0. Since g * (τ 0 ) = 3(1 + √ 3/2), we obtain 0 = (1 + √ 3) k - √ 3 2 k2 -c 1 (1 + √ 3 2 ) which yields k = 1 √ 3 1 + √ 3 -1 + (3 + 2 √ 3)(1 -c 1 ) . (3.25)
That expression holds provided it is compatible with the restriction k ≥ c 1 3 τ 0 , i.e. for c 1 ≥ √ 3/2. Let us now consider the case k < c 1 3 τ 0 . We have P ( k, τ * ) = 0, i.e.

0 = kτ * - √ 3 2 k2 - c 1 3 g * (τ * ) (3.26)
where τ * < τ 0 satisfies (3.24), i.e.

a + b = 3 c 1 k , ab = 9 4c 1 k2 .
It follows that

a = 3 k 2c 1 (1 + √ 1 -c 1 ), b = 3 k 2c 1 (1 - √ 1 -c 1 ) (3.29) for ∈ {-1, 1}. Rewriting the equality cos 2 θ * + sin 2 θ * = 1 in terms of (a, b) leads to k = c 1 (1 + √ 1 -c 1 ) 2 -c 1 + √ 1 -c 1 . (3.30)
The right-hand side of (3.30) is a monotonically increasing function of c 1 and is equal to

c 1 τ 0 /3 for c 1 = √ 3/2.
In summary, the final expression of the bound k on the effective strength of the dielectric composite is thus

k =      c 1 (1 + √ 1 -c 1 ) 2 -c 1 + √ 1 -c 1 for c 1 ≤ √ 3/2, 1 √ 3 1 + √ 3 -1 + (3 + 2 √ 3)(1 -c 1 ) for c 1 ≥ √ 3/2. (3.31)
The bound k given by (3.31) is represented in Fig. 3, along with the Wiener bound k W and the linear comparison bound k LC , respectively given by [START_REF] Ponte Castañeda | Bounds and estimates for the properties of nonlinear heterogeneous systems[END_REF])

k W = c 1 , k LC = c 1 √ 2 -c 1 .
We can observe that k strictly improves on the linear comparison bound. The maximum improvement is about 6%.

Case of isotropic constituents

In general there are substantial difficulties in evaluating the bound in (2.25). For any given T and α, the formula (2.25) indeed requires to evaluate the quantity (W -αV ) * (T , x) defined in (2.9), i.e. to solve a nonlinear optimization problem in R 6 . Since V is nonconvex, local optima do not necessarily coincide with global optima and multiple solutions can exist. In addition, as illustrated in Sect. 3, the 6-dimensional vector T in (2.25) as well as the scalar parameter α should be optimized upon in order to get relevant bounds.

We show that major simplifications occur when the constituent materials of the composite are isotropic, i.e. when the local energy function w(e, x) can be written as where e denotes the norm of the vector e. Under such assumption, we show in Sect. 4.2 that the values of T reaching the supremum in (2.25) are necessarily parallel to the vector N introduced in (2.23). A first consequence is that the optimization with respect to T in (2.25) is reduced from a 6-dimensional to a 1-dimensional problem. As a second consequence, only the values taken by (W -αV ) * (T , x) for T parallel to N are needed for evaluating the optimized bound. For such values of T , the nonconvex optimization problem that defines (W -αV ) * (T , x) can be considerably simplified, as explained in Sect. 4.3. In all the developments presented next, a central role is played by a group of linear transformations that leave V , W and N invariant. Those transformations are first introduced in Sect. 4.1.

w(e, x) = φ(e, x) (4.1) 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 c 1 k eff k W k LC k

Invariance properties

Consider the linear transformations defined on R 6 by

R 1 • (e 1 , e 2 , e 3 ) = ( r( π 3 ) • e 2 , r( π 3 ) • e 3 , r( -2π 3 ) • e 1 ) R 2 • (e 1 , e 2 , e 3 ) = ( r( 2π 3 ) • e 3 , r( -π 3 ) • e 1 , r( -π 3 ) • e 2 ) S 1 • (e 1 , e 2 , e 3 ) = ( s( -2π 3 ) • e 1 , s( π 3 ) • e 3 , s( π 3 ) • e 2 ) S 2 • (e 1 , e 2 , e 3 ) = ( s( 2π 3 ) • e 3 , s( 2π 3 ) • e 2 , s( 2π 3 ) • e 1 ) S 3 • (e 1 , e 2 , e 3 ) = ( s(π) • e 2 , s(π) • e 1 , s(0) • e 3 ) (4.2) • I R 1 R 2 S 1 S 2 S 3 I I R 1 R 2 S 1 S 2 S 3 R 1 R 1 R 2 I S 2 S 3 S 1 R 2 R 2 I R 1 S 3 S 1 S 2 S 1 S 1 S 3 S 2 I R 2 R 1 S 2 S 2 S 1 S 3 R 1 I R 2 S 3 S 3 S 2 S 1 R 2 R 1 I Table 1: Group table of G.
where

r(α) = cos α -sin α sin α cos α
is the two-dimensional rotation with angle α, and

s(α) = cos α sin α sin α -cos α
is the two-dimensional reflection through the axis (cos(α/2), sin(α/2)). The transformations in (4.2) form a group G whose table is shown in Table 1. The action of G on a vector (e 1 , e 2 , e 3 ) is represented in Figure 4.

It can be verified that the function U in (2.19) is invariant by G, i.e. that

G t • M • G = M for all G ∈ G, (4.3)
where the superscript t denotes the transpose operator. There are other linear transformations satisfying (4.3) but those in G have the additional property of leaving N invariant, i.e. we have

G • N = N for all G ∈ G. (4.4)
Also observe that W is invariant under G as the transformations in G have the effect of permuting the norms of the vectors e 1 , e 2 , e 3 .

We further note that all the transformations in G are isometries of R 6 , i.e. they satisfy the relation

G t • G = I. More precisely: • R 1 and R 2 are 'rotations', i.e. det R i = 1. Their invariant space is vect(N , N ) with N = 1 2 (- √ 3, -1, - √ 3, 1, 0, 2). (4.5)
• S 1 , S 2 and S 3 are 'reflections', i.e. S i • S i = I. We have already noted that S i • N = N , and it can be checked that The vector (e 1 , e 2 , e 3 ) shows the transformed of (e1, e2, e3) for each transformation in G. The transformations Ri involve some rotations (with angle in {±π/3, ±2π/3}) and some permutations of the vectors (e1, e2, e3). The transformations Si involve some reflections (through axis shown as dashed lines) and some permutations of the vectors (e1, e2, e3).

S i • N = -N for i = 1, 2, 3. ( 4
For our purpose, an important observation is that

I + R 1 + R 2 + S 1 + S 2 + S 3 = 2N ⊗ N .
(4.7)

We close this section with a justification of the identity (4.7). Take any arbitrary E in R 6 and set

E = (I + R 1 + R 2 ) • E. It follows from Table 1 that R 1 • E = (R 1 + R 2 1 + R 1 • R 2 ) • E = (R 1 + R 2 + I) • E = E , i.e. that E is invariant by R 1 . Hence E ∈ vect(N , N ). This shows that I + R 1 + R 2
takes values in vect(N , N ) and therefore can be written as

I + R 1 + R 2 = N ⊗ a + N ⊗ a (4.8)
for some a and a in R 6 . Taking the transpose of (4.8) and noting from Table 1

that R t 1 = R -1 1 = R 2 , we obtain I + R 1 + R 2 = a ⊗ N + a ⊗ N .
(4.9)

Applying I + R 1 + R 2 to N and observing that N • N = 0 shows that a = N . Similarly, applying I + R 1 + R 2 to N gives a = N . We thus arrive at

I + R 1 + R 2 = N ⊗ N + N ⊗ N .
Now right-multiplying that last equation by I + S 3 gives

(I + R 1 + R 2 ).(I + S 3 ) = N ⊗ (N + S 3 • N ) + N ⊗ (N + S 3 • N ) = 2N ⊗ N
where the relation (4.6) has been used. Expanding (I + R 1 + R 2 ).(I + S 3 ) and using Table 1 shows that (

I + R 1 + R 2 ).(I + S 3 ) = I + R 1 + R 2 + S 1 + S 2 + S 3
, which gives the identify (4.7).

Optimization with respect to T

For a given α, the best bound available from (2.25) is sup Similarly, for a given α, the best bound available from (2.26) is sup

T ∈R 6 ē 3 N • T - √ 3 2 αē 2 -F (T ) (4.10) with F (T ) = 1 3 Ω (W -αV ) * (T , x) dω. ( 4
τ ∈R ēτ - √ 3 2 αē 2 -f (τ ) (4.13) where f (τ ) = F (τ N ). (4.14)
The optimal values of τ in (4.13) are characterized by ē ∈ ∂f (τ ). (4.15)

Observe from the definition of f that

∂f (τ ) = N • ∂F (τ N ). (4.16)
Consequently, if T satisfies the optimality condition (4.12) then the scalar τ = T • N satisfies (4.15). We show in the following that the converse is true: If τ satisfies (4.15) then the vector T = τ N satisfies (4.12). This is not obvious at first sight because (4.15) is only a trace of (4.12) on a particular direction (namely, N ). The key arguments leading to the claimed result lie in the invariance properties studied in Sect. 4.1. A first observation is that the function F in (4.11) is invariant by

G. Since any G ∈ G verifies G t • G = I, we have indeed (W -αV ) * (T , x) = sup E {GE • GT -W (E, x) + αV (E)} = sup Ẽ Ẽ • GT -W (G -1 Ẽ, x) + αV (G -1 Ẽ) (4.17)
where Ẽ = G • E. As W and V are invariant by G, (4.17) gives (W -αV

) * (GT , x) = (W -αV ) * (T , x). It follows that F is invariant by G. A notable consequence is that G • ∂F (E) ⊂ ∂F (G • E) (4.18)
for any E and G ∈ G. From the definition of the subdifferential, any A ∈ ∂F (E) indeed satisfies the distinctive property

F (E ) -F (E) ≥ A • (E -E) (4.19)
where E is arbitrary. Since G t • G = I and F is invariant by G, (4.19) can be rewritten as

F (G • E ) -F (G • E) ≥ G t • G • A • (E -E) = (G • A) • (G • (E -E))
which shows that G • A ∈ ∂F (G • E), thus proving (4.18). Now let τ be a scalar verifying the optimality condition (4.15). From (4.16) we have

ē ∈ N • ∂F (τ N ). (4.20)
Since N is invariant by G, using (4.18) with E = τ N gives

G • ∂F (τ N ) ⊂ ∂F (τ N ) for all G ∈ G. (4.21)
The function F being convex, ∂F (E) is a convex set for all E [START_REF] Rockafellar | Convex Analysis[END_REF]. Hence (4.21) implies that

1 6 (I + R 1 + R 2 + S 1 + S 2 + S 3 ) • ∂F (τ N ) ⊂ ∂F (τ N )
which using the identity (4.7) can be rewritten as

1 3 (N • ∂F (τ N ))N ⊂ ∂F (τ N ). (4.22)
In view of (4.20), we thus obtain that

1 3 ēN ∈ ∂F (τ N ).
This proves that the vector T = τ N satisfies (4.12) and therefore reaches the supremum in (4.10). The conclusion is that the bounds provided by (4.10) and (4.13) coincide. Consequently, the bounds (2.25) and (2.26) are identical. In practice, the prescription (2.26) is simpler to use as it only involves 2 scalar parameters (τ, α) -instead of 7 in (2.25) .

Legendre transform for T parallel to N

In this section we address the evaluation of the quantity (W -αV ) * (τ N , x) in (2.9), which is central for the evaluation of the bound (2.26). With proper rescaling, we can assume that α = 1. For a given x, (W -V ) * (τ N , x) is obtained by the solving the maximization problem sup E∈R 6

F(E) (4.23)
where

F(E) = τ E • N -W (E, x) + 1 2 E • M • E. (4.24)
Eq. (4.23) is a nonlinear optimization problem in R 6 . Except in few special cases, there is no hope to solve the problem (4.23) in closed-form and one needs to resort to numerical techniques. Standard gradient algorithms are able to provide a local optimum but not necessarily a global optimum as needed for our purpose. This issue is especially sensitive as the function F in (4.24) is non convex and expected to have multiple stationary points. Global maximization techniques are likely to overcome such difficulties, but those algorithms are much more costly than gradient algorithms. This is especially critical as the (numerical) optimization with respect to (α, T ) in (2.26) requires numerous evaluations of (W -αV ) * (τ N , x). As a result, applying a global algorithm for solving (4.23) entails prohibitive computational costs when evaluating the bound in (2.26). We show in the following that the problem can actually be reduced to a nonlinear optimization problem in R 2 , which makes the approach tractable.

Under the assumption (4.1) that w(., x) is isotropic, the stationarity condition ∂F/∂E = 0 can be written as

(β 1 e 1 , β 2 e 2 , β 3 e 3 ) -M • E = τ N (4.25)
where

β i = 1 e i ∂φ ∂e (e i , x). (4.26)
Projecting (4.25) on the canonical basis of R 6 gives the 6 equations

β 1 u 1 -v 2 -v 3 = τ 2 , (4.27a) β 1 v 1 + u 2 + u 3 = - √ 3 2 τ, (4.27b) β 2 u 2 + v 1 -v 3 = - τ 2 , (4.27c) β 2 v 2 -u 1 + u 3 = - √ 3 2 τ, (4.27d
)

β 3 u 3 + v 1 + v 2 = -τ, (4.27e
)

β 3 v 3 -u 1 -u 2 = 0. (4.27f)
Since β i in (4.26) is a function of e i , the system (4.27) is nonlinear in e i and cannot be solved in closed form. However, any solution of (4.27) can be written as an explicit function of the two scalar components u 3 and v 3 , as is now explained. We first rewrite the two equations (4.27b-c) as a linear system in (v 1 , u 2 ), namely

β 1 1 1 β 2 v 1 u 2 = - √ 3 2 τ -u 3 -τ 2 + v 3 . (4.28)
Provided that β 1 β 2 -1 = 0, the linear system (4.28) can be inverted to express v 1 and u 2 as functions of (τ, β 1 , β 2 , u 3 , v 3 ). The corresponding expressions are

v 1 = 1 Y (-β 2 p + τ 2 -v 3 ), u 2 = 1 Y (p -β 1 ( τ 2 -v 3 )); (4.29)
where

Y = β 1 β 2 -1, p = u 3 + √ 3 2 τ. (4.30)
Similarly, we obtain from (4.27a) and (4.27d) that

u 1 = 1 Y (β 2 ( τ 2 + v 3 ) -p), v 2 = 1 Y ( τ 2 + v 3 -β 1 p). ( 4
p p -τ 2 + v 3 τ 2 + v 3 β 1 β 2 = Y (β 3 u 3 + τ ) + τ β 3 v 3 Y . (4.32)
Inverting (4.32) with respect to (β 1 , β 2 ) gives, provided that pτ = 0,

β 1 = ( τ 2 + v 3 )(Y (β 3 u 3 + τ ) + τ ) -β 3 v 3 pY τ p , β 2 = ( τ 2 -v 3 )(Y (β 3 u 3 + τ ) + τ ) + β 3 v 3 pY τ p .
(4.33)

The point of those algebraic manipulations is to express β 1 and β 2 as explicit functions of (u 3 , v 3 ). Multiplying the two equations in (4.33) indeed gives a relation in which β 1 and β 2 only appear through the product β 1 β 2 . That equation can be written in terms of

Y = β 1 β 2 -1 as a(e 3 )Y 2 -4b(e 3 )Y + 4c(e 3 ) = 0 (4.34)
where a(e 3 ) = (

β 3 u 3 + τ ) 2 -v 2 3 (2 - √ 3β 3 ) 2 , b(e 3 ) = p 2 -β 3 p τ 2 -(2 - √ 3β 3 )( 1 4 τ 2 -v 2 3 ), c(e 3 ) = 1 4 τ 2 -p 2 -v 2 3 . (4.35)
To simplify the presentation, we only consider the most general situation where a(e 3 ) = 0. Some manipulations on (4.35) show that b 2 (e 3 ) -a(e 3 )c(e 3 ) = p 2 q where

q = p 2 + β 3 (τ u 3 + β 3 u 2 3 + 2 √ 3v 2 3 -2β 3 v 2 3 ).
Therefore, if q ≥ 0, the solutions of (4.34) are

Y = 2 b(e 3 ) + |p| √ q
a(e 3 ) where = ±1. Setting Y = Y in (4.33) gives, after some simplification,

β 1 = p + √ q + β 3 v 3 τ -2v 3 + β 3 (u 3 + √ 3v 3 ) , β 2 = p + √ q -β 3 v 3 τ + 2v 3 + β 3 (u 3 - √ 3v 3 ) .
Replacing in (4.29)-(4.31), we finally arrive at

       u 1 = 1 2 (β 3 v 3 + p - √ q), v 1 = 1 2 (-β 3 u 3 -τ + 2v 3 - √ 3β 3 v 3 ), u 2 = 1 2 (β 3 v 3 -p + √ q), v 2 = 1 2 (-β 3 u 3 -τ -2v 3 + √ 3β 3 v 3 ).
(4.36)

It can be verified that the expressions (4.36) remain valid if pτ = 0 or Y = 0. Through (4.36), the components u 1 , v 1 , u 2 , v 2 are all expressed as fully explicit functions of (u 3 , v 3 ) (we recall that β 3 is related to (u 3 , v 3 ) via (4.26)).

For a vector E = (u 1 , v 1 , u 2 , v 2 , u 3 , v 3 ) of the form (4.36), a direct calculation shows that

F(E) = 1 2 √ 3τ 2 + β 3 u 3 (p + √ q) + √ 3β 2 3 v 2 3 - 3 i=1 w( u 2 i + v 2 i , x).
The analysis so far shows that stationary points of F are necessary of the form (4.36) and therefore can be parameterized by (u 3 , v 3 ). In order to evaluate sup R 6 F, we thus only need to look at the values taken by F over E of the form (4.36), i.e. we are left with solving the two optimization problems sup (u 3 ,v 3 ) F + (u 3 , v 3 ) and sup (u 3 ,v 3 ) F -(u 3 , v 3 ) where

F ± (u 3 , v 3 ) = 1 2 ( √ 3τ 2 + β 3 u 3 (p ± √ q) + √ 3β 2 3 v 2 3 ) - 3 i=1 w( u 2 i + v 2 i , x). (4.37)
At this point the optimization problem (4.23) in R 6 has been reduced to two optimization problems in R 2 . Observe from (4.36) that u 2 + u 3 -u 1 + √ 3τ /2 = √ q. Hence maximizing F + amounts to look for solutions to (4.23) that lie in the half-space H + of R 6 defined by the equation u 2 +u 3 -u 1 + √ 3τ /2 ≥ 0. Similarly, maximizing F -amounts to look for solutions in the complementary half-space H -defined by the equation u 2 + u 3 -u 1 + √ 3τ /2 ≤ 0. Using invariance properties, the problem can further be simplified if w is strongly convex, in the sense that the function φ in (4.1) is convex in e 2 . Observe indeed that the function F in (4.24) is invariant by the group of transformations G introduced in Sect. 4.1. Therefore, if E * reaches the supremum in (4.23) then so does R 1 E * and R 2 E * . Assuming that w is strongly convex, it can be proved that E * , R 1 E * , R 2 E * cannot be in the same half-space H ± : if (say)

E * ∈ H + , then either R 1 E * or R 2 E * is in H -. This implies that sup E∈R 6 F(E) = sup (u 3 ,v 3 ) F + (u 3 , v 3 ) = sup (u 3 ,v 3 ) F -(u 3 , v 3 ) (4.38)
Hence sup E∈R 6 F(E) can be obtained by solving indifferently one of the two optimization problems sup (u 3 ,v 3 ) F + or sup (u 3 ,v 3 ) F -. In most cases, numerical computations are needed to carry out that optimization.

Results for power-law composites

In this section are presented some results for two-phase isotropic composites governed by power-law energy functions. The energy function w i in phase i (i = 1, 2) is taken as

w i (e) = σ i n + 1 e n+1 (5.1)
where n > 0 is a nonlinearity index and σ i > 0 is a nonlinear conductivity parameter. The limiting case σ i → ∞ corresponds to a perfect conductor. The case n = 1 corresponds to a linear behavior, whereas the limiting case n → ∞ corresponds to a perfect dielectric as considered in Sect. 3.2. By Legendre duality, the complementary energy functions w * i is given by w

* i (j) = χ i m + 1 j m+1 (5.2)
where m = 1/n and χ i = σ -1/n i is the nonlinear resistivity in phase i. In (5.2), j denotes the norm of the vector j.

The index n in (5.1) is assumed to be the same for both phases. In such case, one can show that the effective energy w eff of the composite is of the form

w eff (ē) = σ eff n + 1 ēn+1
where σ eff is the effective nonlinear conductivity. Correspondingly, the complementary effective energy w * eff is given by

w * eff () = χ eff m + 1 m+1 (5.3)
where the effective resistivity χ eff is related to the effective conductivity σ eff by

χ eff = σ -1/n eff .
Bounding the energy function w eff (or the complementary energy function w * eff ) thus amounts to bounding the nonlinear conductivity σ eff (or the nonlinear resistivity χ eff ). For n ≥ 1, the proposed methodology delivers a lower bound σ on σ eff that is given by

σ n + 1 = sup α∈R,τ ∈R τ - √ 3 2 α - 1 3 2 i=1 c i (W i -αV ) * (τ N ) (5.4)
where (W i -αV ) * is defined as in (3.12). The normalized bound σ/σ 1 is a function of the volume fraction c i , of the ratio σ 2 /σ 1 and of the nonlinearity index n. The (normalized) linear comparison bound σ LC /σ 1 depends on the same parameters and is used as a reference in the following. The linear comparison bound σ LC can be obtained using either the [START_REF] Talbot | Variational principles for inhomogeneous non-linear media[END_REF] approach, or by combining the variational approach of Ponte Castañeda (1991) with the Hashin-Shtrikman linear bound (3.6).

We first consider the case where phase 2 is a perfect conductor. In such case, the contrast is infinite and the bounds only depend on the two parameters (c 1 , n). The case of finite contrast is next addressed in Sect. 5.2.

Perfectly conducting inclusions in a power-law matrix

Consider a composite mixing a power-law matrix (phase 1) with a perfectly conducting medium (phase 2). The energy w 2 in the perfectly conducting phase 2 is given by (3.8).

The exact expression of the linear comparison bound σ LC can be obtained in closed-form [START_REF] Ponte Castañeda | Bounds and estimates for the properties of nonlinear heterogeneous systems[END_REF] and reads as

σ LC σ 1 = (2 -c 1 ) n+1 2 c n 1 .
(5.5)

The corresponding upper bound χ LC on the effective resistivity is given by χ

LC /χ 1 = c 1 /(2 -c 1 ) n+1 2n . Since (W 2 -αV ) * (τ N ) = 0, the bound (5.4) simplifies as σ n + 1 = sup α∈R,τ ∈R τ - √ 3 2 α - c 1 3 (W 1 -αV ) * (τ N ) .
(5.6)

The bound σ in (5.6) is estimated numerically. As mentioned earlier, the most sensitive issue lies in the calculation of the Legendre transform (W 1 -αV ) * (τ N ). Following the developments of Sect. 4, the value of (W 1 -αV ) * (τ N ) is obtained by maximizing the function F + in Eq. (4.37) over R 2 . Because of the growth behavior of w 1 , the search for the supremum of F + can be restricted to a disk of prescribed radius R (see Appendix B), which is helpful for the numerical optimization.

As an illustration, the function F + is plotted in Fig. 5 in the particular example n = 3, σ 1 = 1, τ = 9.2. The symmetry of the function F + with respect to the axis v 3 = 0 is a direct consequence of the expressions (4.36) and (4.37). As can be observed in Fig. 5, the function F + has several local maxima so that special care must be taken when using a gradient algorithm for finding a global maximum: Depending on the initial guess that is used, a gradient algorithm may converge to a local maximum that does not correspond to a global maximum. To overcome such difficulty, the numerical strategy that we have used consists in discretizing the disk of radius R in a large ( 10000) number of sample points, and looking at the sample point with the largest value of F + . That particular sample is next used as the initial guess in a gradient algorithm. For the case depicted in Fig. 5, the global maximum is indicated by the red arrows and is reached at two values of (u 3 , v 3 ). Some contour plots of the function

f (τ, α) = τ - √ 3 2 α -c 1 3 (W 1 -αV ) * (τ N
) are represented in Fig. 6 for a volume fraction c 1 set to 0.7. Maximizing the function f with respect to the two scalar parameters (τ, α) gives the best bound available from the proposed procedure -see Eq. (5.6). The function f is convex in (τ, α) so there is no substantial difficulty in estimating numerically its maximum. In the example considered, the maximum is reached for (τ, α) (17.9, 1.94). The corresponding bound on σ eff is approximatively equal to 5.58.

Influence of the nonlinearity index n

The bound σ/σ 1 is plotted in Fig. 7(left) as a function of the nonlinearity index n, for a volume fraction c 1 = 0.7. The bound σ is obtained by solving (5.6) numerically using the procedure exposed previously. The linear comparison bound σ LC is shown in dashed lines. As can be observed in Fig. 7(left) the bound σ strictly improves on σ LC -except in the special case n = 1 where both σ and σ LC reduce to the linear Hashin-Shtrikman bound (3.6). The relative improvement (σ -σ LC )/σ LC is found to increase with n and reaches 13% for n = 3. In Fig. 7(right) is plotted the corresponding upper bound χ on the effective resistivity χ eff as obtained through the relation χ /χ 1 = (σ 1 /σ) 1/n . In terms of resistivity, the relative improvement on the comparison bound remains small ( (χ LC -χ )/χ LC 4% for n = 3).

The bound k obtained in Sect. 3.2 can be recovered from the asymptotic behavior of χ as n → ∞. In the limit n → ∞, the power-law energy w 1 in (5.1) indeed reduces to the energy (3.7) of a perfect dielectric with unit strength so that the composite behaves as in Sect. 3.2. From the expressions (3.10) and (5.3) of the effective complementary energy, we thus expect that χ eff -→ As can be observed in Fig. 7(right), the convergence of χ towards the asymptotic value k (shown in dotted line) is relatively slow. The linear comparison bound displays a similar behavior.

5.1.2. Influence of 1 influence of the volume fraction c 1 on the bounds is illustrated in Fig. 8(left) in the case n = 3. The results are presented in terms of resistivity. For any given n > 1, the relative improvement (χ LC -χ )/χ LC decreases with the volume fraction c 1 of the matrix, as shown in Fig. 8(right). Values above 5% of the relative improvement (χ LC -χ )/χ LC can be reached as soon as n ≥ 2.5. Correspondingly, the relative improvement (σ -σ LC )/σ LC can exceed 10% for n ≥ 2.5.

Influence of the contrast

We now consider the situation where the nonlinear conductivities 1 and σ 2 take finite, values. Without loss of generality, we can label the phases so that 2 /σ 1 > 1. The bound σ (5.4) is estimated using a numerical procedure similar to that presented in Sect. 5.1. Contrary to the case considered in Sect.5.1, there is no analytical expression of the linear comparison bound σ LC so the latter has to be evaluated numerically as well. The bounds obtained are shown in Fig. 9 functions of the contrast ratio, in the case n = 3.

In a way similar to Fig. 8, the results in Fig. 9 are presented in terms of the resistivity χ = σ -1/n . For the highest contrast considered in the calculations (σ 1 /σ 2 = 100 i.e. χ 2 /χ 1 0.2), the relative improvement (χ LC -χ )/χ LC is in the range 2-3% (for volume fractions c 1 about 0.5). In terms of conductivity, the corresponding relative improvement (σ -σ LC )/σ LC is the range 7-10%. At small contrast (i.e. σ 2 = σ 1 (1 + r) with r 1), the effective conductivity σ eff is known to admit the expansion (Blumenfeld andBergman, 1989, 1991)

σ eff σ 1 = 1 + c 2 r - n + 1 2(n + √ n) c 1 c 2 r 2 + O(r 3 ).
(5.7)

The small-contrast expansion of the linear comparison bound is given by [START_REF] Ponte Castañeda | Bounds and estimates for the properties of nonlinear heterogeneous systems[END_REF])

σ LC σ 1 = 1 + c 2 r - 1 2 c 1 c 2 r 2 + O(r 3 ).
(5.8) can observed on those expressions, the linear comparison is exact only to first-order in contrast (except in the special cases n = 1 and n = ∞). It can be verified that

σ 1 = 1 + c 2 r + O(r 2 ),
i.e. the bound σ exact to first-order in contrast -just as the linear comparison bound σ LC . Although we have not determined the second-order expansion of the bound σ (the calculations proved to be quite cumbersome), the numerical results shown in Fig. 9 suggest that the bound σ is not to second-order in contrast. For instance, in the case n = 3, c 1 = 0.5, the expression (5.7) gives lim r→0

σ eff /σ 1 -1 -c 2 r r 2 = - n + 1 2(n + √ n) c 1 c 2 -0.106.
(5.9)

The numerical results give

r→0 σ/σ 1 -1 -c 2 r r 2 -0.113
does not match with the value in (5.9). As a comparison, the expression (5.8) gives

lim r→0 σ LC /σ 1 -1 -c 2 r r 2 = -0.125.
proposed bound gives a second-order expansion that is closer to (5.9) but a significant gap remains. 0.2 0.4 0.6 0.8 1 0.5 0.6 0.7 0.8 0.9

1 n = 3 χ 2 /χ 1 χ eff /χ 1 c 1 = 0.5 c 1 = 0.7 0.2 0.4 0.6 0.8 1 0 1 2 3 n = 3 χ 2 /χ 1 (χ LC -χ )/χ LC [%]
c 1 = 0.5 c 1 = 0.7

Bounds on the effective complementary energy

The methodology is now adapted to bound the effective complementary energy w * eff from below (which is equivalent to bounding w eff from above). In the isotropic case, w * eff only depends on  through its norm , i.e. we have w * eff () = φ * eff () for some function φ * eff . Setting J = (j 1 , • • • , j n ), we consider the extended complementary energy

W * (J , x) = n i=1 w * (j i , x). Let J = ( 1 , • • • , n ) be such that 1 = • • • = n = . Reproducing the arguments of Sect.
leads to the lower bound

φ * eff () ≥ 1 n J • T + Ũ ( J ) - Ω (W * -Ũ ) * , x)dω (6.1) where Ũ satisfies 2 Ω Ũ (J )dω ≥ Ũ ( J ) for any J ∈ S( J ) (6.2) with S( J ) = {(j 1 , • • • , j n ) : div j i = 0 in Ω, j i (x) • n = i • n on ∂Ω}.
For two-dimensional problems, potentials Ũ that satisfy (6.2) can easily be generated from quasiconvex potential U , i.e. potentials that satisfy (2.12). For any given potential U defined in R 2n , considered indeed the 'rotated' potential U ⊥ defined by

U ⊥ (J ) = U (r ⊥ J )
where r ⊥ J = (r ⊥ j 1 , • • • , r ⊥ j n ) and r ⊥ is the two-dimensional rotation with angle π/2. For two-dimensional fields, we have

j ∈ S() ⇐⇒ r ⊥ j ∈ K(r ⊥ ). (6.3)
us briefly justify the property (6.3). It is well-known (and easy to verify) that r maps divergence-free fields to curl-free fields. For any j ∈ S(), the rotated field r ⊥ j can thus be written as r ⊥ j = ∇u for function u. From the boundary condition j • n =  • n, we obtain that (∇u -∧ n = 0 on ∂Ω, i.e. that u = (r ⊥ ) • x on ∂Ω (up to some constant that can be taken as 0). This shows that r ⊥ j ∈ K(r ⊥ ). We note that the relation (6.3) also holds for extended fields, i.e. J ∈ S( J ) ⇐⇒ r ⊥ J ∈ K(r J ).

(6.4)

A notable consequence of (6.4) is that U ⊥ satisfies (6.2) whenever U is quasiconvex. For any J ∈ S( J ), we have indeed

Ω U ⊥ (J )dω = Ω U (r ⊥ J )dω ≥ U (r ⊥ J ) = U ⊥ ( J )
the inequality comes from the quasiconvexity of U together with the fact that r ⊥ J ∈ K(r ⊥ J ). now specialize the above relations to n = 3 and to quasiconvex potentials U of the (2.19). Since the two-dimensional determinant is invariant by r ⊥ , the quasiconvex function V in (2.20) satisfies V ⊥ = V . Taking J = N and Ũ = αV in (6.1), we obtain the lower bound on the complementary effective energy:

φ * eff () ≥ sup α∈R,τ ∈R τ - √ 3 2 α 2 - 1 3 Ω (W * -αV ) * (τ N , x) dω .
(6.5) that (6.5) is formally identical to (2.26). Consequently, all the developments of Sect. 4 remain valid by making the formal substitutions W → W * and E → J .

Perfectly insulating inclusions in a power-law matrix

As an illustration, consider a power-law matrix (phase 1) with perfectly insulating inclusions (phase 2). The energy w 1 is of the form (5.1). The energy w 2 (e) is null for all e. By Legendre duality, the complementary energy function w * 2 is given by

w * 2 (j) = 0 if j = 0, +∞ otherwise;
The effective complementary energy w * eff is of the form (5.3). The formula (6.5) gives the inequality χ ≤ χ where χ eff is the effective resistivity and χ is given by χ

m 1 = sup α∈R,τ ∈R τ - √ 3 2 α - c 1 3 (W * 1 -αV ) * (τ N ) .
The bound χ is non trivial m > 1, i.e. for n < 1. Since

χ eff = σ -1/n
eff , the lower bound χ on the effective resistivity translates as an upper bound on the effective conductivity eff . We have indeed σ eff ≤ χ-n . (6.6)

The upper bound given by (6.6) is represented in Fig. 10 along with the linear comparison bound σ LC , given by [START_REF] Ponte Castañeda | Bounds and estimates for the properties of nonlinear heterogeneous systems[END_REF])

σ LC σ 1 = c 1 (2 -c 1 ) (n+1) 2
.

For n = 0.1, the relative improvement on the linear comparison bound σ LC is about 4%.

Concluding remarks

Nonlinear Hashin-Shtrikman type bounds can be obtained by combining the translation method with the idea of embedding, whereby n copies of the original problem are considered in an extended problem of dimension 2n. In contrast with the linear case, the bounds obtained by taking 2 are generally not optimal. In this paper, we investigate the choice n 3, thus considering an extended problem of dimension 6.

the original problem in R 6 to a genuine improvement of the bounds, but it comes at a price: The direct evaluation of the bounds is difficult to perform, even using numerical techniques. For isotropic constituents, invariance properties lead to considerable simplifications, thus making up for the increase in dimensionality that comes with the embedding approach. Ultimately those simplifications make for an efficient and reliable way of evaluating the bounds.

This paper concentrates on two-dimensional conductivity, but it should be mentioned that the results obtained apply to the mathematically equivalent problem of bounding the elastic energy of fiber-reinforced composites under transverse shear. is tempting to further extend the presented approach to three-dimensional conductivity or two-dimensional elasticity. For instance, in the case of three-dimensional conductivity, it can be expected copies of the original problem are necessary to improve on the linear comparison bound, thus leading to an extended problem of dimension 12. The issue of evaluating the corresponding bounds becomes even more critical, but the ideas used in this paper could be possibly be helpful to make the calculations tractable.

The problem (A.1) is easily solved for τ = 0: Since the largest eigenvalue of M is equal to √ 3, the supremum in (A.1) is attained for any eigenvector of M (corresponding to the eigenvalue √ 3) of norm 3 and we have g * (0) = 3 √ 3/2. We note that such an eigenvector if given by 1 2 (1, -1, √ 3, -2, 0).

In the following, we address the general situation where τ = 0. Let E attaining the supremum in (A.1). We write E as E = (u 1 , v 1 2 , v 2 , u 3 , v 3 ) with u 2 i + v 2 i ≤ 1 (i = 1, 2, 3). The stationarity conditions give

τ /2 + v 2 + v 3 = λ 1 u 1 - √ 3τ /2 -u 2 -u 3 = λ 1 v 1 -τ /2 -v 1 + v 3 = λ 2 u 2 - √ 3τ /2 + u 1 -u 3 = λ 2 v 2 -τ -v 1 -v 2 = λ 3 u 3 u 1 + u 2 = λ 3 v 3 (A.3)
where the Lagrange multipliers λ i satisfy

λ i ≥ 0, λ i (u 2 i + v 2 i -1) = 0. (A.4)
The second-order optimality condition reads as .5) for all δE = (δe 1 , δe 2 , δe 3 ) verifying (u i , v i ) • δe i ≤ 0 if u 2 i + v 2 i = 1. Using (A.3) and (A.5), it can be proved that all the Lagrange multipliers λ i are non zero, i.e. that u 2 i +v 2 i = 1 for i = 1, 2, The main argument is that if one of the Lagrange multipliers is equal to 0 then one can construct an eigenvector δE such that M • δE = √ 3 δE and λ i (u i , v i ) • δe i = 0 for i = 1, 2, which is in contradiction with (A.5).

0 ≥ 1 2 δE • M • δE + 3 i=1 λ i (u i v i ) • δe i (A
We now prove that, up to a rotation R i in (4.2), E is of the form E = (cos θ, sin θ, -cos θ, sin θ, , 0) (A.6) for some in {-1, 1}. First assume that λ 1 + λ 2 + λ 3 -λ 1 λ 2 λ 3 = 0. In that case, the linear (A.3) can be inverted to give

        u 1 v 1 u 2 v 2 u 3 v 3         = 1 2 τ λ 1 + λ 2 + λ 3 -λ 1 λ 2 λ 3        Setting x = z + /2 + τ /2
√ 3, (A.13) can be rewritten as 0 = 8z 3 -6z + √ 3τ 3 9 + 2 .

It is necessary at this point to distinguish between the cases = 1 and = -1.

For = -1, the cubic equation (A.13) admits 3 solutions (denoted by x 3 > x 2 > x 1 ) if τ < , 2 solutions (denoted by x 2 > x 1 ) if τ = τ 1 , 1 (denoted by x 1 ) if τ > τ 1 ;

where τ 1 = √ 3 2 2/3 (see Fig. A.11). Let θ be value of θ satisfying x i = cos θ i along with (A.11). We need to find which of the values g -(θ 0 ), g -(θ 1 ), g -(θ 2 ), g -(θ 3 ) is the largest. Setting τ 0 = 1 + √ 3, some manipulation shows that g -= g -(θ 3 ) for 0 < τ ≤ τ 0 , g -(θ 0 ) for τ 0 ≤ τ. (A.14) For = 1, (A.13) only has one solution (denoted by x 1 ) in R. Calculating the second derivative g + shows that x 0 = -/2 corresponds to a local minimum and that x 1 to local maximum of g + . It follows that sup g + = g + 1 )

where θ 1 satisfies x 1 = cos θ 1 along with (A.11). Comparing g + (θ 1 ) with max(g -(θ 0 ), g -(θ 3 ))

shows that sup g + ≤ sup g -. Hence g * (τ ) is equal to sup g -as given by (A.14).

Appendix B. Additional details for a power-law behaviour

For a power-law energy w(e) = σ n+1 e n+1 , the function F introduced in Eq. (4.24) specializes as Combining two preceding relations, we get

F(E) = T • E + 1 2 E • M • E - σ n + 1
3 i=1 e n+1 i ≥ 3 -n+1 2 E n+1 .
Substituting in (B.3) yields

F(E) ≤ T E - 1 2 A( E ) E 2 .
where

A( E ) = 3 -n+1 2 2σ n + 1 E n-1 - √ 3.
Since σ > 0 and n > 1, there exists R 0 > 0 such that A( E ≥ 1 for any E ≥ 0 . Such a value of R 0 is by

R 0 = 3 n+1 2(n-1) (n + 1)(1 + √ 3) 2σ 1 n-1
. any E verifying E ≥ 0 , we obtain

F(E) ≤ T E - 1 2 E 2 = 1 2 E (2 T -E ).
Hence F(E) ≤ 0 for any E such that E ≥ max(R 0 , 2 T ). A value of R satisfying (B.2) is thus given by R = max(R 0 , 2 T ).

  Consider a two-dimensional inhomogeneous electric conductor occupying a domain Ω of unit volume. The electric field e and the current density j are related by the local constitutive

Figure 1 :

 1 Figure 1: The three directions (ē1, ē2, ē3) constituting the vector Ē used in the calculation of the bounds.

Figure 2 :

 2 Figure 2: Plots of the function g * and of its subdifferential ∂g * . The function g * is depicted in Fig. 2(left). It follows from (3.16)-(3.18) that g * is differentiable for τ = τ 0 , with a derivative given by dg * dτ (τ ) = cos θ -√ 3 sin θ + 1 for 0 ≤ τ < τ 0 , 3 for τ > τ 0 . (3.19)

  27), cos θ * and sin θ * are obtained by setting τ = τ * in (3.17). Using (3.27) and the expression (3.16) of g * , Eq. (3.26) simplifies as 0 cos θ * + 1) sin θ * . (3.28) Let a = 1 + cos θ * and b = -√ 3 sin θ * . From (3.27) and (3.28) we have

Figure 3 :

 3 Figure 3: Upper bounds on the effective strength k eff of a two-phase dielectric composite: Wiener bound kW (dotted line), linear comparison bound kLC (red dashed curve), proposed bound k (solid curve).

Figure 4 :

 4 Figure4: Action of G on a vector (e1, e2, e3) (shown in upper left). The vector (e 1 , e 2 , e 3 ) shows the transformed of (e1, e2, e3) for each transformation in G. The transformations Ri involve some rotations (with angle in {±π/3, ±2π/3}) and some permutations of the vectors (e1, e2, e3). The transformations Si involve some reflections (through axis shown as dashed lines) and some permutations of the vectors (e1, e2, e3).

  .11) The definition (2.9) shows that (W -αV ) * (T , x) is the supremum of linear functions in T . Hence (W -αV ) * (T , x) in convex in T . Consequently, F is convex and the function T -→ ē 3 N • T -F (T ) that appears in (4.10) is concave. It follows that the optimal values of T in (4.10) are fully characterized by the relation 1 3 ēN ∈ ∂F (T ).(4.12)

Figure 5 :

 5 Figure 5: Plot of F+(u3, v3) for n = 3, σ1 = 1, τ = 9.2.

Figure 6 :

 6 Figure 6: Contour plots of f (τ, α) for c1 = 0.7, n = 3, σ1 = 1.

Figure 7 :

 7 Figure 7: Bounds on the effective properties of a power-law composite, as a function of the nonlinearity index n. Case of infinite contrast, with c1 = 0.7. Solid line: proposed bound. Dashed line: linear comparison bound.

Figure 8 :

 8 Figure 8: Left: Upper bounds on the effective resistivity of a power-law composite, as a function of the volume fraction c1. Case of infinite contrast, with n = 3. Solid line: proposed bound. Dashed line: linear comparison bound. Right: relative improvement as a function of the volume fraction c1 for n = 2, 3, 5.

Figure 9 :

 9 Figure 9: Left: Upper bounds on the effective resistivity of a power-law composite, as a function of the contrast ratio. Solid line: proposed bound. Dashed line: linear comparison bound. Right: relative improvement as a function of the contrast ratio.

Figure 10 :

 10 Figure 10: Upper bounds on the effective conductivity of a power-law composite, as a function of the nonlinearity index n. Case of infinite contrast, with c1 = Solid line: proposed bound. Dashed line: linear comparison bound.

FFigureF

  Figure Stationary points of g-.Since M is symmetric with eigenvalues 0 and ± √ 3, we haveE • M • ≤ 3 E 2 for any E. Therefore

  21)Since P is continuous in (k, τ ), the function sup 0≤τ ≤τ 0 P (k, τ ) is continuous in k.

				In such
	condition, (3.21) implies that
	sup	P ( k, τ ) = 0.	(3.22)
	0≤τ ≤τ 0
	Noting that τ → P ( k, τ ) is concave, the values τ * reaching the supremum in (3.22) are
	fully characterized by the relation 0 ∈ ∂P ( k, τ ), i.e.
	3 c 1	k ∈ ∂g * (τ * ).	(3.23)
	As a consequence of (3.19) and (3.20), Eq. (3.23) can be rewritten as

except in the unusual case where the constitutive phases are governed by non strongly convex energyfunctions[START_REF] Ponte Castañeda | The effective behavior of nonlinear composites: a comparison between two methods[END_REF].

distinctive property (6.2) enters into the framework of A-quasiconvexity[START_REF] Fonseca | A-quasiconvexity, lower semicontinuity, and young measures[END_REF] and can be interpreted as quasiconvexity in divergence-free fields. To avoid any confusion, use of the word 'quasiconvexity' is restricted in this paper to the standard definition (2.12).

( √ 3 -λ 2 ) √ 3λ 2 λ 3 -2λ 2 -λ 3 -λ 3 ( √ 3 -λ 1 ) √ 3λ 1 λ 3 -2λ 1 -λ 3 2λ 1 λ 2 -√ 3λ 1 -√ 3λ 2 λ 1 -λ 2

Appendix A. Calculation of g * for dielectric breakdown

In this Appendix we address the calculation of the function g * in (3.14). We have

where C = {E = (e 1 , e 2 , e 3 ) : e i ≤ 1} and

Using the expressions (A.7), the relation

where P is the cubic polynomial defined by

Calculating the discriminant of P shows that P only has one root R. we deduce from (A.8) and (A.10) that two of the multipliers λ i are necessarily equal. By applying one of the transformation R i in (4.2), we can assume that λ 1 = λ 2 . The expressions (A.7) rightly show that

The above analysis needs to be modified if λ 1 +λ 2 +λ 3 -λ 1 λ 2 λ 3 = 0. In that case, some manipulations on (A.3) (that we do not report here) show that λ 1 , λ 2 , λ 3 are necessarily all equal to √ 3 and that E parallel to N (and therefore of the form (A.6)).

The developments so far show that we only need to optimize the function in (A.2) respect to E of the form (A.6). For such E, a direct calculation gives

Hence we are left with solving the scalar optimization problems sup θ g -(θ) and sup θ g + (θ) with g + (θ) = 2(cos θ -1) sin θ + τ (cos θ -√ 3 sin θ -1), g -(θ) = 2(cos θ + 1) sin θ + τ (cos θ -√ 3 sin θ + 1).

The stationarity condition g ± (θ) = 0 gives .11) with = 1 for g + and = -1 for g -. Taking the square of (A.11) leads to the polynomial equation

where x = cos θ. Observe from (A.11) that 0 = -/2 is a solution of (A.12). Factorizing the polynomial in (A.12) by (x + /2), we obtain that the solutions of (A.12) (distinct from -/2) are the roots of the cubic equation .13)


Fatima El Salahdine

Naima El Kaabouch
email: naima.kaabouch@engr.und.edu

Hassan El Ghazi
email: elghazi@inpt.ac.ma

Fatima Salahdine
email: fatima.salahdine@und.edu

Bayesian Compressive Sensing with Circulant Matrix for Spectrum Sensing in Cognitive Radio Networks

Keywords: Cognitive radio networks, Wideband spectrum sensing, Compressive sensing, Bayesian models, Circulant matrices, 1 norm minimization

    Bayesian compressive sensing with circulant matrix for spectrum sensing in cognitive radio networks

I. INTRODUCTION

Spectrum sensing is one of the main processes of the cognitive radio cycle [START_REF] Kaabouch | Software-Defined and Cognitive Radio Technologies for Dynamic Spectrum Management[END_REF][START_REF] Manesh | A Comparative Study of Spectrum Sensing Techniques for Cognitive Radio Systems[END_REF][START_REF] Reyes | A Spectrum Sensing Technique Based on Autocorrelation and Euclidean Distance and its Comparison with Energy Detection for Cognitive Radio Networks[END_REF][START_REF] Subramaniam | Spectrum occupancy measurement: an autocorrelation based scanning technique using USRP[END_REF][START_REF] Salahdine | Matched filter detection with dynamic threshold for cognitive radio networks[END_REF]. In order to sense the wideband radio spectrum, communication systems must use multiple RF frontends simultaneously, which results in long delays, high hardware cost, and computational complexity [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF]. To address these problems, fast and efficient sensing techniques are needed. Compressive sensing has been proposed as a low cost solution to speed up the scanning process and reduce the computational complexity. It involves three main processes: sparse representation, encoding, and decoding. During the first process, the signal, S, is projected in a sparse basis. During the second process, S is multiplied by a sampling matrix, , of MxN elements to extract M samples from N of the signal, S, where M << N. In the last process, the signal is reconstructed from the few M measurements [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF][START_REF] Salahdine | A survey on compressive sensing techniques for cognitive radio networks[END_REF][START_REF] Sharma | Application of compressive sensing in cognitive radio communications: A survey[END_REF][START_REF] Donoho | Compressed sensing[END_REF].

For the encoding process, a number of sampling matrices have been proposed in the literature, including random matrix [START_REF] Abo-Zahhad | Compressive sensing algorithms for signal processing applications : A Survey[END_REF][START_REF] Rudelson | Sparse reconstruction by convex relaxation: Fourier and Gaussian measurements[END_REF], Circulant matrix [START_REF] Valsesia | Compressive signal processing with Circulant sensing matrices[END_REF][START_REF] Yin | Practical compressive sensing with Toeplitz and Circulant matrices[END_REF], Toeplitz matrix [START_REF] Yin | Practical compressive sensing with Toeplitz and Circulant matrices[END_REF], and deterministic matrix [START_REF] Nguyen | Deterministic sensing matrices in compressive sensing: A survey[END_REF]. Because of their simplicity, more interest has been paid to random matrices. These matrices are randomly generated with independent and identically distributed (i.i.d) elements such as Gaussian and Bernoulli distributions [START_REF] Donoho | Compressed sensing[END_REF][START_REF] Abo-Zahhad | Compressive sensing algorithms for signal processing applications : A Survey[END_REF]. In general, compressive sensing requires that the sampling matrix satisfies the Restrict Isometry Property (RIP) condition [START_REF] Bubacarr | Improved bounds on restricted isometry constants for Gaussian matrices[END_REF]. RIP is a characteristic of orthonormal matrices bounded with a Restrict Isometry Constant (RIC), which is a positive number between 0 and 1 that respects the RIP condition [START_REF] Cai | New bounds for restricted isometry constants[END_REF]. This condition guarantees the uniqueness of the reconstructed solution, , during the decoding process. For random matrices, the matrix satisfies the RIP condition for small RIC [START_REF] Abo-Zahhad | Compressive sensing algorithms for signal processing applications : A Survey[END_REF][START_REF] Rudelson | Sparse reconstruction by convex relaxation: Fourier and Gaussian measurements[END_REF]. However, these matrices require a great deal of processing time and high memory capacity to store the matrix coefficients [START_REF] Valsesia | Compressive signal processing with Circulant sensing matrices[END_REF][START_REF] Yin | Practical compressive sensing with Toeplitz and Circulant matrices[END_REF]. Because of the randomness, the results are uncertain, which makes the signal reconstruction inefficient.

Unlike random matrices, Circulant matrices are efficient, fast in terms of signal acquisition, and require fewer measurements [START_REF] Valsesia | Compressive signal processing with Circulant sensing matrices[END_REF]. A Circulant matrix is a structured matrix determined using a predefined vector by a cyclic permutation [START_REF] Yin | Practical compressive sensing with Toeplitz and Circulant matrices[END_REF]. This matrix satisfies the RIP condition for a small number of measurements [START_REF] Rauhut | Restricted isometries for partial random Circulant matrices[END_REF]. Unlike random matrices, Circulant matrices are not universal. Universality means the sampling matrix can be used to compress a signal sparse in any domain. In addition, Circulant matrices have been used only with the ℒ 1 norm minimization technique [START_REF] Valsesia | Compressive signal processing with Circulant sensing matrices[END_REF][START_REF] Yin | Practical compressive sensing with Toeplitz and Circulant matrices[END_REF][START_REF] Abo-Zahhad | Compressive sensing algorithms for signal processing applications : A survey[END_REF].

For the decoding process, a number of algorithms that exploit the sparsity feature of signals have been proposed in the literature [START_REF] Abo-Zahhad | Compressive sensing algorithms for signal processing applications : A survey[END_REF][START_REF] Xu | An improved Toeplitz measurement matrix for compressive sensing[END_REF][START_REF] Tropp | Signal recovery from random measurements via orthogonal matching pursuit[END_REF][START_REF] Ji | Bayesian compressive sensing[END_REF][START_REF] Babacan | Bayesian compressive sensing using Laplace priors[END_REF][START_REF] Baron | Bayesian compressive sensing via belief propagation[END_REF][START_REF] Needell | Stochastic gradient descent, weighted sampling, and the randomized Kaczmarz algorithm[END_REF][START_REF] Blumensath | Iterative hard thresholding for compressed sensing[END_REF][START_REF] Narayanan | Greedy Pursuits Assisted Basis Pursuit for Compressive Sensing[END_REF][START_REF] Tropp | Signal recovery from random measurements via orthogonal matching pursuit[END_REF][START_REF] De | Sparse representation for blind spectrum Sensing in cognitive radio: A compressed sensing approach[END_REF]. A sparse signal can be estimated from a few measurements by solving the underdetermined system using three different types of algorithms: Iterative relaxation [START_REF] Abo-Zahhad | Compressive sensing algorithms for signal processing applications : A survey[END_REF], Greedy [START_REF] Xu | An improved Toeplitz measurement matrix for compressive sensing[END_REF], and Bayesian models [START_REF] Ji | Bayesian compressive sensing[END_REF]. The iterative relaxation category includes techniques that solve the underdetermined system using linear programing. Some techniques classified under this category are ℒ 1 norm minimization known as basis pursuit [START_REF] Abo-Zahhad | Compressive sensing algorithms for signal processing applications : A survey[END_REF], gradient descent [START_REF] Needell | Stochastic gradient descent, weighted sampling, and the randomized Kaczmarz algorithm[END_REF], and iterative thresholding [START_REF] Blumensath | Iterative hard thresholding for compressed sensing[END_REF]. Greedy algorithms consist of selecting a local optimal at each step in order to find the global optimum, which corresponds to the estimated signal coefficients. Examples of techniques classified under this category are matching pursuit [START_REF] Narayanan | Greedy Pursuits Assisted Basis Pursuit for Compressive Sensing[END_REF], orthogonal matching pursuit [START_REF] Tropp | Signal recovery from random measurements via orthogonal matching pursuit[END_REF][START_REF] De | Sparse representation for blind spectrum Sensing in cognitive radio: A compressed sensing approach[END_REF], and stage wise orthogonal matching pursuit [START_REF] Donoho | Sparse solution of underdetermined systems of linear Equations by stagewise orthogonal matching pursuit[END_REF]. Bayesian compressive sensing algorithms consist of using a Bayesian model to estimate the unknown parameters in order 978-1-5090-1496-5/16/$31.00 ©2016 IEEE to deal with uncertainty in measurements. Examples of techniques classified under this category are: Bayesian model using relevance vector machine learning [START_REF] Abo-Zahhad | Compressive sensing algorithms for signal processing applications : A Survey[END_REF], Bayesian model using Laplace priors [START_REF] Ji | Bayesian compressive sensing[END_REF], and Bayesian model via belief propagation [START_REF] Babacan | Bayesian compressive sensing using Laplace priors[END_REF]. All these Bayesian based algorithms were used only with random matrices.

Iterative relaxation algorithms are more accurate compared to Greedy algorithms, but they are complex, uncertain, require high number of measurements, and, thus high processing time. On the other hand, Greedy algorithms are fast and require low processing time; however, they are inefficient, uncertain, and require more measurements for the reconstruction process. Bayesian based techniques combine the strengths of both categories. They are fast, accurate, require less measurements for a high recovery rate, and can deal with uncertainty. In this paper, we propose an approach that combines the strengths of both Circulant matrices and Bayesian models to address the previously mentioned problems during the encoding and decoding processes.

The paper is organized as follows. Section II presents the methodology followed for both encoding and decoding processes as well as the performance evaluation. Section III discusses the simulation results of the proposed approach based on specific metrics. Finally, a conclusion is given at the end.

II. METHODOLOGY

A. Bayesian compressive sensing

As previously explained, compressive sensing involves three processes: sparse representation, sensing matrix, and reconstruction, as shown in Fig. 1. In our approach, for the sparsity representation process, we assumed the signal to be sparse. For the sensing matrix process, we used the Circulant matrix and for the reconstruction process, we used the Bayesian model. The Circulant matrix, C, is generated from a vector c, where

( , ) = ( ) ( ) for i, j = 0…N. It is expressed as = … ⋯ ⋯ ⋮ ⋮ ⋮ ⋮ ⋮ … (1)
Where c is a vector given as ( ,,… , ).

The values of c are chosen randomly according to a suitable probability distribution to reduce the amount of randomness of the sensing matrix. During the encoding process, each column of the matrix C is obtained by right cyclic shifting the previous column. Then, a MxN partial Circulant matrix, M c , is defined as the submatrix of C and considered for the signal sampling, with M<<N. The signal, S, is then multiplied by for signal compression. This multiplication is fast because of the reduced number of random coefficients in the Circulant matrix [START_REF] Rudelson | Sparse reconstruction by convex relaxation: Fourier and Gaussian measurements[END_REF][START_REF] Valsesia | Compressive signal processing with Circulant sensing matrices[END_REF].

For the reconstruction process, we used the Bayesian model, which is a probabilistic approach that requires a prior knowledge of parameters to calculate the posterior distribution of the unknown parameters. The Bayesian compressive sensing process consists of finding the sparse solution of a regression problem by exploiting the probabilistic distributions. It solves the underdetermined system and finds the accurate solution by estimating efficiently the unknown parameters using the information that we have about the system. It is based on two main elements: the knowledge about the linear relationship between the signal measurements and the original signal, and the knowledge about the fact that the original signal is k-sparse. Under the Bayesian model, the k-sparse signal is acquired through a product with the Circulant matrix. A noise, W, is added to the signal measurements, which includes the noise measurements and the sparse representation error. Fig. 2 illustrates the simulation methodology of our proposed model. The noisy measurements can be formulated as

= + (2)
According to the Theorem of Central Limit for N >> M, W can be approximated as a zero mean Gaussian noise with unknown variance δ w , which can be expressed as N (0, δ w ). The signal to be approximated can be considered as a Gaussian variable with S = (S 1 , S 2 ,… , S N ). Therefore, the Bayesian model implies that the noisy measurements, R, is an i.i.d Gaussian and depends on the unknown S and δ w . This is expressed as

| , ~ ( , ) (3)
Fig. 3 illustrates the proposed Bayesian model in which the unknown signal, noise, and vector to generate the Circulant sampling matrix are parents of the noisy measurements. Noise variance δ w , signal mean μ S , and signal variance δ s are the parameters of the noise and the signal that need to be estimated. The model specifies the conditional probabilities of the measurements P(R/S) and the noise ( | ). P(R/S) and P( | ) present the probability densities of R, as functions of the values taken by the signal and the noise variance respectively. The conditional probability of the signal to be estimated given the measurements can be expressed through the Bayes' rule as

( / ) = ( / ) ( ) ∑ ( / ′) ( ′) (4)
where S' represents the other alternative solution of the underdetermined system. Consequently, the Gaussian likelihood of the noisy measurements can be expressed from the previous conditional probability as

P( | , ) = (2π ) exp( - ‖ - ‖ ) (5)
Given a vector c, the Circulant sampling matrix coefficients are generated to be used for the Bayesian reconstruction process. This matrix should be the same as the one used for compression to reduce the randomness. With prior knowledge that S is sparse and is known, S and δ w are the two quantities to be estimated. The prior density of the unknown signal P(S/k) in terms of sparsity can be expressed as

( / ) = ( /2) ( ∑ ) (6)
where k is a parameter that represents the sparsity of S = (S 1 , S 2 ,…,S N ). Using the Bayesian model, the underdetermined system in (2) becomes a linear problem with S sparse and (2) can be reformulated as

= arg min ‖ - ‖ + ‖ ‖ ( 7
)
where z is a positive scalar. The objective of our Bayesian model is to find the posterior probabilistic distribution for S and δ w taking into account the known evidences. The maximum posterior probability corresponds to the sparsest solution of the undetermined system presented in equation [START_REF] Salahdine | A survey on compressive sensing techniques for cognitive radio networks[END_REF].

The algorithm calculates the joint probability distribution of all unknown parameters and computes the prior distribution of each element of S with the hyper parameters a and b. The hyper parameter a = (a 1 , a 2 ,…, a N ) represents the initial posterior of the signal variance and the hyper parameter b represents the initial posterior of the noise variance. The prior distribution of S given the hyper parameters a and b can be expressed as the product of the conjugate prior of signal variance ( / , ) and the likelihood function of S i, which is defined as a zero mean Gaussian prior for each signal coefficients ( /0, ). It is also called the marginal likelihood for Bayes estimation [START_REF] Cai | New bounds for restricted isometry constants[END_REF]. This probability of the signal given a and b is expressed as

( / , ) = ∏ ( /0, ) ( / , ) (8)
The algorithm optimizes the hyper parameters for the Gaussian process in an iterative loop, estimates new values of a and b, and then maximizes the marginal likelihood using the new estimated values of a and b. The algorithm is based on the previous results for learning and searching for the new values of the hyper parameters a and b. Taking into account the assumption about the knowledge of a and b in addition to and R, the posterior probabilistic distribution of S can be expressed as a Gaussian distribution S ~ N (μ S , δ S ) with a mean μ S and variance δ S which are given by

= = ( + ) (9)
where A = diag(a 1 , a 2 , …, a N ). The last estimated value of b represents the noise variance. At the end of the algorithm, the signal is approximated and the uncertainty is reduced.

B. Performance evaluation

In order to evaluate the efficiency of our approach, the results of this model have been compared to the results of the basis pursuit technique [START_REF] Yin | Practical compressive sensing with Toeplitz and Circulant matrices[END_REF] using several metrics. These metrics are: Mean Square Error, reconstruction error, correlation coefficients, processing time, recovery time, and sampling time.

The reconstruction error is a metric that calculates the norm of the difference between the expected signal and the original signal divided by the norm of the original signal. It is expressed as

= - ‖ ‖ (10)
Mean Square Error is a metric that measures the average magnitude of the squared difference between the reconstructed signal and the original signal. It corresponds to one of the loss functions used for error estimation. It is given by

= ∑ ( -) (11)
Correlation measures the similarity between the original signal, S, and the reconstructed signal, . The measure of correlation is known as the correlation coefficient C c , which is a scalar quantity that takes values between -1 and 1. It is expressed as

= ∑( ) -(∑ ) ( ∑ ) (∑ ) -(∑ ) (∑ ) -∑ (12)
When C c is positive and less than 1, it means the two signals are positively correlated and the strength of the correlation is expressed with a percentage value. When C c is null, it means there is no relationship between the two signals. When C c is negative and greater than -1, it means the two signals are negatively correlated and the strength of the correlation is expressed with a percentage value.

Recovery time is the time required by the reconstruction process to reconstruct the signal. Sampling time is the time required by the sampling matrix process in order to compress the signal using a specific matrix. Finally, the processing time is the time required to perform all processes.

III. RESULTS AND DISCUSSION

The two algorithms, Bayesian compressive sensing and ℒ1 norm minimization with the Circulant matrix, were implemented and extensively tested. Their efficiencies were compared using the metrics previously mentioned (R e, MSE, C c , t p , t r and t s ). In this performance evaluation, we investigated the efficiency of the Circulant matrix in sampling signals and compared its results with those of random matrices. We also investigated the performance of our Bayesian model and compared its results with the basis pursuit technique.

Examples of the results are shown in Figs. 4 to 7. Fig. 4(a) shows an example of original signals with 200 samples and 15 spikes. The noise was added to the original signal and fed to the two algorithms. Fig. 4(b) shows the output signal after applying the Circulant sampling matrix and the Bayesian technique to the signal with added noise. Fig. 4(c) represents the output signal after applying random sampling and basis pursuit technique to the original signal with added noise. As one can see in Fig. 4(c), the output signal has more fluctuations than the output signal shown in Fig. 4(b). These fluctuations correspond to the null coefficients that are not reconstructed as zero coefficients. Thus, the reconstruction with the Circulant matrix is more efficient compared to the reconstruction with random matrix. The sampling time, t s , was computed for each technique. The results show that random sampling matrix requires a great deal of time to process compared to the Circulant matrix. For example, for the signal seen in Fig. 4(a), the t s of Circulant matrix is 0.06 ms while the t s of random matrix is 0.30 ms. This result shows that dense matrices are slow in terms of computation because of the high required number of measurements and the randomness of their coefficients.

Figs. 5(b) and (c) show the reconstructed signal using basis pursuit and Bayesian model with Circulant matrix technique, respectively. As one can see, for the Bayesian model, the output signal is similar to the original signal and the spikes are completely recovered. However, the output signal of basis pursuit presents more fluctuations than the output signal of the Bayesian model, as shown in Fig. 5(b). All the results show that the Bayesian reconstruction is more efficient than basis pursuit reconstruction. In addition, the sparsity level of the output signal is 14 for the Bayesian technique and 200 for the basis pursuit technique. The basis pursuit technique cannot estimate the exact value of each coefficient of the original signal, and it estimates the zero values as non-zero values with low magnitude. Moreover, the number of measurements to recover the signal needed by each technique is 15 for the Bayesian technique while it is 200 for the basis pursuit technique. Thus, the Bayesian technique is more efficient in reconstructing the original signal and it also requires fewer measurements than the basis pursuit algorithm. For the other metrics, Table I gives an example of results of the comparison performance. As shown in this table, for the reconstruction error R e , the Bayesian technique has an average of 0.77% of reconstruction error. However, the basis pursuit has 6.7% of reconstruction error. Thus, Bayesian is 80 times more efficient than the basis pursuit with the same matrix. Unlike the basis pursuit, Bayesian technique permits to reconstruct the signal with a very small number of errors, which can be explained by the fact that the technique is able to deal with the uncertainty. Basis pursuit reconstructs the signal with high error level, which can be explained by the fact that this technique cannot handle uncertainty due to the noisy measurements.

This table also shows that for the correlation metric, the Bayesian technique presents an average of 100% of correlation while the basis pursuit technique presents an average of 82.87%. Thus, both techniques present a high correlation with values close to 100%, which indicates that the two signals are positively correlated. However, the Bayesian technique presents a better correlation coefficient.

For the recovery time, the Bayesian technique requires an average of 0.90 ms to recover the original signal, but the basis pursuit technique requires 7.20 ms. This represents 12 times the recovery time of the Bayesian technique, thus slower. For the processing time, Bayesian technique requires an average of 0.96 ms to process while the basis pursuit requires an average of 7. 26 ms. Thus, the Bayesian is 7 times faster than the basis pursuit. These examples of results show that the Bayesian technique with Circulant matrix is more accurate, faster, and deals with uncertainty. In addition, our model requires less measurements, less sampling time, less recovery time, and less processing time. It also estimates the original signal with a low sparsity level, high correlation, and low Mean Square Error as well as handles uncertainty during the encoding and decoding processes. Thus, our proposed approach includes the strengths of both Bayesian reconstruction and Circulant sampling matrix.

CONCLUSIONS

In this paper, we have proposed a Bayesian approach that combines the Circulant sampling matrix with the Bayesian model. This approach reduces the randomness and deals with uncertainty during the compressive sensing processes. It also speeds up the spectrum sensing and scanning of the radio spectrum in cognitive radio networks. The results of the proposed approach are discussed and compared to those of basis pursuit with Circulant and random matrix techniques. The performance comparison involved several metrics that evaluate most aspects of the characteristics of compressive sensing approaches, including reconstruction success, speed, robustness, efficiency, memory, and certainty. The results show that the Bayesian based algorithm with Circulant matrix is more efficient and faster in compressing and recovering signals than the Bayesian with random matrix as well as the basis pursuit with Circulant and random matrices.

Fig. 1 .

 1 Fig. 1. Block diagram of compressive sensing model.

Fig. 2 .

 2 Fig. 2. The model of the Bayesian compressive sensing with Circulant matrix.

Fig. 3 .

 3 Fig. 3. Graphical model of the Bayesian compressive sensing technique.

Fig. 4 .

 4 Fig. 4. Example of input signals an outputs after applying the reconstruction techniques (a) Input signal; (b) Output signal using Bayesian combined with Circulant matrix technique; (c) Output signal using Bayesian combined with random matrix technique.

Fig. 5 .

 5 Fig. 5. Example of input signals and outputs after applying the reconstruction techniques (a) Input signal; (b) Output signal using Bayesian technique combined with Circulant matrix; (c) Output signal using basis pursuit technique combined with Circulant matrix.

Fig. 6

 6 Fig.6shows the mean square error as a function of the number of samples N for the two reconstruction techniques, Bayesian and basis pursuit. As expected, for both techniques the MSE decreases with the increase of the number of samples. For N from 0 to 100, the Bayesian technique has lower MSE than the basis pursuit algorithm. For higher values of N, MSE of both techniques are slightly similar.

Fig. 6 .

 6 Fig. 6. Mean Square Error as a function of number of samples N.

Fig. 7

 7 Fig.7shows an example of results of the mean square error as a function of the sparsity level for the two techniques, Bayesian and basis pursuit. As can be seen, the MSE values corresponding to the Bayesian technique and those corresponding to the basis pursuit technique are similar and they increase with the increase of the sparsity. This figure also shows that the more the number of non-zero elements of the signal increases, the more the reconstruction becomes inefficient. One can conclude that the two techniques minimize the MSE with the same way with the increase of sparsity level.

Fig. 7 .

 7 Fig. 7. Mean square reconstruction error as a function of sparsity level k.

TABLE I .

 I TECHNIQUES COMPARISON BASED ON METRICS

		Re (%)	Cc (%)	tr (ms)	tp (ms)
	Our technique	0.77	100.00	0.90	0.96
	Basis Pursuit with Circulant matrix [6]	61.72	82.87	7.20	7.26