Dielectric and Elastic Characterization of Nonlinear Heterogeneous Materials

This review paper deals with the dielectric and elastic characterization of composite materials constituted by dispersions of nonlinear inclusions embedded in a linear matrix. The dielectric theory deals with pseudo-oriented particles shaped as ellipsoids of revolution: it means that we are dealing with mixtures of inclusions of arbitrary aspect ratio and arbitrary non-random orientational distributions. The analysis ranges from parallel spheroidal inclusions to completely random oriented inclusions. Each ellipsoidal inclusion is made of an isotropic dielectric material described by means of the so-called Kerr nonlinear relation. On the other hand, the nonlinear elastic characterization takes into consideration a dispersion of nonlinear (spherical or cylindrical) inhomogeneities. Both phases are considered isotropic (actually it means polycrystalline or amorphous solids). Under the simplifying hypotheses of small deformation for the material body and of small volume fraction of the embedded phase, we describe a theory for obtaining the linear and nonlinear elastic properties (bulk and shear moduli and Landau coefficients) of the overall material.


Introduction
The central problem of considerable technological importance is to evaluate the effective physical properties (dielectric or elastic) governing the behavior of a composite material on the macroscopic scale, taking into account the actual microscale material features [1,2]. At present, it is well known that it does not exist a universal mixing formula giving the effective properties of the heterogeneous materials (permittivity or elastic moduli) as some sort of average of the properties of the constituents. In fact, the details of the morphology or micro-geometry play a central role in determining the overall properties, particularly when the crystalline grains have highly anisotropic or nonlinear behavior or when there is a large difference in the properties of the constituent materials. The primary aim in the study of materials is to understand and classify the relationship between the internal micro-structure and the physical properties. Such a relationship may be used for designing and improving materials or, conversely, for interpreting experimental data in terms of micro-structural features. A great number of theoretical investigations have been developed in order to describe the behavior of composite materials when a specific microstructure is considered. On the other hand, a different class of theories does not assume a given microstructure, searching for general results of broad applicability. The most important properties are the classical Hashin-Shtrikman variational bounds [3,4], which provide an upper and lower bound for composite materials properties, and the expansions of Brown [5] and Torquato [6,7] which take into account the spatial correlation function of the constituents. Moreover, in [8], a functional unifying approach has been applied to better understand the intrinsic mathematical properties of a general mixing formula.
Dispersions or suspensions of inhomogeneities in a matrix are examples of widely studied heterogeneous materials: these media have been extensively analyzed both from the electrical and the elastic point of view. One of the first attempts to characterize dielectric dispersions of spheres was developed by Maxwell [9,10], who found out a famous formula valid for very diluted suspensions. The first papers dealing with mixtures of ellipsoids were written by Fricke [11,12] dealing with the electrical characterization of inhomogeneous biological tissues containing spheroidal particles: he found some explicit relationships that were simply an extension of the Maxwell formula to the case with ellipsoidal inclusions. In current literature, Maxwell relation for spheres and Fricke expressions for ellipsoids are the so-called Maxwell-Garnett Effective Medium Theories (MG-EMT) [13,14]: both theories hold on under the hypothesis of very low concentration of the dispersed component. A better model has been provided by the differential scheme [15,16]. In this case the results maintain the validity also for less diluted suspensions [17].
Some other types of microstructures have been taken into consideration. For example, the problem of the mixture characterization has been exactly solved in the case of linear and nonlinear random mixtures, that is, materials for which the various components are isotropic, linear and mixed together as an ensemble of particles having random shapes and positions (in this case there is not a material having the role of a matrix and all the media have the same importance in defining the overall properties) [18][19][20][21][22].
Recent progresses in this field concern dielectrically linear and nonlinear spheroidal inhomogeneities with geometric factors probabilistically distributed [23]. The size-dependent Bruggeman theory, which considers the effective particle dimension for non dilute dispersions, has been introduced as well [24]. A wide survey of mixture theory applications to metamaterials can be found in [25]. Finally, the dielectric (focusing or defocusing) Kerr nonlinearity [26] has been utilized to explore the importance of the particle shape [27].
On the other hand, dealing with the elastic characterization of dispersions, a similar line of research has been developed [28,29]. A famous result exists for a material composed by a very dilute concentration of linear spherical inhomogeneities dispersed in a linear solid matrix [30]. To adapt this theory to the case of any finite volume fraction, the differential method is also applied to the elastic theories for spherical [31], cylindrical [32] and ellipsoidal particles [33]. Recent works focus on microstructures that can be characterized as continuous matrices containing inhomogeneities of diverse shapes, properties and orientations [34,35]. The evaluation of the effective elastic properties of a body containing a given distribution of cracks belongs to the field of homogenization techniques as well [36]. Recent investigations consider the effects of the orientational statistical distribution of cracks in a given material [37,38].
The aim of the present review paper is to describe the dielectric and elastic characterization of composite materials constituted by dispersions of linear or nonlinear inclusions embedded in a linear matrix. The complete dielectric theory deals with pseudo-oriented particles shaped as ellipsoids of revolution: it means that we are dealing with mixtures of inclusions of arbitrary aspect ratio and arbitrary non-random orientational distributions. Moreover, the nonlinear elastic characterization takes into consideration a dispersion of nonlinear (spherical or cylindrical) inhomogeneities. Under the simplifying hypotheses of small deformation for the material body and of small volume fraction of the embedded particles, we describe a theory for obtaining the linear and nonlinear elastic properties of the overall material.
The paper is structured as follows: in Section 2 we present a brief outline of the most important results describing the dielectric homogenization techniques for linear and nonlinear dispersions. In particular, in Section 2.1. we describe the methodologies for linear and nonlinear two-phases materials (i.e., linear or nonlinear homogeneous inclusions embedded in a linear homogeneous matrix). Moreover, in Section 2.2. we describe the methods applied to investigate the three-phases materials (i.e., coated or core-shell inclusions embedded in a given matrix). We have unified all the specific results in a single framework, which is suitable for both two-and three-dimensional systems. These results are well-known in scientific literature and they are introduced here for sake of completeness. Therefore, all the relevant references have been accurately quoted in order to facilitate the interested reader.
In Section 3 we describe a microstructure constituted by pseudo-oriented ellipsoids, which is important since mimics several real materials and exhibits an overall behavior depending on the combination of two different aspects: the degree of ordering of the system (i.e., the degree of orientational distribution of the ellipsoidal particles inside the medium) and the aspect ratio (controlling the shape of the inclusion ranging from oblate to prolate spheroids) of the ellipsoids embedded in the matrix. In the following Sections 4-7 the corresponding nonlinear electric homogenization is introduced and developed. The final equations obtained at the end of the procedure are original achievements of the present work. However, several intermediate results can be found in earlier literature and will be discussed in detail in order to provide a complete review. More precisely, in Section 4 we discuss the results concerning the electrical behavior of a single nonlinear ellipsoid embedded in a linear matrix: the general theory has been presented for an arbitrary nonlinearity and the application has been performed for the so-called Kerr nonlinear constitutive equation. In Section 5 the electric field induced inside the ellipsoidal particle has been averaged over all the possible orientations of the particle itself. This calculation is original and it has been developed for a Kerr nonlinearity. Then, in Section 6 we consider a dilute suspension of randomly oriented ellipsoidal inclusions and we obtain the dielectric linear and nonlinear characterization in terms of the degree of order of the system and the aspect ratio of the particles. Finally, in Section 7 we show an example of application of the previous theory and we discuss the behavior of the effective permittivities and of the nonlinear susceptibilities.
In Section 8 we introduce the homogenization for nonlinear elastic composite materials. In particular we describe the microstructures analyzed in the following sections and we introduce the applied methodologies. This second part of the paper represents a detailed review of recent results in the field of the nonlinear homogenization (the relevant references will be properly quoted). In Section 9 we discuss the nonlinear constitutive equations adopted to model the embedded particles (for the nonlinear elastic homogenization scheme we consider for simplicity spherical or cylindrical inhomogeneities). In Section 10 we introduce the nonlinear generalization of the Eshelby theory, which allows us to determine the elastic fields induced inside an inhomogeneity. Finally, in Sections 11 and 12 we develop the procedures for dealing, respectively, with dispersions of spheres and cylinders.

Introductory Remarks on Linear and Nonlinear Dielectric Homogenization
In order to introduce the standard homogenization techniques for linear and nonlinear heterogeneous structures, we consider a two-phases material (spheres or cylinders embedded in a given matrix) and a three-phases material (coated spheres or cylinders embedded in a given matrix). In both cases we analyze the electric behavior of the single particle and of a given assembly of inclusion.

Two-phases materials
We describe the linear electric behavior of a given spherical (d = 3) or cylindrical (d = 2) particle with permittivity 2 and radius R embedded in a matrix with permittivity 1 . For the spherical object we adopt a system of spherical coordinates (ρ, θ, φ) (x = ρ sin θ cos φ, y = ρ sin θ sin φ, z = ρ cos θ) and for the cylindrical particle we use a system of cylindrical coordinates (ρ, θ, y) (x = ρ sin θ, y = y, z = ρ cos θ). We suppose to apply a remote electric field E ∞ along the z-axis in both cases. The electric potentials inside and outside the inclusion, φ i and φ e , are given by [14,15] Therefore, the electric field induced inside a particle is This internal field is called Lorentz field and it is always uniform: it is possible to prove that Figure 1 one can observe the lines of the electric field for an inclusion with 2 > 1 and for the case 2 < 1 . In both cases the uniformity of the internal field and the dipolar character of the external one are evident. Figure 1. Field lines of the electric field for an inclusion with 2 > 1 and for another one with 2 < 1 .
We consider now a dispersion of particles in a given region. We assume that the volume fraction c of inclusions is very low (c 1): it means that Equation (3) is valid for each particle and that the average value of the electric field over the entire structure can be obtained as E = cE i + (1 − c) E ∞ . In similar way it is possible to calculate the average value of the displacement vector as D = c ( 2 − 1 ) E i + 1 E . We define the concept of effective permittivity ef f through the relation D = ef f E . The above expressions allow us to obtain the following important result (Maxwell formula) [9] This is the main result concerning the linear homogenization of an assembly of spherical (d = 3) or cylindrical (d = 2) particles.
We consider now the nonlinear behavior of a single inclusion. This methodology has been developed in [39] and it has been utilized both for a single particle and an assembly of nonlinear spheres. The benefit of this approach is that of considering an arbitrary nonlinearity describing the electric behavior of the inhomogeneities embedded in the matrix. The nonlinear constitutive relation for the particle can be written as D i =˜ 2 (E i )E i where˜ 2 (E i ) is the field-dependent permittivity. The important Equation (3) is still valid but now it is an implicit equation giving the actual internal fieldẼ ĩ Therefore, the balance equations for a dilute dispersion of nonlinear particles are the following From the first two expressions in Equation (6) we can eliminate the remote field E ∞ , by obtaining an explicit relationship between the internal fieldẼ i and the average value E Now, it is necessary to solve Equation (7) with respect toẼ i for different values of the average electric field E (nonlinear microdosimetry). This point can be accomplished with different numerical methods depending on the complexity of the nonlinearity˜ 2 (Ẽ i ). If we substitute Equation (7) in the last expression given in Equation (6) we obtain the nonlinear effective constitutive equation Therefore, the effective nonlinear dielectric constant (depending on E ) is given by The explicit form ofẼ i in terms of E (which can be obtained numerically or analytically from Equation (7)) must be substituted in Equation (9) in order to obtain the nonlinear effective behavior. This method can be applied to the case of a Kerr nonlinearity for the particles [27] where the dots indicate the phasors of the corresponding quantities (we assume a sinusoidal permanent regime). In this case both 2 and α can be considered complex numbers. We may substitute Equation (10) in Equation (7), obtaininġ where we have definedĖ i = X (the main unknown). This equation assumes the simple form The solution can be approximately obtained through the following expansion truncated after the third odd term By using Equation (13) in Equation (12) we obtain the values for the coefficients Finally, adopting the value ofĖ i = X in Equation (9) we obtain the nonlinear homogenization where [27] Here we have used the definition A = (˙ 2 − 1 ) (1 − c) + d 1 . The first expression in Equation (16) represents the Maxwell formula describing the linear behavior of the composite system. The second relation in Equation (16) represents the Kerr nonlinear coefficient of the dispersion: we observe that it is directly proportional to the volume fraction c and to the Kerr parameterα of the inclusions. When c → 1 the value of˙ ef f approaches˙ 2 andα ef f approachesα, as expected. Finally, the coefficientβ ef f represents the fourth nonlinear term of the heterogeneous structure. We note that it is proportional to c (1 − c): therefore, it is zero for c = 0 (no inclusions in the system) and for c = 1 (the inclusions fill the whole space). Anyway, we remember that all the results are valid for dilute dispersions (c 1). Similar results and some generalizations can be found in [40][41][42].

Three-phases materials
We consider a spherical (or cylindrical) structure constituted by a matrix with permittivity 1 where a coated particle is embedded. The particle is composed by a core with radius a and permittivity 3 and a shell contained between the radii a and b with permittivity 2 . On such a system an electric field E ∞ = E ∞ e z is remotely applied. The electric potentials in the three regions can be calculated as [39,43] As before, these expression are valid both in cylindrical symmetry (d = 2) and in spherical symmetry (d = 3). The coefficients A, B, C, D can be obtained by imposing the continuity of the electric potential and of the normal component of the electric displacement. These conditions allows us to obtain the system of equations The exact solutions have been eventually obtained [39] To begin with the homogenizing schemes, we search for a homogeneous dielectric material (with per-mittivity˜ ) in the region 0 < ρ < b, which is equivalent to the structure above described. It means that the substitution of the regions with 2 and 3 with a single homogeneous particle with permittivity˜ does not alter the values of the electric potential and of the electric field in the external region ρ > b. In the original heterogeneous structure the external electric potential is given by φ 1 reported in Equation (17).
In the equivalent case with a single uniform inhomogeneity this electric potential is described by the same expression where A = A | 2 =˜ , 3 =˜ since the theory is still correct for 2 = 3 =˜ . By equating the values A ( 1 , 2 , 3 ) = A ( 1 ,˜ ,˜ ), we obtain an equation for the unknown˜ . The solution is We remark that this value of the effective permittivity does not depend on the permittivity of the matrix. In fact, it is a characteristic quantity of the heterogeneous structure (constituted by the core and the shell) embedded in the homogeneous material with permittivity 1 . Moreover, it depends on the ratio between the radius of the core and the radius of the shell. Therefore, we can define the volume fraction of the core of permittivity 3 in the shell of permittivity 2 asc = a b d . Therefore, Equation (20) assumes the very simple form˜ This is an exact result. It is interesting to note that Equation (21) is formally identical to the Maxwell formula, which describes the effective permittivity of a dilute dispersion of inclusions of permittivity 3 in a hosting material of permittivity 2 . Nevertheless, the Maxwell formula is an approximate solution which is valid only for low values of the volume fractionc of the spheres embedded in the matrix. We have obtained a uniform electric field in the core (with intensity −D) and a more complicated spatial behavior in the shell. For the following development it is useful to calculate the average value of the electric field in the shell [39]. More precisely, we search for the average value of E 2 = − ∇φ 2 over the region a < ρ < b where Ω is the region corresponding to the shell (or coating). When d = 2, we solved the double integral Moreover, for d = 3 we obtained It means that, in each case, the relation E 2 = −Be z is always fulfilled. Now, we consider a dispersion of coated particles in the linear electric regime. We utilize the property which homogenizes each particle, in order to apply the Maxwell formula in a second step. We define 1 where V is the total volume of the material, V d is the total volume of the effective inclusions with permittivity˜ and V f is the volume of the external region with permittivity 1 . At the end of this procedure we obtain where˜ is given in Equation (21). The optical particles of core-shell particle composites has been also analyzed for the confocal ellipsoids configuration: interested readers can consult [44,45] for more information. This procedure can be generalized to the case with a nonlinear core in the particles [39,46,47]. This model has been applied to investigate the optical bistability of suspensions of nonlinear coated nanoparticles [48,49] and the orientation of core-shell inclusions in an electric field [50]. We suppose that each core can be described by a Kerr nonlinearity˜ 3 = 3 + χ|E 3 | 2 , where E 3 is the electric field for ρ < a (for simplicity we omit the dots indicating the phasors). From the above discussion we obtain E 3 = −D where D is given in Equation (19). In nonlinear regime, the coefficient D is not constant since it depends on˜ 3 and, therefore, on the core field E 3 . Explicitly By recalling the definitionc = a b d we can write the denominator of the previous relation as where we have defined the constants α and β as follows Now, Equation (26) can be written in the simple form At this point it is useful to obtain the average value of the electric field inside the composite material. We develop our calculations along the z direction. The value of E z can be obtained as follows where we have used the property E 2 = −B above discussed. Here, B is not constant since it depends on˜ 3 . We have (note that D and B have the same denominator) From Equation (29) we can evaluate˜ 3 in terms of α and β Therefore, the coefficient −B assumes the simpler form Now, we can obtain E ∞ in terms of E 3 from Equation (29) Substituting the previous expression in Equation (33) we obtain the average value of the electric field in the coating shell The average value of the electric field over the entire structure is then obtained from Equation (30) where we have defined Now, Equation (36) is in the form E z =ãE 3 +b|E 3 | 2 E 3 , which is similar to Equation (12). It can be simply solved with respect to E 3 by considering the first two terms (of order one and three) We are searching for the effective dielectric constant and, therefore, for the relationship between the average electric field and the average electric displacement. Then, it is important to evaluate the average electric displacement (in direction z) where D 1 , D 2 and D 3 are the displacement vectors in the three phases. We use Equation (38), giving the value of E 3 , in order to approximate the quantity |E 3 | 2 . We obtain Therefore, the average displacement field is given by where we have defined the effective quantities The results given in Equations (41) and (42) represent the complete nonlinear characterization of a dilute dispersion of coated particles with nonlinear core. It is not difficult to verify that the result for ef f is coincident with Equation (25), as expected [39]. Moreover, the second relation for χ ef f represents a complementary result, describing the nonlinear behavior of the heterogeneous structure.
In this introductory section we have summed up the most important methodologies and results regarding simple geometries of heterogeneous structures. In the following we will use similar techniques to analyze more complicated microstructures.

Nonlinear Electric Homogenization for Pseudo-Oriented Ellipsoids
In recent material science development, considerable attention has been devoted to electromagnetically nonlinear composite structures due to their applications, for instance, to integrated optical devices (such as optical switching and signal processing devices) [51][52][53]. More specifically, intrinsic optical bistability has been extensively studied theoretically as well as experimentally with the help of mixture theory [54,55]. In all of these cases, a linear medium containing spherical or spheroidal inclusions has been considered. Important results concerning a dispersion of dielectrically nonlinear and graded parallel cylinders have been achieved by Wei et al. [56]. Our aim is to extend previous works and the techniques discussed in the previous section and to explore the importance of the orientational distribution and of the inclusion shape in this context. To do this, we consider a dispersion of dielectrically nonlinear spheroidal particles (ellipsoids of revolution), pseudo randomly oriented in a (dielectrically) linear matrix and we then develop a mathematical procedure to perform the needed averages of the electric quantities over all possible orientations of the inclusions. This analysis leads to the nonlinear anisotropic constitutive equation connecting the macroscopic electric displacement to the macroscopic electric field. A particular attention is devoted to the analysis of the effects of the orientational distribution of the particles inside the composite material. The limiting cases of the present theory are represented by all the particles aligned with a given direction (perfect order) and all the particles randomly oriented (complete disorder). We take into account all the intermediate configurations between order and disorder with the aim to characterize a material with particles partially aligned. In Figure 2 one can find some orientational distributions between the upon described limiting cases. Structure of a dispersion of pseudo-oriented ellipsoids. One can find some orientational distributions ranging from order to disorder. The two-phase material is described by the electric response of each phase, by the state of order and by the volume fraction of the inclusions.
To define the geometry, we consider a given orthonormal reference frame and we take as preferential direction of alignment the z-axis. Each particle embedded in the matrix is not completely random oriented. The orientation is described by the following statistical rule: the principal axis of each particle forms with the z-axis an angle ϑ, which follows a given probability density f Θ (ϑ) symmetrically distributed in [0, π]. The symmetry of the density can be written as f Θ (ϑ) = f Θ (π − ϑ). We assume that the orientation of each particle is statistically independent from the orientation of other particles. If f Θ (ϑ) = (1/2)(δ (ϑ) + δ (ϑ − π)) (where δ is the Dirac delta function) we have all the particles with ϑ = 0 (or ϑ = π, which corresponds to the same orientation) and, therefore, they are all oriented along the z-axis. If f Θ (ϑ) = sin (ϑ) /2 all the particles are uniformly random oriented in the space over all the possible orientations. Any other symmetric statistical distribution f Θ (ϑ) defines a transversely isotropic (uniaxial) material with principal axis aligned with the z-axis. For example, if f Θ (ϑ) = δ (ϑ − π/2), all particles have the principal axis orthogonal to the z-axis. In the following sections we develop a complete analysis of the combined effects of the shape (aspect ratio or eccentricity) of the particles and their state of order/disorder. This analysis allows us to evaluate the overall electric properties of the heterogeneous material. In particular, from the point of view of the shape of the particles, the so-called depolarization factor L is the parameter that intervenes to characterize the medium. We verified that the state of order acts on the overall linear and nonlinear dielectric properties by means of two parameters that are defined as follows: C 2 = cos 2 ϑ and C 4 = cos 4 ϑ . They correspond to the average values of cos 2 ϑ and cos 4 ϑ, computed by means of the density probability f Θ (ϑ). The results may be applied to describe the physical behavior of heterogeneous materials starting from the knowledge of the physical properties of each medium composing the mixture as well as of the structural properties of the mixture itself, i.e., shape of the inclusions and state of order of the orientations (L, C 2 and C 4 ). It is worth pointing out, as it frequently occurs in this field, that the presented results have been derived under electrostatic assumption, but they hold valid also in the low frequency regime, as long as the wavelength is much larger than the largest dimension of the inclusions. The analysis performed in the following has immediate application to the field of the liquid crystals. Actually, our microstructure describes a material positionally disordered, but with partial orientational order, which corresponds to a nematic phase in liquid crystals [57,58]. The level of ordering is reflected in the macroscopic properties. Some previous works have been devoted to an analysis similar to that developed in this work but only from a dielectrically linear point of view [59][60][61][62]. So, the following development can be considered as a nonlinear extension of such previous ones.

Field Perturbation Due to One Single Nonlinear Ellipsoidal Inclusion in a Uniform Field
Here we present a general solution to the problem of a nonlinear ellipsoidal particle embedded in a linear material. The theory is based on the following result derived for the linear case, which describes the behavior of one electrically linear ellipsoidal particle of permittivity ε 2 in a linear homogeneous medium of permittivity ε 1 . Let the axes of the ellipsoid be l x , l y and l z (aligned with axes x, y, z of ellipsoid reference frame) and let a uniform electric field E 0 = (E 0x , E 0y , E 0z ) be applied to the structure. Then, according to Stratton [43], the electric field E s = (E sx , E sy , E sz ) inside the ellipsoid is uniform and it can be expressed as follows Here, and throughout the paper, the index i takes the x, y and z values. The expressions for the depolarization factors L i in the case of generally shaped ellipsoid can be found in the literature [17]. They can be expressed in terms of elliptic integrals. The condition L x + L y + L z = 1 is always fulfilled. Let's now generalize Equation (43) to the case where a dielectrically nonlinear ellipsoid is embedded in the linear matrix. A nonlinear isotropic and homogenous ellipsoid can be described from the electrical point of view by the constitutive equation D = ε (E) E [63]. Here, D is the electric displacement inside the particle, E is the electric field and the function ε depends only on the modulus E of E. This latter property accounts for the fact that the medium inside the ellipsoid is isotropic and homogenous. The main result follows. The electric field inside the inclusion is uniform even in the nonlinear case and it may be calculated by means of the following system of equations [27] where, as before, E 0 is a uniform electric field applied to the structure and E s , the unknown in the nonlinear system (44), is a uniform field as well. This property holds true due to the following reason: if a solution of (44) exists, due to self-consistency, all the boundary conditions are fulfilled and the problem is completely analogous to its linear counterpart, treated by Stratton [43], provided that ε 2 = ε (E s ).
In order to simplify the following analysis we will adopt ellipsoids of revolution. Thus we consider l x = l y and we define the aspect ratio as e = l z /l x = l z /l y . The depolarization factors for ellipsoids of revolution may be computed in closed form as follows and the results depend on the shape of the ellipsoid [17]. It is prolate (of ovary or elongated form) if e >1 and oblate (of planetary or flattened form) if e <1 where p = √ e 2 − 1 and q = √ 1 − e 2 . The relation 2L x +L z = 1 holds on and therefore we will consider L = L z as main geometric parameter of the system. An interesting aspect related to the problem faced in this section shows up when one considers the nonlinear Equation (44) and tries to solve it iteratively [27]. This means that, in order to solve for E s , one starts with a given initial value E 0 s , and one uses the successive approximations described by the iteration rule The following sufficient convergence criterion has been verified [27]: the iteration rule given by Equation (47) is convergent to the exact internal electric field if the nonlinear material of the ellipsoid fulfils the condition E ε ∂ε ∂E < 1. In a general context, one can describe isotropic nonlinear dielectric materials by means of the so-called Kerr nonlinearity relation, often adopted in metamaterials study which assumes that ε 2 and α are constant. The Kerr nonlinearity is called focusing or defocusing according to the fact that α > 0 or α < 0, respectively [63]. It is straightforward to verify that the convergence condition E ε ∂ε ∂E < 1 is always verified for defocusing Kerr nonlinearity and is verified only if E 2 s < ε 2 /α (here E s is the modulus of the actual electric field inside the inclusion) in the case of focusing nonlinearity [27].

Average Electric Field Inside a Single Pseudo-Random Oriented Inclusion
Now, our aim is to find an explicit version of Equation (44), which is valid when the nonlinear permittivity is given by Equation (48). To begin the analysis, we substitute Equation (48), holding for a single ellipsoid, in Equation (44) This is an algebraic system of degree nine with three unknowns, namely E sx , E sy , and E sz . It might be hard, if not impossible, to be solved analytically, but we are interested, for our purposes, in just the first terms of a series expansion for the solution. To obtain it, we may adopt the ansatz and solve for k i and h i . Alternatively, we may use the iterative scheme given in Equation (47), in literal form, adopting only the first iterations. For sake of brevity, we omit here the simple but long calculation, which leads to the solution We observe that the first term represents the classical Lorentz field appearing in a dielectrically linear ellipsoidal inclusion. The second term is the first nonlinear contribution, which is directly proportional to the inclusion hyper-susceptibility α. To simplify the expressions, from now on, we will use the notation: To derive the mixture behavior, we need to calculate the electric field in a single nonlinear ellipsoidal inclusion arbitrarily oriented in space and embedded in a homogeneous medium with permittivity ε 1 . In order to do this, we shall express Equation (50) in the global framework of reference of the mixture. We define three unit vectors, indicating the principal directions of each ellipsoid in space: they are referred to asn x ,n y andn z , and they correspond to the axes of the ellipsoid. By using Equation (50), we may compute the electric field induced by a given external arbitrary uniform electric field inside the inclusion (from now on we will omit the additional higher order terms) We shall now average it over all the possible orientations of the particle. The expression for the internal electric field in Equation (51) can be rewritten component by component as follows where n jk is the k-th component of the unit vectorn j , (j = x, y, z) and we have considered the implicit sums of i, j, l, q and pover 1, 2 and 3. For the following derivation, we are interested in the average value of the electric field E s over all the possible orientations of the ellipsoid itself and then we have to compute the following We may use Euler angles representation (ψ, ϕ and ϑ) to write down the explicit expressions for the components of the unit vectorsn x ,n y andn z     n x = (cos ψ cos ϕ − sin ψ sin ϕ cos ϑ, − cos ψ sin ϕ − sin ψ cos ϕ cos ϑ, sin ψ sin ϑ) n y = (sin ψ cos ϕ + cos ψ sin ϕ cos ϑ, − sin ψ sin ϕ + cos ψ cos ϕ cos ϑ, − cos ψ sin ϑ) n z = (sin ϕ sin ϑ, sin ϑ cos ϕ, cos ϑ) (54) In order to obtain the average value of the electric field E s , we need to calculate the average value of the quantities defined in Equation (53). This is done by the following integral over the Euler angles In order to represent the above described orientation of the particles, the angles ϕ and ψ are uniformly distributed over the entire range [0 2π] and the angle ϑ follows the given probability density f Θ (ϑ) over the range [0 π]. So, by performing the integration described in Equation (55) and by using Equation (54), we may obtain that the average value of E sk depends on the two following parameters, defined by means of the density probability f Θ (ϑ) These two parameters completely characterize the effects of the pseudo-orientation of the particles inside the medium. Some particular values follow. When we are in the case of perfect order we have f Θ (ϑ) = (1/2)(δ (ϑ) + δ (ϑ − π)) and the corresponding values are C 2 = 1 and C 4 = 1. If we are in the case of complete disorder we have f Θ (ϑ) = sin (ϑ) /2 and we obtain C 2 = 1/3 and C 4 = 1/5. Finally, when the particles all have the principal axis orthogonal to the z-axis we have f Θ (ϑ) = δ (ϑ − π/2) and the values of the parameters are C 2 = 0 and C 4 = 0.
Since we are dealing with ellipsoids of revolution, in performing the integration of Equation (55) we use the simplified notation a 1 = a 2 and L 1 = L 2 = (1 − L)/2, L 3 = L. The factor L assumes some characteristic values in correspondence to special shapes of the particles: for spheres L = 1/3, for cylinders L = 0 and for lamellae or penny shaped inclusions L = 1. Summing up, we verified, after a very long by straightforward integration, that the following simple relation gives the final result of the averaging process Here, the sum on the index l is implied and the parameters γ k and µ kl can be organized as follows The explicit results for the parameters γ k are Moreover, the explicit results for the parameters µ kl are This is a first analytical result, which will play a crucial role in the following development of the theory. It is interesting to observe that if a 1 = a 3 and L = 1/3 (we are dealing with spherical inclusions), then the terms containing C 2 and C 4 completely disappear in Equations (59)- (64): it is correct since the orientation is not important for an isotropic spherical object.

Averaging Process in a Dilute Mixture
From now on we analyze the dispersion of pseudo-oriented nonlinear ellipsoids. The permittivity of the inclusions is described by the isotropic nonlinear relation ε (E) = ε 2 + αE 2 [see Equation (48)] and the linear matrix has permittivity ε 1 ; the overall electrical behavior of the dispersion is expected to be anisotropic because of the pseudo-random orientation of the particles. This is true because the z-axis has a special character induced by the partial alignment of the particles. Therefore, the equivalent electric constitutive equation can be expanded in series with respect to the averaged electric field components: where the coefficients ε eq (the superscript eq pointing out the equivalent character of the term) and χ eq are tensors depending on various parameters of the mixture such as the aspect ratio e of the ellipsoids, the volume fraction c of the included phase, the density probability f Θ (ϑ) describing the orientational distribution, the permittivities ε 1 , ε 2 and the Kerr susceptibility α of the inclusions. The homogenization procedure should provide the structure of the entries of the tensors ε eq and χ eq in terms of the mentioned parameters. In the technical literature, the coefficients α and χ eq (the first nonlinear terms of the expanded constitutive equations for inclusions and mixture, respectively) are often called hyper-susceptibilities [27].
The main achievement of this work is the derivation of a closed form expression for the hypersusceptibility ratio χ eq /α. These quantities are of interest as much as they represent the amplification of the composite material nonlinear behavior with respect to that of the inclusions. In particular we are interested in the dependence of these parameters on the state of order/disorder of the system, which is well described by the density probability f Θ (ϑ). In other words, it means that we will write our results in terms of the order parameters C 2 and C 4 . Moreover, we may describe the dependence of the hypersusceptibility ratio χ eq /α on the aspect ratio of the embedded particles, i.e., on the parameter e or L. The final expressions are derived under the assumption that the constitutive equation of the composite medium is of the form D k = ε eq kj E j +χ eq kjil E j E i E l , which neglects higher order terms. All the computations are carried out under the same hypothesis underlying the linear Maxwell-Garnett theory, that is, low concentration c of the dispersed phase. So, if we consider a mixture with a volume fraction c << 1 of pseudo-randomly oriented, dielectrically nonlinear ellipsoids embedded in a homogeneous matrix with permittivity ε 1 , we can evaluate the average of the electric field over the space occupied by the mixture. It can be done via the following relationship This means that we do not take into account the interactions among the inclusions because of the very low concentration: each ellipsoid behaves as an isolated one. Once more, to derive Equation (65), we assume an approximately uniform average electric field E 0 in the space outside the inclusions. To evaluate the equivalent constitutive equation, we compute the average value of the displacement vector inside the random material. The region V is defined as the space occupied by the mixture, Ve as the region occupied by the inclusions, and Vo as the remaining space (so that V = Ve∪Vo). The average value of D ( r) = ε E ( r) is evaluated as follows ( D and E represent the local fields, D and E their macroscopic counterparts) Here |V |is the measure of the region V . It can be noted that the average value given by the expression ε E s − ε 1 E s is not available from the previous computations. We consider a single ellipsoidal nonlinear inclusion and we search for the average value of the quantity ε E s − ε 1 E s over all the possible orientations of the particle. From Equations (44) and (50) we obtain therefore, in vector notation By taking the k-th component in the global reference framework, we may write and averaging, after some straightforward computation At this point, the explicit average value can be found by taking into consideration the expressions of the unit vectorsn x ,n y andn z , given in Equation (54), and performing the integration in a similar way to that shown in Equation (55). The effects of the pseudo-orientation of the particles inside the medium are described, as before, by the order parameters C 2 and C 4 . Of course, performing the averaging in Equation (70) we have used again the simplified notation a 1 = a 2 and L 1 = So, the final result for the quantity D sk − ε 1 E sk is given by the following simple relation Here the sum on the index l is implied and the parameters γ k have been defined in Equations (58), (59) and (60). Moreover, the parameters λ kl can be arranged in following matrix notation The explicit values of the relative entries have been calculated as follows Once again, we observe that if a 1 = a 3 (we are dealing with spherical inclusions), then the terms containing C 2 and C 4 completely disappear in Equations (73)-(75): the orientation is not relevant for an isotropic spherical object. At this point we have all the balance equations needed to describe the overall electrical behavior of the pseudo random dispersion. These relationships have been summarized in Equation (76). The first relation corresponds to Equation (65) and furnishes the average electric field over the mixture volume in terms of the applied field and the average internal field (equation deduced under the hypothesis of low concentration). The second relation is taken from Equation (57) and gives the explicit value of the average electric field inside an inclusion (it accounts just for the first nonlinear terms). The third relation [see Equation (66)] furnishes the average value of the displacement vector over the entire mixture volume (this formula is exact). Finally, the fourth equation is taken from Equation (71) and gives the average value of the quantities D sk − ε 1 E sk in terms of the applied electric field (as before, it accounts just for the first nonlinear terms). The complete set of the balance equations follows We may observe that the nonlinear terms, appearing in the second and in the fourth equations, are simply proportional to the hypersusceptibility α of the inclusions. By substituting the second equation in the first one and the fourth relation in the third one we obtain a simpler system Here, we have three vector fields involved: the average electric field over the entire mixture, the average electric displacement and the external applied electric field. In order to find out the effective constitutive equation for the whole composite material, we should obtain a relationship among the components D k and the components E h . So we have to eliminate the external field E 0k in Equation (77). Therefore, now we need to solve the first relation in Equation (77) with respect to E 0k : for our purposes it is sufficient to obtain a series solution with two terms and thus we let substitute it in the first relation in Equation (77) and we solve for the unknown coefficients g k and m kl . The result is The final achievement is obtained by substituting Equation (78) in the second equation of system (77) and neglecting the powers of E k greater than three This is the first form of the constitutive equation of the overall dispersion. This result can be further simplified by defining the following quantities, which better describes the transversely isotropic (uniaxial) character of the composite material: The symbol // indicates the components aligned with the principal axis (the z-axis) and the symbol ⊥ indicates the components orthogonal to the principal axis. With such conventions Equation (79) may be rearranged as follows The linear permittivities ε ⊥ and ε // and the nonlinear hyper-susceptibilities χ ⊥,⊥ , χ ⊥,// , χ //,⊥ and χ //,// can be derived by comparison with Equation (79) and the relative explicit expressions are given below All these quantities are the main parameters describing the nonlinear electrical behavior of the overall dispersion. We may observe that the nonlinear susceptibilities χ ⊥,⊥ , χ ⊥,// , χ //,⊥ and χ //,// are proportional to the susceptibility α of the inclusions and depend on the factors γ k , µ kl and λ kl defined in Equations (59)- (64) and Equations (73)- (75). It is worth pointing out that, as it is expected, the results are explicitly written in terms of the depolarizing factor L z = L of the inclusions, which directly depends on the aspect ratio e [see Equation (45)], and in terms of the order parameters C 2 and C 4 that define the state of orientational order/disorder, which depends on the probability density f Θ (ϑ) (see Equation (56)). Finally, the particular cases of spherical inclusions (a 1 = a 3 and L = 1/3) and of ellipsoidal inclusions with isotropic distribution (C 2 = 1/3 and C 4 = 1/5) provide results in perfect agreement with previous investigations [27].

Example of Application
In order to show some results of the previous procedure we choose a particular probability density f Θ (ϑ) that depends on one parameter a. This probability density is particularly useful because when the parameter a varies from −∞ to +∞ the orientational distribution of the inclusions assumes all the interesting possibilities. More precisely, when a → −∞ we are in the case of perfect order and all the particles are aligned with the z-axis, when a = 0 we are in the case of complete disorder (all the particles uniformly random oriented in the space) and when a → +∞ all the inclusions have the principal axis orthogonal to the z-axis. The expression of the normalized probability density over the range [0, π] follows The function is symmetrical with respect to ϑ = π/2. The above statement can be also formulated as follows: if a → −∞ one can verify that f Θ (ϑ) = (1/2)(δ (ϑ) + δ (ϑ − π)), where δ is the Dirac delta function (perfect order); if a = 0 we obtain f Θ (ϑ) = sin (ϑ) /2 (complete disorder); finally, if a → +∞ it is possible to show that f Θ (ϑ) = δ (ϑ − π/2) and all the particles have the principal axis orthogonal to the z-axis. In Figure 3 one can find the shape of this probability density in correspondence to three different values of the parameter a (a = −4, a = 0 and a = 4). One can observe that, for negative values of a, we obtain a bimodal density, for a = 0 we obtain the sine shaped function f Θ (ϑ) = sin (ϑ) /2 and for positive value of a, we have a unimodal behavior. This probability density is particularly useful also because it allows computing the order parameters C 2 and C 4 in a closed form The previous expressions furnish the following special values: if a → −∞ we have C 2 = C 4 = 1, if a = 0 we have C 2 = 1/3 and C 4 = 1/5 and, finally, if a → +∞ we obtain C 2 = C 4 = 0. In Figure 4 one can find the behavior of the coefficients C 2 and C 4 versus the parameter a. Figure 4. Behavior of the order parameters C 2 and C 4 in terms of the coefficient a as described by Equations (83) and (84), respectively.
We have written a software code that implements the complete procedure summed up in Equation (81), in order to obtain the macroscopic linear and nonlinear features of the composite material in terms of the aspect ratio e of the ellipsoids and of the parameter a controlling the state of order, as above described. A first series of results concerns the case with ε 1 = 1, ε 2 = 10 and c = 0.25. They are shown in Figures 4-9 versus a and Log 10 (e). More precisely, one can find: in Figure 4 the permittivity ε ⊥ , in Figure 5 the permittivity ε // , in Figure 6 the susceptibility amplification Log 10 (χ ⊥,⊥ /α), in Figure 7 the susceptibility amplification Log 10 χ ⊥,// /α , in Figure 8 the susceptibility amplification Log 10 χ //,⊥ /α and, finally, in Figure 9 the susceptibility amplification Log 10 χ //,// /α . A second series of results, concerning the case with ε 1 = 1, ε 2 = 0.1 and c = 0.25 can be found in Figures 10-15. As before, they have been represented in terms of a and Log 10 (e) and we have adopted the same ordering for the plots. We may draw a comparison among the results obtained when ε 2 /ε 1 = 10 and the results obtained when ε 2 /ε 1 = 1/10: in both cases, as for the longitudinal and transversal permittivities, we may observe that the effect of the order/disorder has opposite behavior for prolate and oblate particles. Moreover, the complex behavior of the susceptibility amplifications is inverted moving from the case with ε 2 /ε 1 = 10 to the case with ε 2 /ε 1 = 1/10. The plots exhibit a very complex scenario for the macroscopic properties of the nonlinear material, strongly dependent on the state of order and on the geometric features of the embedded ellipsoids of revolution (prolate or oblate).
We draw some conclusions about the homogenization procedure described. We have analyzed the nonlinear dielectric effects of the orientational order/disorder of non-spherical particles in composite or heterogeneous materials. As result of this analysis we have found the correct definition of two order parameters (C 2 and C 4 ) in such a way to predict the macroscopic electric properties as function of the state of microscopic order. In particular, we have found out explicit relationships that allow us the calculation of the linear permittivity tensor and the nonlinear susceptibility tensor in terms of the shape of the embedded particles and the order parameters. We have outlined and applied a complete procedure which takes into account any given orientational distribution of ellipsoids in the matrix. The theory can find many applications to real physical situations ranging from technological aspects of composite materials to optical characterization of nematic liquid crystals and to tissues modeling in biophysics.

Nonlinear Elastic Homogenization
The theoretical approaches utilized to analyze nonlinear elastic composite materials are typically based on rigorous variational principles which, in addition to possessing mathematical rigor, have the advantage of leading to bounds and relatively accurate estimates for the mechanical properties. Such variational principles allow one to obtain estimates of the effective energy densities of nonlinear materials in terms of the corresponding information for linear composites with the same microstructure. These methodologies can be found on an excellent review by Ponte Castañeda and Suquet [64]. From the historical point of view, the variational procedure of Hashin and Shtrikman [3,4] is the first important result concerning the linear behavior of electric and elastic heterogeneous materials. This variational procedure provides lower and upper bounds on the elastic moduli and elastic tensors for isotropic composites (reinforced by randomly positioned particles). A generalization of the Hashin-Shtrikman variational principles, suitable for nonlinear materials, was developed by Talbot and Willis [65,66]. This extension can be used to obtain improved bounds (depending on a two-point statistical information) for nonlinear composites. Variational methods for deriving improved bounds and estimates for the effective properties of nonlinear materials, utilizing the effective modulus tensor of suitably selected linear-elastic comparison materials, were introduced by Ponte Castañeda [67][68][69][70] for materials with isotropic phases and by Suquet [71] for composites with power-law phases. Moreover, a hybrid of the Talbot-Willis and Ponte Castañeda procedures, using a linear thermoelastic comparison material, was proposed by Talbot and Willis [72]. An important advantage of the variational procedures that involve linear comparison materials is that, they can not only produce the nonlinear Hashin-Shtrikman-type bounds of the Talbot-Willis procedure directly from the corresponding linear Hashin-Shtrikman bounds, but also yield higher-order nonlinear bounds, such as Beran-type bounds [73], as well as other types of estimates.
The general nonlinear elastic features are relevant in many materials science problems. For example, in biomechanics, transient elastography has shown its efficiency to map the nonlinear properties of soft tissues and it can be used as diagnostic technique [74,75]. In material science the linear theory is incapable of fully capturing all fracture phenomena and hyperelasticity plays a governing role in the dynamics of fracture [76,77]. The quantum dots growth, ordering and orientation (occurring during processing) are largely affected by elastic phenomena, even beyond the linear regime [78,79]. Finally, many problems of fracture mechanics in composite materials do contain nonlinear features like, e.g., the interaction between a crack and a fiber (or, more generally, an inclusion) [80].
The aim of the remaining part of the present work is to review the elastic properties of dispersions of nonlinear elastic inclusions embedded in a linear elastic hosting matrix. Here, we do not apply methodologies based on variational principles since we can obtain, in this particular case, the estimate of the effective elastic behavior by means of the direct calculations of the elastic fields. Therefore, we describe a procedure similar to that utilized, in the previous sections, for the dielectric homogenization. In particular, we will study the elastic fields induced in a single nonlinear particle and then we use such results to homogenize complex dispersions. It is known that the concept of nonlinearity can be introduced in the theory of elasticity in two different ways [81]. A first nonlinearity can be taken into account by means of the exact relation for the strain (not limited for small deformation) and the exact equilibrium equations for a volume element of the body (this first aspect is referred to as geometrical nonlinearity since it is related to the equations not depending on the material under consideration). Secondly, another nonlinear effect can be considered through the arbitrariness of the (generically not Hookean) stress-strain constitutive relation (this aspect is referred to as physical nonlinearity since it is related to properties of the material under consideration). Therefore, by combining the two previous contributions, it follows that there are four different types of problems in the theory of elasticity [82] • those having both physical and geometrical linearity; • those which are physically nonlinear but geometrically linear; • those linear physically but nonlinear geometrically; • those nonlinear both physically and geometrically.
The problems of the first type are the subject of the (classical) theory of elasticity (small deformation in Hookean materials). In this review, we adopt the second conceptual framework. The angles of rotation can be neglected in determining changes in dimensions in the line elements and in formulating the conditions of equilibrium of a volume elements: therefore, the balance equations are based on the standard small-strain tensor and on the Cauchy stress tensor (typically introduced in the problem of the first type). However, the elongations exceed the Hookean limit of proportionality (between stress and strain) and this requires a nonlinear stress-strain relationship. This conceptual framework is sometimes referred to as hypoelasticity: it is intended to model perfectly reversible nonlinear stress-strain behavior but restricted to infinitesimal strains. Such a description has been already adopted in the past in order to model nonlinear cubic polycrystals with perturbative and self-consistent methods [83].
In the following we outline a complete homogenizing procedure for two nonlinear composite materials that paradigmatically represent most features of the above described examples. Firstly, we consider a dispersion of nonlinear (but isotropic, i.e., amorphous or polycrystalline) spheres embedded in a linear homogeneous and isotropic matrix. The nonlinearity of the spheres can be described, at most, by four parameters (the so-called Landau coefficients) measuring the deviation from the linearity. Since the overall behavior of the heterogeneous structure will be elastically nonlinear, the key point is the evaluation of the effective nonlinear properties of the composite material. Secondly, a similar procedure has been developed for a distribution of parallel (nonlinear) cylinders embedded in a (linear) matrix. In both geometries, the most important methodological aspect is given by a useful generalization of the Eshelby theory [84] to nonlinear inhomogeneities.

Nonlinear Elastic Constitutive Equations
In geometrically linear elasticity, the balance of linear and angular momentum hold for all materials, regardless of their constitution.
The balance of the linear momentum leads to the equation of motion in the form where the T ij are the components of the Cauchy stress tensorT , b j are the components of the externally applied body force b, ρ is the mass density and u i are the components of the displacement u [81]. The balance of the angular momentum leads to the symmetry of the stress tensor (T ij = T ji ). However, these relations are generally insufficient to determine the elastic fields produced by given boundary conditions and body forces. They need to be supplemented by a further set of relations, referred to as constitutive equations, which characterize the constitution of the body. The convenient starting point is a set of relations in which the stress components are regarded as single-valued functions of the strain components where the functions f ij are chosen so that f ij = f ji in order to satisfy the stress symmetry andˆ is the small-strain tensor with components ij = 1 2 ∂u i ∂x j + ∂u j ∂x i [81]. For example, in physically linear elasticity the tensor relationship is taken into consideration in order to describe any kind of anisotropy. The stiffness tensor C ijkh must fulfill the following constraints: • C ijkh = C jikh and C ijkh = C ijhk , in order to preserve the symmetry of the stress tensor and of the strain tensor; • C ijkh = C khij , derived by the existence of an elastic energy density (Green hypothesis).
For isotropic materials the previous considerations lead to the notable stress-strain relationT = 2µˆ + λTr (ˆ )Î, which is based on the two independent Lamè constants λ and µ [85,86]. When considering nonlinear elastic material, Equation (86) can be generalized by means of higher order elastic moduli taking into account the deviation from the stress-strain proportionality [87,88] whereĈ N L (ˆ ) is the nonlinear (strain dependent) stiffness tensor. It can be noticed that the tensorĈ has 21 independent entries, while the second order tensorL has 56 independent components. Tables for the values of C ijkh and L ijkhnm can be found in literature [83]. These values can be obtained by experimental procedure [89,90] and by computational techniques (e.g., molecular dynamics [91] or first-principles calculations [92]). For the following purposes we are interested in the isotropic nonlinear constitutive equations expanded up to the second order in the strain components. In order to introduce these forms of physical nonlinearities we can take into account two different approaches, as described below.

Cauchy elasticity
The Cauchy approach to the constitutive equations is the less restrictive starting point for the elasticity theory since it does not consider the strain energy function. It is simply based on the Equation (85). To develop this approach in an isotropic context an assumption must be made concerning the behavior of Equation (85) for all proper orthogonal tensorR representing the rotation. A function satisfying the previous identity is known as an isotropic tensor function, and it can be represented in the form [81] whereÎ is the identity operator and q 1 , q 2 and q 3 are scalar functions of the invariants Tr(ˆ ), Tr(ˆ 2 ) e Tr(ˆ 3 ) of the strain tensorˆ q α = q α Tr(ˆ ), Tr(ˆ 2 ), Tr(ˆ 3 ) The development of Equation (89), up to the second order in the powers ofˆ , provides the following constitutive equation where µ and λ are the standard Lamè moduli concerning the linear contribution and A, B, C and D are the coefficients describing the nonlinear behavior of the material.

Green elasticity
The Green elasticity is based on Equation (85) with an additional hypothesis: we suppose that the stress power, in a given deformation, is absorbed into a strain energy function U (ˆ ), representing the density of elastic potential energy. The existence of such a function and the consideration of energy balance in the continuum, lead to the evolution equation [85] dU (ˆ ) dt = T ij (ˆ ) d ij dt (92) affirming that the function U (ˆ ) is an exact differential form such that So, if a function U (ˆ ) exists, the (arbitrarily nonlinear) constitutive equation for a given material can be determined by Equation (93) [85,86]. From the thermodynamics point of view, the strain energy function can be identified with the internal energy per unit volume in an isentropic process, or with the Helmholtz free energy per unit volume in an isothermal process. Such an approach can be further developed for isotropic media: in this case, the function U (ˆ ) must satisfy the relation [81] for any rotation tensorR. Equation (94) represents the scalar counterpart of the tensor relation Equation (88). If Equation (94) is true then it follows that the function U (ˆ ) can depend only on the principal invariants of the strain tensor We may expand Equation (95) up to the third order in the strain components, obtaining [86] U Finally, performing the derivatives indicated in Equation (93), we obtain the nonlinear isotropic constitutive equation (within the Green approach) expanded up to the second order in the strain tensor It is evident by comparison of Equation (91) and Equation (97) that the Green elasticity is more restrictive than the Cauchy elasticity: we obtain the Green formulation from the Cauchy formulation by imposing D = 2B. We use four independent parameters (A, B, C and D) in the Cauchy elasticity and three independent parameters (A, B and C) in the Green elasticity. These parameters are called Landau coefficients [86].

Eshelby Theory for Nonlinear Inhomogeneities
A nonlinear isotropic and homogenous ellipsoid can be generically described by the relationT = C (2) (ˆ )ˆ (see Equation (87)). Let us now place this inhomogeneity in a linear matrix characterized by a stiffness tensorĈ (1) (see Figure 17) and let us calculate the strain field inside the particle when a uniform fieldT ∞ =Ĉ (1)ˆ ∞ is remotely applied to the system.
If the particle were linear, withĈ (2) independent from the strain, we would have, inside the ellipsoid, a uniform strain fieldˆ s given by the Eshelby theory [93,94] where, we have introduced the Eshelby tensorŜ, which depends only on geometrical factors of the ellipsoid (the semi-axes a 1 , a 2 and a 3 ) and on the Poisson ratio of the host matrix [84]. Conversely, if the ellipsoid were nonlinear, it is easy to prove that the internal uniform field must satisfy the equation Figure 17. Scheme of an ellipsoidal inhomogeneity.
x 1 x 2 If a solutionˆ s * exists for a given ∞ , it means that the nonlinear inhomogeneity could be replaced by a linear one with constant stiffnesŝ , without modifications of the elastic fields at any point. Therefore, if a solution exists, then Equation (99) exactly describes, through self-consistency, the elastic behavior of the nonlinear anisotropic inclusion. This is not a trivial result: for instance, such a generalization of Equation (98) is not valid if a nonlinear behavior is assumed for material 1 (matrix). The calculation of the internal strain field from Equation (99) is very complicated and it strongly depends on the kind of nonlinearitŷ T =Ĉ (2) (ˆ )ˆ . This task will be accomplished in the following, dealing with a sphere or a cylinder described by physical nonlinearities as those in Equation (91) (Cauchy) or Equation (97) (Green).
To conclude, we have proved the following general statement: if the linear elastic space with a single inhomogeneity of ellipsoidal shape is subjected to remote uniform loading, the stress field inside the inhomogeneity will be uniform independent of the constitutive law for the inhomogeneity, provided that both the matrix and the particle are homogeneous bodies. Some similar properties can be found in earlier literature [95][96][97].
When the Green approach is considered it is also possible to verify the existence and the uniqueness for the solution of Equation (99) [98,99]. The proof follows.

Nonlinear Eshelby theory within green elasticity
We adopt here, from the energetic point of view, the Green formulation of the elasticity theory. A strain energy function U (ˆ ) defines the constitutive equationT (ˆ ) = ∂U (ˆ ) ∂ˆ of the inhomogeneity, which is equivalent toT (ˆ ) =Ĉ (2) (ˆ )ˆ . In these conditions, the existence and uniqueness of a solution for Equation (99) can be exactly proved under the sole hypothesis of convexity for the strain energy function U (ˆ ) [98,99]. To prove this statement, we rearrange Equation (99) as follows Now, the first linear term can be converted to the gradient of a quadratic form and the second constant term can be converted to the gradient of a linear form. At the end we observe that the internal strain field must satisfy the following relation [98,99] ∂ ∂ˆ which is exactly equivalent to Equation (99). The first term represents a symmetric and positive definite quadratic form inˆ (see Appendix A) while the second term is a linear function ofˆ . Therefore, the sum of these two terms is a convex functional with relative minimum at Î −Ŝ ˆ ∞ . This value represents the strain field in a void (Ĉ (2) (ˆ ) = 0 in Equation (99) or U (ˆ ) = 0 in Equation (101)) embedded in the matrix with stiffnessĈ (1) . If U (ˆ ) is a convex functional (with U (0) = 0) the brackets in Equation (101) contain the sum of two convex terms: they result in an overall convex functional with a unique minimal extremum atˆ s .

Elastic Dispersion of Nonlinear Spherical Inhomogeneities
We consider an assembly of spherical inhomogeneities (see Figure 18) described by a Cauchy constitutive relation randomly embedded in a linear matrix with stiffness tensorĈ (1) (moduli λ 1 and µ 1 ). We also introduce the bulk moduli K 1 = λ 1 + 2 3 µ 1 and K 2 = λ 2 + 2 3 µ 2 . If needed, we can easily move to the Green elasticity by assuming D = 2B. We suppose that the volume fraction c of the embedded phase is small (dilute dispersion). Since the elastic interactions can be neglected, each sphere behaves as an isolated one under the effect of a remote loadT ∞ =Ĉ (1)ˆ ∞ . The starting point for the evaluation of the induced internal strainˆ s is Equation (99), which can be usefully rearranged as followŝ Here, we have introduced the internal stress given by the relationT s =Ĉ (2) (ˆ s )ˆ s . The result of the application of Ĉ (1) −1 over the stress tensorT s can be easily written in explicit form Tr T s Î (104) Figure 18. Scheme of a dispersion of nonlinear spheres embedded in a linear matrix.
Moreover, the explicit expression of the Eshelby tensor for a sphere is reported in literature [28,84] We can evaluate the effect of S ijkh over an arbitrary strain s kh , getting Now, the Poisson ratio ν 1 of the matrix can be written in terms of the bulk modulus K 1 and the shear modulus µ 1 through the standard relation ν 1 = 3K 1 −2µ 1 2(3K 1 +µ 1 ) , obtaininĝ Sˆ s = 6 5 In order to find a single equation for the internal strainˆ s , we can substitute Equations (102), (104) and (107) in Equation (103). A long algebraic calculation leads to the important equation [98,99].
which completely defines the internal strain induced in a nonlinear sphere by the uniform remote defor-mationˆ ∞ . The parameters L, M, N, O, P and Q have been written in terms of the shear moduli, bulk moduli and nonlinear coefficients as follows L = 1 + 6 5 At this point we take into consideration the actual dispersion of spheres. We define V as the total volume of the composite material, V e as the volume corresponding to the spheres and V o as the volume of the matrix (V = V o ∪ V e , see Figure 18). Since we are working under the hypothesis of small volume fraction c, we can consider the average value of the strain in the matrix equal to the externally applied strainˆ ∞ . Therefore, the average value of the strain in the overall system is given by On the other hand, the average value of the stress over the entire structure can be calculated as follows In order to obtain the macroscopic characterization of the material, we search for the relationship between T and ˆ , given in Equations (115) and (116), respectively.
By substituting Equation (108) in Equation (115), we obtain the average strain ˆ in terms of the internal strainˆ s and by substituting the constitutive relations in Equation (116), we obtain the average stress T in terms ofˆ s The last two expressions, although in implicit form, define the macroscopic constitutive equation relating T and ˆ . In fact, we may obtainˆ s in terms of ˆ from Equation (117) and this result can be replaced in Equation (118), leading to the final characterization. In order to follow this scheme, we rewrite Equation (117) in a simpler form where we have used the definitions Starting from Equation (119), we can straightforwardly calculate the quantities Tr ˆ , ˆ 2 , ˆ Tr ˆ , Tr ˆ 2 and [Tr ˆ ] 2 in terms of the internal strainˆ s (by using the relation Tr(Î) = 3). These set of relations can be written neglecting the terms of order greater than two inˆ s , since we are interested in the characterization of the nonlinear elastic properties of the dispersion up to the second order. Therefore, this set of equations can be arranged in a matrix form, as follows The elements of the matrixŨ have been written in terms of the parameters defined in Equations (120)-(125) Finally, by using Equation (118) and by inverting Equation (126), we may obtain the matrix form of the complete constitutive relation

Results
The constitutive equation in the form of Equation (128) can be written in terms of the effective linear and nonlinear elastic moduli as follows As for the linear elastic moduli, we obtain and, as for the nonlinear elastic moduli, we have [98,99] If we use the definitions of the parameters L and M , given in Equations (120) and (121), we obtain for the effective shear and bulk moduli the explicit expressions [98,99] These two expressions are coincident with those obtained for a linear dispersion of elastic spheres [30]. However, they are completed by Equations (132) 2. The nonlinear elastic moduli A, B, C and D influence the effective nonlinear moduli of the composite material following the universal scheme showed in Figure 19. Therefore, there is a complicated mixing of the nonlinear elastic modes induced by the heterogeneity of the structure. The results for the nonlinear effective parameters, obtained with a single coefficient (A, B, C or D) different from zero, are reported in Appendix B, in form of series expansions in the volume fraction up to the first order. They are coherent with the scheme shown in Figure 19. 3. If the linear elastic moduli of the matrix and of the spheres are the very same (K 1 = K 2 and µ 1 = µ 2 ), we simply obtain K ef f = K 1 , µ ef f = µ 1 and the following special set of results for the nonlinear components It means that the nonlinearity of the overall system is simply proportional to the nonlinearity of the spheres.

4.
We have developed our procedure under the hypothesis of Cauchy nonlinear elasticity for the spheres embedded in the linear matrix. If we let D = 2B we move from the Cauchy elasticity to the Green elasticity, assuming the existence of a strain energy function for the inhomogeneities. It is important to remark that the following property holds: if D = 2B then the relation D ef f = 2B ef f is true for the effective nonlinear moduli. It can be verified by direct calculation and it means that our approach is perfectly consistent with the energy balance of the composite material. In other words, we have verified that if a strain energy function exists for the embedded spheres, then an overall strain energy function exists for the whole composite structure.
5. If we consider the special value of the Poisson ratio ν 1 = ν 2 = 1/5 (both for the matrix and the spheres) and different values for the Young moduli E 1 = E 2 , we obtain another interesting result: the effective Poisson ratio assume the same value ν ef f = 1/5, the effective Young modulus E ef f assumes the value and the effective nonlinear elastic moduli can be calculated as follows where the symbol X represents any modulus A, B, C or D (the four effective parameters exhibit the same behavior). Therefore, we can say that the special value ν 1 = ν 2 = 1/5 uncouples the behavior of the nonlinear elastic modes (described at the point 2), generating a direct correspondence among the nonlinear moduli of the spheres and the effective nonlinear moduli. Furthermore, if we add the condition E 1 = E 2 , we get back to the point 3. The special value 1/5 for the Poisson ratio comes out in several issues considering a dispersion of spheres. For example, for linear porous materials (with spherical pores) and for linear dispersions of rigid spheres the value 1/5 is a fixed points for the Poisson ratio: if ν 1 = 1/5, then we have ν ef f = 1/5 for all spheres concentrations [33,100]. Moreover, there is another interesting behavior of the effective Poisson ratio for high volume fraction of pores or rigid spheres: in both cases for c → 1 the effective Poisson ratio converges to the fixed value ν ef f = 1/5, irrespective of the matrix Poisson ratio [33, 100-102].
6. Finally, we analyze the properties of the dispersion when incompressible material is utilized for the embedded spheres: the constitutive relation Equation (102) describes an incompressible medium in the limit λ 2 → ∞ (or, equivalently, K 2 → ∞ since K 2 = λ 2 + 2µ 2 /3); by inverting Equation (102), writing the strain tensor in terms of the stress tensor and performing such a limit, we obtain (up to the second order) s Tr T s which describes a nonlinear isotropic and incompressible material. We remark that only the nonlinear modulus A intervenes in defining such a constitutive equation and that Equation (144) imposes Tr (ˆ s ) = 0, as requested by the incompressibility. In this limiting condition, as for the effective linear moduli, we observe that Equation (136) for µ ef f remains unchanged and Equation (137) leads to On the other hand, the nonlinear elastic moduli have been eventually found as One can observe that, as expected, the effective nonlinear elastic moduli depend only on the modulus A describing the nonlinearity of the spheres, as shown in Equation (144). Moreover, we remark that a single modulus A for the spheres can generate four different effective nonlinear moduli, as predicted by the scheme in Figure 19.
To conclude, we present some numerical results obtained by the implementation of Equations (132)-(137). In Figure 20 we have considered Green nonlinear elasticity and the mixture parameters: µ 1 = 1, µ 2 = 4, K 1 = 7, K 2 = 1, A = 2, B = 3, C = 5, D = 2B in arbitrary units. In Figure 21 we have considered Cauchy nonlinear elasticity and the mixture parameters: µ 1 = 1, µ 2 = 4, K 1 = 10, K 2 = 1, A = 2, B = −3, C = −5, D = 4 in arbitrary units. The results have been presented in terms of the volume fraction c of the spheres. In both cases we may observe a consistent amplification of the nonlinear effective modulus C ef f . We have verified that such a phenomenon is always exhibited when , when the matrix is much more incompressible than the spheres) and that the higher values of C ef f appear for small values of the volume fraction c, belonging to the range of applicability of the present theory.
As it is well known, simple limitations for the values of the linear effective moduli are well established The lower bounds in Equations (152) and (153) are referred to as the Voigt bounds, and the upper bounds are designated as the Reuss bounds [29]. Unfortunately, these bounds are of no practical value, but more refined bounds, with realistic applications, have been derived by Hashin and Shtrikman [4]. From our numerical results, shown in Figures 20 and 21, we may observe that the nonlinear properties, contrary to the linear ones, are not bounded by some given values and in certain conditions exhibit a strong amplification, which leads to nonlinear effective moduli much greater than those of the constituents. This point is important in the topic of designing materials with desired properties and functions.

Elastic Dispersion of Parallel Nonlinear Cylindrical Inhomogeneities
We now take into consideration an assembly of parallel cylinders, as represented in Figure 22, described by an arbitrary Cauchy constitutive relation [see Equation (102)]. As before, when needed, we can easily move to the Green elasticity by assuming D = 2B. The cylindrical inhomogeneities are randomly embedded in a linear matrix with elastic moduli K 1 and µ 1 . This is a simple but complete way for modeling a nonlinear fibrous material. In earlier works the linear analysis for a parallel distribution of fibers has been developed by means of the Eshelby methodology and of the differential effective medium theory [32,103]. Moreover, the mechanical response of elastic and inelastic fiber-strengthened materials has been investigated, also with self-consistent models [104][105][106]. Here, in order to deal with the nonlinear properties, we suppose that the volume fraction c of the embedded phase is small (dilute dispersion). It means that each cylinder can be considered isolated in the space (not interacting with other inhomogeneities) and subject to the same external loading. In order to simplify the modeling and considering that the system shows a transverse isotropic symmetry (otherwise said uniaxial symmetry), we assume the plane strain condition on an arbitrary plane π (see Figure 22) orthogonal to the cylinders. It means that we are dealing with a problem belonging to the two-dimensional elasticity. Moreover, in plain strain condition, it is a common choice to introduce the two dimensional elastic moduli µ 2D = µ and K 2D = K+µ/3, where K and µ are the customarily used three-dimensional moduli [103]. Throughout this section we indicate for brevity K and µ alluding to the two-dimensional version of the elastic moduli. It means that the linear matrix is described bŷ and the cylindrical inhomogeneities are described by the Cauchy constitutive relation where any strain or stress tensor is represented by a square matrix of order two, working in the framework of the two-dimensional elasticity. Now, we remark that Equation (99) or, equivalently, Equation (103) are correct for any geometry and, therefore, they can be directly used in the present analysis. Nevertheless, in order to use Equation (103) we need to consider some ingredients: the result of the application of the compliance tensor of the matrix on the stress tensorT s can be written as Moreover, the effect of the Eshelby tensorŜ for a cylinder over an arbitrary strain tensorˆ s is given by [84]Ŝˆ s = 1 2 We follow a procedure similar to that described in Section 4. We use again Equations (116) and (119) for the average values of the stress and the strain over the whole composite material. At this point, starting from Equation (119), we obtain the system given in Equation (126) (by using the two-dimensional relation Tr(Î) = 2). It is defined by the matrixŨ where the elements depend on the parameters defined in Equations (120)-(125) and calculated by means of Equations (159)-(164) Finally, by inverting Equation (126), we may obtain the matrix form of the complete constitutive relation

Results
The constitutive equation in the form of Equation (166) can be written in terms of the effective linear and nonlinear elastic moduli as follows As for the linear elastic moduli, we obtain [99] It is important to remember that the bulk modulus K ef f represents the two-dimensional version, as above defined. Moreover, the two linear results given in Equations (168) and (169) are perfectly coincident with earlier literature [107]. As for the effective nonlinear elastic moduli, we have the following final results [99] A ef f = Ac They represent the complete nonlinear characterization of the random dispersion of parallel cylinders. It is interesting to observe that all the properties described in the previous section for the dispersion of spheres (points 1-6) can be easily verified also for the present case [99]. In particular, the scheme represented in Figure 19 remains valid. In Appendix C we have reported the explicit results giving the first order expansions of the nonlinear elastic moduli with respect to the volume fraction, corresponding to the simple cases where only one nonlinear parameter of the cylinders is different from zero. We analyze the case corresponding to the point 5 of the previous section: we consider the special value of the three-dimensional Poisson ratio ν 1 = ν 2 = 1/4 (corresponding to the two-dimensional Poisson ratio ν 2D = ν 3D /(1 − ν 3D ) = 1/3 [103]) and different values for the three-dimensional Young moduli E 1 = E 2 . In this case, the effective 3D Poisson ratio assume the value ν ef f = 1/4 and the effective 3D Young modulus E ef f assumes the value Moreover, the effective nonlinear elastic moduli can be calculated as follows where the symbol X represents any modulus A, B, C or D (the four effective parameters exhibit the same behavior). Therefore, as before, we can say that the special value ν 1 = ν 2 = 1/4 uncouples the behavior of the nonlinear elastic modes, generating a direct correspondence among the nonlinear moduli of the spheres and the effective nonlinear moduli. Finally, we have numerically implemented Equations (168)-(173) in order to show some explicit results. In Figure 23 we have considered Green nonlinear elasticity and the mixture parameters: µ 1 = 1, µ 2 = 5, K 1 = 10, K 2 = 1, A = −8, B = −2, C = −1, D = 2B in arbitrary units. In Figure 24 we have considered Cauchy nonlinear elasticity and the mixture parameters: µ 1 = 1, µ 2 = 5, K 1 = 10, K 2 = 1, A = 8, B = −2, C = −1, D = 6 in arbitrary units. As in the previous section, we may observe a consistent amplification of the nonlinear effective modulus C ef f . We have also verified that such a phenomenon is exhibited when K 1 K 2 (i.e., when the matrix is much more incompressible than the spheres) and that the higher values of C ef f appear for small values of the volume fraction c, belonging to the range of applicability of the present theory.

Conclusions
In the previous sections we have considered the linear and nonlinear elastic behavior of a composite material. In particular we have taken into account a dispersion of isotropic nonlinear inhomogeneities (spheres or parallel cylinders) embedded into a linear isotropic host matrix. The nonlinearity of the inhomogeneities has been described either by the Cauchy model (four parameters) or by the energybased Green approach (three parameters).
We have introduced two simplifying hypotheses: the small volume fraction of the embedded particles and the small deformations of the whole solid body. Nevertheless, we have described useful results both for analyzing the mechanical properties of a given heterogeneous structure and for designing a composite material with a desired linear and nonlinear elastic behavior.
The main concept introduced to homogenize the heterogeneous structures is a generalization of the linear Eshelby methodology developed for extending its applicability to nonlinear materials. This approach has been analytically applied to perform a linear and nonlinear micromechanical averaging in the composite structure and, therefore, to develop a complete homogenizing procedure yielding the mechanical behavior of the solid body at the macro-scale.
As for the linear properties, we have obtained a series of results in perfect agreement with earlier researches on this subject. This point can be considered as a check of the mathematical procedure. As for the nonlinear properties, firstly, we have obtained the expressions of the four effective elastic moduli of the composite medium with inhomogeneities described by the Cauchy constitutive equations, which represent the less restrictive way to model the nonlinear elasticity. Then, we have considered, as a particular case, the Green elasticity to describe the nonlinear behavior of the particles. In this case we have verified that if a strain energy function exists for the inhomogeneities, then an overall strain energy function exists for the whole composite structure. This point confirms the perfect coherence between our micromechanical averaging procedure and the thermodynamics of the composite material.
Moreover, we have observed that the nonlinear effective elastic moduli, contrary to the linear ones, are not subject to specific bounds that limit their values when the behaviors of the constituents are chosen. We have indeed found some strong amplifications of the nonlinear behavior in certain given conditions. More specifically, for example, we have observed that the nonlinear modulus C ef f can assume values much greater than C if the matrix is much more incompressible than the inhomogeneities. This is a crucial point that can be applied in analyzing and designing composite materials with a given microstructure.
Finally, some special values of the Poisson ratio of the materials have been found in order to obtain a direct correspondence among the nonlinear moduli of the inhomogeneities and the effective moduli of the composite structure. It means that, under the above conditions, we can realize a perfect scaling of the nonlinear properties (see Equation (143) or (175)) modulated by the ratio E 1 /E 2 between the Young moduli of the constituents.
Appendix A Symmetry and Positive Definiteness of the TensorĈ (1) Ŝ −1 −Î We briefly outline the concepts of inclusion and linear inhomogeneity in order to present the adopted notation and to recall the most important equations of the Eshelby theory [84,93,94].

Concept of inclusion
We consider an infinite medium with stiffness tensorĈ (1) ; moreover, we consider an embedded ellipsoidal inclusion V described by the constitutive equationT =Ĉ (1) (ˆ −ˆ * ). The strainˆ * is called eigenstrain (or stress-free strain). In these conditions the following relations describe the strain inside and outside the inclusion [84]ˆ whereŜ is the internal Eshelby tensor andŜ ∞ is the external Eshelby tensor.

Concept of inhomogeneity
We now consider an infinite medium with stiffness tensorĈ (1) in 3 V (matrix) andĈ (2) in the ellipsoidal region V (inhomogeneity). We remotely load the system with a uniform strainˆ ∞ or, equivalently, with the uniform stressT ∞ . Of course we haveT ∞ =Ĉ (1)ˆ ∞ . This configuration can be analyzed by means of the Eshelby equivalence principle [93]. The system can be described by the superimposition of two simpler cases (see Figure 25) [84]. The first situation A concerns a medium with stiffnessĈ (1) Figure 25. Scheme of an ellipsoidal inhomogeneity and the Eshelby equivalence principle.
ijkh kh (without inclusions or inhomogeneities) uniformly deformed by means of the remote loadsˆ ∞ or T ∞ . The second situation B is represented by an inclusion embedded in a medium, characterized everywhere byĈ (1) and having an eigenstrainˆ * in V. The situation B is without remote loads. The eigenstrain must be imposed searching for the equivalence between the original inhomogeneity problem and the superimposition A + B. The following relation hold on inside the region V (s means inside V) In the inhomogeneity we haveT s =Ĉ (2)ˆ s and thereforê The following relations can be finally obtained for the eigenstrain and for the actual strain in V * = Î − Ĉ (1) −1Ĉ ( We consider the same inclusion V with two different values for the eigenstrainˆ * andˆ * * embedded in the material defined byĈ (1) . The symmetry of the tensor can be established by means of a revised version of the Betti's reciprocal theorem [85]. We defineT * =Ĉ (1)ˆ * andT * * =Ĉ (1)ˆ * * . The first situation is described by the fieldsT ,ˆ , u and the second one byT ,ˆ , u everywhere in the space.
The preliminary symmetry of the tensorŜ Ĉ (1) −1 is proved. We begin by considering the following relation (V is the inclusion volume , Σ its boundary and n its external normal unit vector) At the interface Σ we haveT n| Σ − =T n| Σ + (sign + indicates the external side of Σ and sign − indicates its internal side). Recalling the definition of inclusion we simply obtainĈ (1) (ˆ −ˆ * ) n| Σ − =Ĉ (1)ˆ n| Σ + and finally we getĈ (1)ˆ n| Σ − −Ĉ (1)ˆ n| Σ + =Ĉ (1)ˆ * n. We use it in Equation (183), obtaining On Σ − we haveT =Ĉ (1) (ˆ −ˆ * ) and on Σ + we haveT =Ĉ (1)ˆ , therefore We have now obtained a symmetric form (sinceĈ (1) is symmetric). Therefore, the following dual relation is valid and it can be verified as above . We consider two similar situations as described in Figure 26. The first deals with an homogeneous medium with displacement prescribed on the boundary, while the second case considers the addition of an inhomogeneity without changing the fixed displacements on the external surface. No body forces are present in both schemes. We begin searching for the difference between the elastic energy stored in the two cases We simply verify that where ∂ u a ∂ xT a = ∂ u aTa ∂ x since ∂T a ∂ x = 0 at equilibrium and similarly ∂ u b ∂ xT a = ∂ u bTa ∂ x . Therefore, we obtain Since the stiffness tensorĈ (1) is symmetric, we obtain the following general expression for the energy difference We now suppose that the prescribed displacement on Σ imposes a uniform strain in the first case of Figure 26; therefore, the second situation can be described by the Eshelby solution. With this additional It is interesting to observe that all the results given in Appendix A can also be exactly applied to an anisotropic and homogeneous ellipsoidal inhomogeneity embedded in an anisotropic and homogeneous matrix. In this case, the Eshelby tensorŜ depends on the geometry and onĈ (1) [84,98,99].

Appendix B First Order Expansions for a Dispersion of Spheres
In this Appendix we present the first order expansions in the volume fraction of the effective nonlinear moduli A ef f , B ef f , C ef f , and D ef f for a dispersion of spheres. In particular we consider four different cases where only one nonlinear modulus of the spheres (A, B, C or D) is different from zero. These solutions are coherent with the scheme represented in Figure 19.