On the effectiveness of the Moulinec–Suquet discretization for composite materials

Moulinec and Suquet introduced a method for computational homogenization based on the fast Fourier transform which turned out to be rather computationally efficient. The underlying discretization scheme was subsequently identified as an approach based on trigonometric polynomials, coupled to the trapezoidal rule to substitute full integration. For problems with smooth solutions, the power of spectral methods is well‐known. However, for heterogeneous microstructures, there are jumps in the coefficients, and the solution fields are not smooth enough due to discontinuities across material interfaces. Previous convergence results only provided convergence of the discretization per se, that is, without explicit rates, and could not explain the effectiveness of the discretization observed in practice. In this work, we provide such explicit convergence rates for the local strain as well as the stress field and the effective stresses based on more refined techniques. More precisely, we consider a class of industrially relevant, discontinuous elastic moduli separated by sufficiently smooth interfaces and show rates which are known to be sharp from numerical experiments. The applied techniques are of independent interest, that is, we employ a local smoothing strategy, utilize Féjer means as well as Bernstein estimates and rely upon recently established superconvergence results for the effective elastic energy in the Galerkin setting. The presented results shed theoretical light on the effectiveness of the Moulinec–Suquet discretization in practice. Indeed, the obtained convergence rates coincide with those obtained for voxel finite element methods, which typically require higher computational effort.

a computational method to solve the cell problem of homogenization which is based on the fast Fourier transform (FFT) as an efficient alternative to more classical numerical approaches, for example, based on a finite element discretization.
In its original form, FFT-based computational micromechanics involved essentially three ingredients-a class of models to be treated, a discretization and a solver. The original Moulinec-Suquet method was designed to handle quasistatic mechanical problems at small strains. In particular, inelastic problems were in the focus right from the beginning. In a continuous effort, researchers extended the domain of applicability, for example, towards problems at finite strains [6][7][8] or coupled problems. [9][10][11] We refer to the review articles [12][13][14] for a glimpse at the numerous applications of this class of numerical methods.
These developments were complemented by extensions of the original solution scheme of Moulinec-Suquet which they termed the basic scheme. The latter was based on the weak-contrast expansion well-known in micromechanics 15,16 which itself derives from series representations 17,18 of effective properties in the material contrast established with the help of the Lippmann-Schwinger equation. [19][20][21] A key advantage of the basic scheme is its low memory demand-there are implementations which require only a single strain field in memory. 22,23 Over the years, faster solution methods were introduced, usually at the cost of a higher memory footprint. Eyre and Milton 24 introduced a numerical method similar to the basic scheme but with higher convergence rate. To be more precise, the method has a link to the strong-contrast expansion 25,26 and may be written as a series expansion in the square root of the contrast, which is responsible for the faster convergence behavior. The original method was restricted to linear problems. Subsequently, modifications of the Eyre-Milton method were introduced [27][28][29] which applied to inelastic materials in a hardening regime 30,31 or problems of fracture mechanics. [32][33][34] The caveat of these so-called polarization methods is the need to evaluate the constitutive law in a nonstandard way, and only for specific cases tricks are known to reduce the inherent effort considerably. 35,36 For this reason, also more traditional computational techniques were made available in the context of FFT-based methods, like Newton's method, 6,7,37 the (linear and nonlinear) conjugate gradient method 38,39 and a list of Quasi-Newton methods. [40][41][42] In particular, users of FFT-based methods may select a solver appropriate to their task at hand depending on the available resources and their performance demands, see the review article 14 for background and recommendations.
The extensions of applicability and solver technology were complemented by developments in discretization methods. To be more precise, there were essentially two lines of thought which were followed. On the one hand, the original Moulinec-Suquet method led to stress and strain fields which showed characteristic ringing artifacts in the vicinity of material interfaces. This so-called Gibbs phenomenon 43 is well-known in Fourier methods and reflects the limited ability of classical Fourier series to represent discontinuous functions. To remove these oscillations, different discretization methods were exploited, for example, finite difference, [44][45][46] finite element, [47][48][49] and finite volume 50,51 methods. The underlying key observation is that translation-invariant stencils of finite difference operators may be block-diagonalized in Fourier space, thus giving rise to a simple preconditioning strategy reminiscent of the original Moulinec-Suquet method. In particular, the developed solver technology could be transferred to this extended setting essentially without change. As an added benefit, it was observed that problems involving pores [52][53][54] or rigid inclusions 55,56 could be handled by dedicated discretization schemes which was difficult for the original Fourier-type discretization due to an intrinsic numerical instability related to the global support of the ansatz functions. 57 The second line of thought on discretization methods involved understanding the discretization underlying the original Moulinec-Suquet method. This was motivated, on the one hand, by the idea that understanding is key to improvement, and, on the other hand, by the desire to create methods which provide upper and lower bounds on the effective properties. To obtain bounds, the classical strategy to minimize the average strain energy over compatible strain fields leads to the Fourier-Galerkin method [58][59][60][61] when using trigonometric polynomials as the ansatz space. A similar strategy may be applied to minimize the complementary energy over self-equilibrated stress fields, 62 leading to computable lower bounds 59,61 of the effective elastic energy. Alternatively, the variational principle of Hashin and Shtrikman 63-65 may be used to compute upper and lower bounds to the effective elastic properties. Brisard and Dormieux 22,66 pioneered such an approach based on voxel-wise constant strains, later extended to B-spline ansatz functions. 67 Both variational principles were used as the point of departure for interpreting the original Moulinec-Suquet discretization. Brisard and Dormieux regarded the original scheme as a nonconforming version of their Galerkin method of Hashin-Shtrikman type and gave the first convergence proof 66 for the former. Their proof relied upon voxel-wise homogeneous stiffnesses that were obtained by a specific averaging method that inspired the composite voxel technology. [68][69][70] Vondřejc et al. 71 introduced another interpretation of the Moulinec-Suquet discretization as an underintegrated Fourier-Galerkin method where the underintegration arises from applying the trapezoidal rule. They also provided a convergence proof which-in contrast to Brisard and Dormieux 66 -required the coefficients of elasticity to be sufficiently smooth. These rather restrictive assumptions were subsequently removed. More precisely, convergence was established for merely Riemann integrable, that is, measurable and bounded, coefficients of elasticity. Moreover, a class of nonlinear elastic constitutive laws 72 was possible to be treated, all without the need for any averaging close to the interface.
Several years after these results, Ye and Chung 73 provided a fresh perspective on the convergence analysis of FFT-based computational micromechanics. They established convergence of the Moulinec-Suquet discretization, the discretization on a rotated staggered grid 46 and for trilinear finite elements. 47,48 Their analysis of the Moulinec-Suquet discretization involved an appropriate smoothing operator, allowing them to incorporate essentially arbitrary composite-voxel methods. [68][69][70] In the same paper, they provided a convergence proof for the discretization on a rotated staggered grid pioneered by Willot 46 based on an ingenious coordinate transformation in Fourier space which emphasizes the connections to the non-conforming P1-quadrilateral elements introduced by Park and Sheen. 74 In addition to the mentioned qualitative convergence statements, Ye and Chung 73 also provided quantitative estimates for the convergence rate of a finite element discretization on a regular grid. More precisely, they showed that the strain field converges as N −1∕2 in L 2 where N denotes the number of voxels per edge. A similar convergence rate is known when approximating outer boundaries with Cartesian finite elements, 75 so that the rate N −1∕2 did not come as a surprise for methods which approximate inner boundaries on regular grids. However, Ye and Chung 73 also provided the insight that the effective stresses converge as 1∕N, that is, with the quadratic rate of the local fields.
Based on Ye-Chung's observation, it was shown that this phenomenon of higher order convergence of effective stresses also holds for inelastic problems and Galerkin discretizations. 76 Moreover, it was shown that a clever averaging procedure may increase the accuracy of the computed effective stresses significantly 77 by using the insights of Ye and Chung 73 appropriately.

Contributions
The starting point of this article was the Ye-Chung superconvergence result 73 for the effective stress and its inelastic generalization. 76 In its simplest form, the higher-order rate for the effective stresses follows directly from Galerkin orthogonality. It was empirically observed by many authors that non-conforming discretizations, in particular based on underintegration, typically outperform Galerkin discretizations. However, their analysis is more involved, as the effect of inexact integration needs to be accounted for, see, for example, Ciarlet's monograph 78( §4.1) . As already discussed, the Moulinec-Suquet discretization 4,5 arises quite naturally as a Fourier-Galerkin method with trapezoidal integration. 71 The quadrature error of the trapezoidal rule is closely linked to the error arising from trigonometric interpolation. 79 Even simpler than trigonometric interpolation is trigonometric approximation, that is, the best approximation of a given function or field by trigonometric polynomials of a given degree in an appropriate norm. 80 Therefore, it seems reasonable to understand the approximation of the solution to the cell problems first.
The difficulty is the following. Spectral methods, in particular involving trigonometric polynomials, are extremely effective for problems with smooth solutions, as trigonometric polynomials approximate such functions rapidly. 81 Much less is known for problems with discontinuous solutions like the strain field solving the cell problem with discontinuous coefficients of elasticity.
We will briefly describe the idea to leverage the theory developed for the smooth setting in the following. The impatient reader may skip these details.
Suppose we wish to estimate the following quantity || − P N || L 2 , that is, the difference of the strain field and the trigonometric best approximation, encoded by the trigonometric projector P N . Suppose the strain field is sufficiently smooth in the different phases of the considered multi-phase material. Moreover, assume that the strain field is bounded and satisfies an H 1 bound in each phase, that is, the strain gradient is square integrable. Introduce a suitable smooth cutoff function N (see Section 2.2 for details) which is zero close to the interface and equal to unity away from the interface. Then, we may decompose the approximation error For the first term, which measures the cut-off error, we observe with a generic constant C. For the second term, by classical convergence rates for trigonometric approximation 81(Lemma 8.5.1) , we deduce with a generic constant C that may change from line to line. Here, we used the product rule and that the term ||∇ N || L 2 grows as N 1∕2 . Moreover, H 1 ± stands for the "broken" norm over all phases. With the sketched argument, we showed that the strain in the Fourier-Galerkin discretization converges as N −1∕2 provided the assumed regularity assumptions on the solution of the cell problem hold (see Section 3.1 below). By the Ye-Chung argument, 73 valid for a Galerkin discretization, we thus obtain an 1∕N-convergence rate for the effective stiffness.
Even this argument for the convergence behavior of the Fourier-Galerkin discretization is not published to the best of our knowledge. However, the article at hand is more ambitious: we show an N −1∕2 convergence rate for both the local strain and stress fields in L 2 and a 1∕N convergence rate for the effective stiffness and the Moulinec-Suquet discretization under practical assumptions on the discontinuous stiffness distribution.
Let us give a brief overview of the challenges and the strategies to overcome them. Classical convergence estimates for trigonometric interpolation 81 require sufficiently high Sobolev regularity. Indeed, these assumptions ensure Sobolev embedding into a Hölder class, ∈ H k for k > d∕2, where d denotes the spatial dimension, and evaluating the Fourier series pointwise makes sense. Unfortunately, the classical regularity estimates are not sufficient in the physically relevant dimensions 2 and 3, as they only provide ∈ H 1 in the phases. There are more sophisticated convergence estimates for trigonometric interpolation which apply to general Riemann integrable functions. 82,83 These estimates are formulated in terms of so-called averaged moduli of smoothness. 79 Translated back into the scale of Sobolev spaces, we would be led to Lebesgue norms of higher mixed derivatives (up to order d in d spatial dimensions), which are known to be sharp. Thus, the basic problems remains also for these finer tools.
To overcome the problem with higher derivatives, we use Bernstein estimates which relate the L p -norm of the gradient of a trigonometric polynomial to the L p -norm of the polynomial for 1 ≤ p ≤ ∞ at the expense of a prefactor that depends on the order of the trigonometric polynomial, see Nikol'skii 84( §2.5) for an exposition. By a careful analysis, we are able to discard higher-order derivatives as necessary and retain only the derivatives where bounds are available in the continuous setting. Based on this strategy, we provide the desired convergence-rate estimate for the local strain and stress fields in Sections 3.3 and 3.4, respectively. Due to the interpolatory nature of the Moulinec-Suquet discretization, the convergence rate for the stresses does not immediately follow from the corresponding result of the strains, and an extra argument is required.
To provide a sharp convergence analysis of the effective stiffness, in addition to the Ye-Chung argument, 73 we must control the quadrature error. However, when inspecting the error of the trapezoidal rule, we run into essentially the same problem as for trigonometric interpolation: elementary estimates, see Section 2.3, show that mixed derivatives of high order for ∈ N d with || || ∞ = 1 must be integrable. Using Bernstein estimates allows us to handle this case as well, see Section 3.5.
We close this section by brief remarks on organization, notation and other simplifications. Section 3 represents the core of this article, as it contains the arguments for the announced convergence rates. To increase readability, the necessary mathematical tools were out-sourced into Section 2, where characteristics of trigonometric polynomials, the delicate construction of the cut-off function and an elementary error estimate for the trapezoidal rule are collected. The sharpness of the developed theory is demonstrated via computational experiments, see Section 4.
Concerning notation, we use the ≲-symbol to replace the expression ≤ C with an implicit constant C to increase readability of the article. The constant depends on the spatial dimension, the lower and upper bounds (55) of the stiffness and the interface. Sometimes, we will also make the dependence on norms of the solution field implicit. Last but not least, let us discuss simplifications made for the sake of exposition. We will restrict to a two-phase decomposition of the microstructure, although any finite number of phases may be treated. Also, we assume the stiffness to be phase-wise constant, although a sufficiently smooth variation could be handled, as well. We will not take into account composite voxel methods, although they do not deteriorate the convergence rate as long as they respect the bounds on the stiffness. We will only consider cubic cells and suppose that the number of voxels N in each direction is identical. Cuboidal cells and different boxel sizes 85 in different directions can be handled at the expense of simplicity of notation and exposition.
Last but not least let us remark that we restrict to odd voxel numbers N to avoid discussing issues with the Nyquist frequency. This is actually only critical when bounds are desired in the Fourier-Galerkin method. 59 In practice, one may either set the Nyqvist frequency of the strain or the stress field to zero. Again, to declutter the presentation, we will omit these details.

Fourier series and trigonometric polynomials
We are concerned with the cubic cell in d spatial dimensions. In this manuscript N will always refer to an odd natural number, which controls the level of discretization. The "discrete" version of the cell (4) will be denoted as where the typical element x I has the form By a trigonometric polynomial of order N on the cell Y d we denote a function f N ∶ Y d → R on the cell Y d of the form with suitable coefficientsf ( ) that satisfy the relationŝ which ensure that the function f N in the representation (7) is real-valued. To keep the notation concise, we introduced the set of restricted Fourier frequencies.
We denote the (finite-dimensional) vector space of trigonometric polynomials (7) by  N (Y d ). In addition to scalar-valued trigonometric polynomials (7), we will frequently need vector-and tensor-valued trigonometric polynomials, as well. The statements in this section extend to trigonometric polynomials valued in a finite-dimensional Euclidean space V in a straightforward way. We denote the latter space by  N (Y d ; V), and will consider V = R d and V = Sym(d), the space of symmetric d × d-tensors, in this work.
Parseval's identity 86( §7) asserts that the equality holds for any trigonometric polynomial f N ∈  N (Y d ), where − ∫ denotes the mean integral, that is, the integral multiplied by the inverse of the d-volume of the set. By polarization, the identity (10) also holds for the product of two trigonometric that is, the L 2 inner product between two trigonometric polynomials may either be computed by integrating in real space or by summing in Fourier space. Notice that we will always scale Lebesgue norms in such a way that mean integrals are used.
Last but not least it is also possible to describe trigonometric polynomials via interpolation. Suppose discrete values are given. Then, there is a unique polynomial f N ∈  N (Y d ) interpolating these values, that is, the condition holds for all indices I and where the points x I are defined in Equation (5). The solution to the problem (13) is constructive and given by the discrete Fourier transformf As a direct consequence, a discrete version of the product formula (11) follows immediately. Combining the two formulas (11) and (15) yields the statement which asserts that products of trigonometric polynomials of order N may be integrated exactly by the trapezoidal rule, see Section 2.3 below.
In the article at hand, we are interested in the approximation quality of trigonometric polynomials. Any square-integrable function f ∈ L 2 (Y ) admits a Fourier-series expansion 87 (17) to be understood in the L 2 -sense. Let us denote by P N f the trigonometric polynomial which arises from truncating the Fourier coefficients. The trigonometric polynomial P N f represents the unique L 2 -best approximation of the function f by elements of the space  N (Y d ) of trigonometric polynomials. The association f  → P N f gives rise to a bounded linear operator P N on the space L 2 (Y d ) which satisfies the projector identity P N P N = P N . A classical result of Fourier theory 87( §4.2.4) ensures that any function f may be approximated by the trigonometric polynomial P N f in L 2 , that is, the condition The rate of convergence may be quantified provided the function f satisfies additional regularity assumptions. More precisely, the estimate 81(Thm. 8.2.1) for any k ≥ 1 with a constant depending only on the cell Y d , the dimension d and on the order k, provided the right-hand side is finite. Although the trigonometric projection P N f of a given function f ∈ L 2 (Y d ) is rather natural in the L 2 -scale of functions, its use with other norms is rather limited. For instance, the projection P N f does not converge to the function f in the L 1 -norm, in general. Also, the truncations P N f are not bounded uniformly in N for an essentially bounded function f . For this reason it is convenient to consider the Féjer approximation see Weisz 88( §12) . The formula (21) may be rewritten in real space in the form of a periodic convolution with the non-negative and integrable kernel Then, it is immediate to see that F N f converges to f in any For the manuscript at hand, we will need two properties of the Féjer means. For a start, the convolution representation (22) and Young's convolution inequality 89 imply the estimate As the integral of the Féjer kernel K N is bounded independently of N, for example, readily observed from its Fourier representation, we thus obtain that the Féjer mean F N f of a bounded function f ∈ L ∞ (Y d ) is also bounded uniformly in N, that is, holds independently of N. Moreover we take a look at the approximation quality of the Féjer means. As mentioned previously, the trigonometric projection P N f represents the L 2 -best approximation to a given function f . Fortunately, the additional error introduced by using Féjer means is not too large. More precisely, the estimate This can be seen as follows. By the triangle inequality, and the estimate (20) for k = 1, it suffices to show By Parseval's identity (10) and the definition (21) of the Féjer means, we observe The term in squares is given as a sum of products of the form | |∕M | | for multi-indices ∈ N d whose components are not all equal to zero. Using the estimate | i | ≤ M implied by the definition of the set Z d N , see Equation (9), we are led to the bound with a constant depending only on the dimension d, which readily implies the statement (28), whose validity implies the convergence rate assertion (26). We will conclude this section with two further results about trigonometric polynomials. The first concerns convergence rates of trigonometric interpolation. For a given bounded function f , denote by Trigonometric interpolation is more subtle than trigonometric approximation as the former is not well-defined on Lebesgue spaces. Indeed, for interpolation to be reasonable, the function values to be interpolated need to be finite. Moreover, Fourier series need not converge pointwise, and arguments based on Fourier series may not be applicable. These problems may be overcome provided the function to be interpolated is sufficiently smooth. Indeed, Hölder continuity is a sufficient condition for a Fourier series to converge pointwise, and Hölder continuity is in turn implied by Sobolev embedding provided the sufficiently high derivatives of the function under consideration belong to suitable Lebesgue spaces. Then, the following result may be shown, see, for example, Vondřejc et al. 71(Lemma 4) , where the index k satisfies the condition k > d∕2 required for the Sobolev embedding theorem 90(Ch.V, Thm. 2(iii)) to hold. Thus, the convergence rate of trigonometric interpolation (31) coincides with the rate of trigonometric approximation (20), at least for sufficiently regular functions. Unfortunately, the regularity of the solutions to the micromechanical problems considered in this article is not sufficient to ensure the interpolation estimate (31) to hold. Therefore, it is convenient to use inverse estimates to trade (higher-order) derivatives for powers of the discretization parameter N. Of course, such a strategy is only feasible on a suitable finite-dimensional space of functions like the space of trigonometric polynomials. The prototypical result which we need is the following Bernstein estimate valid for any f N ∈  N (Y d ) and multi-index ∈ N d . In particular, for any trigonometric polynomial concerns one-dimensional trigonometric polynomials f N ∈  1 ([0, 2 ]), see Nikol'skii 84( §2.5) for a proof. Applying this estimate to a single coordinate j ∈ {1, 2, … , d} and noticing that the function which implies the general assertion (32) by treating a single derivative at a time.
We will also need the L 2 -version of the Bernstein estimate (32) in valid for any f N ∈  N (Y d ) and multi-index ∈ N d . The estimate (35) follows from the formulâ for the Fourier coefficients of derivatives, Parseval's identity (11) and the fact that the nonzero frequencies of a trigonometric polynomial (7) are concentrated in the set Z d N , see Equation (9), s.t., | j | ≤ N∕2 holds.

Cutting off smoothly at the interface
We are interested in composite materials, that is, those whose microstructure consists of a number of distinct materials, for example, a matrix-inclusion composite. Due to the different material behavior in the phases, the mechanical fields are not smooth across the interfaces. Rather, the solutions to the governing equations are only contained in broken Sobolev spaces.
To make this precise, we suppose that the microstructure Y d admits the non-overlapping decomposition into two open subsets Y d ± whose boundary Y d ± , considered as a Y d -periodic object, is sufficiently smooth. More precisely, the interface should be in the Hölder class C 1, for positive for the Li-Nirenberg estimates 91 to hold.
For the subsequent discussion, let denote the periodic distance to the interface Y d ± . The distance function (38) is continuous. Moreover, due to the smoothness of the interface, the distance function is also smooth in the vicinity of the interface, that is, there is a natural number N 0 , s.t. the function is sufficiently smooth at all points x ∈ Y d which are 1∕(2N 0 )-close to the interface, that is, which satisfy (x) ≤ 1∕(2N 0 ). We will tacitly assume the condition N ≥ N 0 to hold in the subsequent manuscript.
For the following construction we fix a smooth, that is, infinitely often continuously differentiable, function ∶ R ≥0 → [0, 1], s.t. the conditions hold. Such a function is easily constructed, for example, by mollifying a piece-wise linear function. We fix once and for all such a function , satisfying the conditions (39). Properties of the function, in particular upper bounds on the partial derivatives, will enter our estimates as universal constants.
With the function at hand, we define the cutoff function where N is a positive integer. For sufficiently large parameter N, the function N is smooth and permits to patch together a function f ∶ Y d → R whose restrictions f ± ∶ Y d ± → R are smooth to a smooth function on the entire cell. More precisely, the function N f coincides with the original function f for points x whose distance to the interface exceeds 2∕N. 1∕N-close to the interface, the function N f vanishes identically.
In the remainder of the section, we quantify the error introduced by this smoothing procedure in Sobolev norms. By construction (39) of the function , the function N satisfies the bound Moreover, the function N satisfies the condition and the estimate with a constant that depends on the interface, the dimension and the function . In particular, for any multi-index ∈ N d ⧵ {0} and every exponent p ∈ [1, ∞), we obtain the estimate which implies the bound Notable special cases include We will also need the elementary estimate

The trapezoidal rule for rough functions
The Moulinec-Suquet discretization 4,5 was identified as a Fourier-Galerkin discretization with numerical integration by Vondřejc et al. 71 More precisely, the trapezoidal rule (apparently in use since more than 2000 years 92 ) was used as the quadrature rule in question. The trapezoidal rule is rather simple, but turns out to be quite powerful for smooth periodic functions. 93 In one spatial dimension, the classical analysis of the trapezoidal rule shows that using N integration points leads to an error which decays as N −2 with a constant that involves the maximum of the second derivative of the function to be integrated 94(eq.(5.1.7)) . For the manuscript at hand, the L ∞ -norm is not suitable, see Section 2.2, and estimates in Lebesgue norms are required. Our subsequent discussion is based on the following one-dimensional result, which may be used to analyze the convergence behavior of the trapezoidal rule of rather rough functions, see Cruz-Uribe and Neugebauer. 95 For any positive integer N ∈ N there is a function w N ∶ [0, 2 ] → R, bounded uniformly in N, s.t., for any continuous and periodic function f ∶ [0, 2 ] → R with integrable weak derivative f ′ the representation holds, where − ∫ refers to the mean integral. The result in well known. For the reader's convenience, a derivation is included in Appendix A. In one dimension, equation (48) can be used to bound the error decay for the trapezoidal rule with a constant independent of the function f and the parameter N.
Applying the representation formula (48) for each dimension individually, we arrive at an estimate for the error of the trapezoidal rule in d dimensions.
For any dimension d ∈ N, there is a constant C, s.t. for any continuous periodic function f ∶ [0, 2 ] d → R whose derivatives are integrable up to order d, the estimate Indeed, applying the identity (48) for each dimension individually we arrive at the representation see Verlinden et al. 96( §5) , with the short-hand notation Due to the boundedness of the weight function w N , uniformly in N, the estimate (50) follows from a repeated application of the triangle inequality and Hölder's inequality in view of the N-independent bound on the functions w N . Apart from the gradient of the function f , the estimate (50) involves mixed derivatives up to the order of the dimension, for example, the quantity 1 2 3 f in dimension d = 3.

The homogenization problem
Suppose that the microstructure Y d is decomposed according to equation (37) with a sufficiently smooth interface s.t. the constructions of Section 2.2 are valid. Let stiffness tensors C ± be given, that is, linear operators on the space Sym(d) of symmetric d × d-tensors, which are symmetric in the sense that with the double-contraction operator defined by 1 ∶ 2 = tr( 1 2 ), and positive definite. Denote by and positive constants, s.t. the estimates hold for both C + and C − in terms of the norm || || = ( ∶ ) 1∕2 . Define the heterogeneous stiffness field C ∶ Y d → L(Sym(d)) via For prescribed macroscopic strain ∈ Sym(d), we seek the periodic displacement field u ∈ H 1 # (Y d ; R d ), the closure of smooth periodic and mean-free fields on the periodic cell Y d w.r.t. the H 1 -norm, solving the equation where ∇ s refers to the symmetrized gradient and the operator ⟨⋅⟩ Y is a short-hand notation for the mean integral over the unit cell Y d . Based on the inequalities of Korn and Poincaré, the Riesz representation theorem is readily applied to show existence and uniqueness of solutions to the problem (57).
Li and Nirenberg 91 showed that the solution u to the problem (57) has a uniformly bounded strain, that is, the estimate and the restriction of the displacement fluctuation field u to the sets Y d ± lies in the Sobolev space H 2 , that is, the inequalities on the decomposing domains Y d ± . The two higher regularity estimates (58) and (59) will be used to derive convergence estimates for the Moulinec-Suquet discretization. We are not aware of similar results beyond the linear (thermo-)elastic setting, which is why we restrict to this setting. Indeed, a-priori estimates and convergence results (without explicit rates) are available for the more general setting of strongly monotone and Lipschitz continuous stress operators. 57 Higher regularity results are lacking, compare for instance Neff and Knees 97 for work in this direction (treating boundary regularity instead of interface regularity).
Let us emphasize the dependence of the solution u on the macroscopic strain by a subindex, that is, we write u . Then, the effective stiffness C eff ∈ L(Sym(d)) of the microscopic stiffness distribution C is defined by Showing that this construction is well-defined and gives rise to a symmetric stiffness respecting the bounds (55) is standard, see, for example, the arguments in Schneider 72(Corollary 2.2) . For our purposes, we need the representation formula which is a direct consequence of the validity of the cell problem (57) and the definition (60).

The a-priori estimate
For any odd positive integer N, the Moulinec-Suquet discretization 4,5 seeks a displacement field is satisfied for all test fields v N ∈  N (Y d ; R d ). As observed by Vondřjec et al., 71 the Moulinec-Suquet discretization proceeds by restricting the space of ansatz and test fields to trigonometric polynomials and replaces the integrals appearing in the cell problem (57) by the trapezoidal rule. For the analysis, it is convenient to recast the Equation (62) in terms of the trigonometric interpolation operator Q N . 71,72 Indeed, due to Parseval's identity (11), Equation (62) may be rewritten in the form

Introducing the stress field
, the equation (63) is, in turn, equivalent to the strong form of the quasi-static balance of linear momentum. This property contrasts with traditional finite element discretizations, where both the compatibility and the material law are evaluated exactly, but the equilibrium equation (64) is approximated. Indeed, the Moulinec-Suquet discretization works with a compatible strain field N = + ∇ s u N ∈  N (Y d ; Sym(d)) and a self-equilibrated stress field N ∈  N (Y d ; Sym(d)), but the constitutive law is valid only in the collocation (or quadrature) points. Under the hypothesis (55), the equation (62) admits a unique solution u N among the trigonometric polynomials  N (Y d ; R d ) with vanishing mean. This can be seen as follows. 71 Writing the problem (62) in the suggestive form existence and uniqueness follows from the Riesz representation theorem 87(Thm.4.17) . Indeed, the right-hand side is apparently linear in the test field, whereas the left-hand side represents a symmetric bilinear form in the fields u N and v N , respectively. Due to the bound (55), we may estimate where we used Parseval's identity, Korn's inequality 98 and Poincaré's inequality 99( §5.8.1) on the flat torus. After recalling the argument 71 for existence and uniqueness of solutions to the discretized cell problem (62), we provide a stream-lined argument for a-priori estimate 71(Proof of Prop.8) which is essential for the analysis in Section 3.3. Denote by N = + ∇ s u N the strain associated to the solution u N of the cell problem (62) in response to the imposed strain ∈ Sym(d). We denote bỹN any compatible competition strain. Then, the difference between these two strains admits the representation In particular, we observe in view of the discretized cell problem (63). Arguing in a reverse order but using the continuous cell problem (57) yields the chain of arguments Applying the Cauchy-Schwarz inequality yields the desired a-priori estimate This inequality served as the basis of previous convergence proofs 71-73 which either required higher regularity of the solution field u or led to convergence results without explicit rates. The mentioned strategies selected N = P N with = + ∇ s u, that is,̃N = + ∇ s P N u (73) in terms of the trigonometric approximation operator (18).

Convergence rate for the strain field
The goal of this section is to show the convergence-rate estimate for the strain field in the Moulinec-Suquet discretization, where the constant depends on the dimension, the interface, the ellipticity ratio, the L ∞ -bound on the strain field and the H 1 -bounds on the strain field on the subdomains Y d ± . To analyze the convergence behavior of the strain field, our strategy utilizes Féjer means (21) instead of the trigonometric projection (18). More precisely, instead of equation (73), we set By the triangle inequality and the a-priori estimate (72), we deduce Thus, the discretization error is naturally split into an approximation error and an interpolation error Thus, to derive convergence rates for the discretization error, we analyze these two contributions individually. For the approximation error (77), we would like to apply classical convergence rate estimates for trigonometric approximation (20). However, the strain field lacks the necessary regularity, and we rely upon the smoothing approach introduced in Section 2.2. We write where we used that the Féjer-mean operator gives rise to a bounded linear operator on L 2 with operator norm 1. For the first term, we have, by the estimate (47) for p = 2, whereas, for the second term, we observe where we used the convergence-rate estimate for Féjer approximation (26) and the growth estimate (46) for the gradient of the cut-off function N . Thus, the approximation error (77) decays as N −1∕2 .
Let us turn our attention to the interpolation error (78). We would like to apply the error estimate (31) of trigonometric interpolation. However, the integrand is not smooth. As for the approximation error, we use the cut-off function N to split the interpolation error (78) into the contributions For the first term, we observe where we used the decay estimate (47) for p = 2 and the crucial property (25) of the Féjer means. Notice that this would not hold for the Fourier projection P N . Using a similar reasoning in view of Parseval's identity (11) permits to bound the third term in the estimate (82) To treat the middle term in the estimate (82), we invoke the error estimate (31) for trigonometric interpolation We notice that higher-order derivatives appear on the right-hand side, that is, second derivatives in the physically relevant dimensions d = 2 and d = 3. However, we may only work with a uniform bound (58) on the strain and the H 1 -bounds (59) for the phase strains. To deal with this situation, we utilize the Bernstein estimates (32) and (35). Then, the estimate can be shown, see Appendix B for details, which in turn implies Together with the estimates (83) and (84), we are thus led to the following bound on the interpolation error (82) In combination with the bound (81) on the approximation error, we are, in view of the decomposition (76), finally led to the estimate (74).

Convergence rate for the stress field
The goal of this section is to show the convergence-rate estimate for the stress field when approximating the continuous stress field = C ∶ . Due to the special nature of the Moulinec-Suquet discretization, the convergence estimate (89) does not immediately follow from the established convergence rate (74) for the strain field, as the discretized stress field (90) involves an additional trigonometric interpolation step (13). However, to establish the estimate (89) we may heavily re-use the machinery developed in the previous section. Indeed, recasting the stress error in the form where we used Parseval's identity (11) and the bound (55) on the stiffness, we recognize the approximation error (77), the interpolation error (78) and the total error (76). For each of these terms, the convergence rate N −1∕2 was established in the previous section. This shows the validity of the statement (89).

Convergence rate for the elastic energy
In this section, we wish to establish a convergence estimate of the form that is, the effective stiffness of the Moulinec-Suquet discretization converges to the continuous effective stiffness (60) with twice the rate of the strains (74). In the Galerkin setting, the estimate (92) follows from the estimate for the strain (74) by Galerkin orthogonality, see Ye and Chung 73(Thm. 5) or Schneider and Wicht 76( §2.2) . The Moulinec-Suquet discretization 4,5 requires additional arguments to treat the underintegration. To proceed, we fix a macroscopic strain ∈ Sym(d) and recall the representation (61) solves the cell problem (57). For any competition field we expand the square where we used that the difference −̃= ∇ s (u −ũ) is a symmetrized gradient and that the field satisfies the equilibrium equation (57). We may rewrite this identity in the form Completely analogous arguments may be applied to the discretized setting for a competition field based on the discrete equivalent of Equation (61) to yield Subtracting the identity (97) with̃=̃N from the previous line (100) provides the exact expression Taking norms thus leads us to the estimate where we used Parseval's identity (11) and the upper bound (55) on the stiffness. Taking Féjer means the results of the previous section, that is, the estimate (74), permit us to conclude that the estimate holds. Therefore, in view of the goal (92), it is necessary to study the quadrature error in the estimate (102), where we used the fact (take Equation (11) with g N ≡ 1) that trigonometric polynomials of degree N are integrated exactly by the trapezoidal rule of order N.
To proceed, let us introduce the field which is a trigonometric polynomial of order 2N with values in fourth-order tensors, and rewrite the quadrature error (105) in the slightly shorter form We wish to apply the quadrature estimate (50). However, the integrand C ∶∶Ẽ N is discontinuous due to the discontinuous stiffness C. Therefore, we invoke the splitting in terms of the cut-off function N , introduced in Equation (40). For the first term, we observe where we used Hölder's inequality, the bound on the stiffness (55), the definition (106) of the fieldẼ N , the property (25) of Féjer means and the inequality (47) for p = 1. With a completely analogous argument, the estimate is readily seen to hold. To treat the quadrature error in the middle of the right hand side of equation (108), we make use of the quadrature estimate (50) To get an idea how to handle the terms on the right-hand side, we compute the first (j = 1, 2, … , d) and the second derivative (j, k = 1, 2, … , d) explicitly. From the expression (112), we estimate where the critical ingredient is the L 1 -bound (46). The first term on the right hand side of Equation (113) may be handled as follows by the estimate (45) for p = 1 and | | = 2. To treat the second term in the identity (113), we notice where, apart from the uniform L 1 -bound (46) the Bernstein estimate (32) is crucial to remove the derivative from the L ∞ -estimate for the trigonometric polynomial̃N. The other terms in Equation (113) can be estimated with similar ideas to obtain the bound Actually, completely analogous arguments may be used to derive the estimate for higher derivatives, as well. Thus, in view of the bound (111), we are led to the estimate Together with the estimates (109) and (110), the bound on the quadrature error (105) emerges. Combined with the quadratic estimates (104), the derivation of the desired statement (92) is complete.

COMPUTATIONAL INVESTIGATIONS
The goal of this section is to provide computational results showing that the obtained convergence rates for the local strain and stress fields as well as for the effective stiffness are actually attained upon grid refinement. Indeed, our considerations could be suboptimal, requiring more elaborate mathematical techniques to derive optimal rates. We relied upon an in-house FFT-based computational homogenization code, 41 run on a workstation with two AMD EPYC 7642 with 48 physical cores each. To handle discretizations with an even number of voxels in certain coordinate directions, the values of the strain field are set to zero at the Nyquist frequencies. For every considered mesh size, the equilibrium equation (64) was solved to an accuracy of 10 −5 w.r.t. the "natural" convergence criterion 14( §3.6) . Then, the effective stresses were computed by volume averaging.
We consider a microstructure with a single spherical inclusion placed at the center of the microstructure, see Figure 1A, with a volume fraction of 6.54%. The isotropic linear elastic material parameters are listed in Table 1, corresponding to a polyamide matrix and an E-glass filler. The considered microstructure satisfies the conditions presented in Section 2.2, in particular the smoothness of the interface, required for the Li-Nirenberg estimates to hold, see Section 3.1. Thus, the theory developed in this article applies. We load the structure via a shear in the x-y-plane.
The highest considered resolution, 512 3 voxels, is taken as the reference. Then, the relative errors of the strain and stress field, both measured in L 2 (Y d ), and the effective stress are shown in Figure 1B. We observe that the local solution fields, that is, the strain and the stress fields, show an N −1∕2 convergence rate, as predicted by the theory, see estimates  (74) and (89). The effective stress, on the other hand, converges much faster, and shows roughly an N −1 -rate, as predicted by our derivations, see the estimate (92). The previous example was admittedly simple in order to enable going to fine resolutions and satisfying all the requirements of the theory.
We retain the material models and parameters of the previous example, but investigate a more complex short-fiber reinforced microstructure. More precisely, we take a look at a material reinforced by 20% E-glass fibers with an aspect ratio of ten and a second-order fiber-orientation tensor 101,102 which describes a small deviation from a planar isotropic fiber-orientation state. As the fibers are described by cylinders, the interface between matrix and inclusions does not satisfy the C 1, -condition required for the Li-Nirenberg estimates, 91 and the theory developed in this article does not apply. The microstructure has dimensions 256 × 256 × 64 m 3 , shown in Figure 2A and contains 106 fibers. It was generated by the SAM 103 algorithm, employing the exact closure approximation. 104 We consider resolutions N × N × N∕4 with discretization parameters N in the set {64,128, 256, … , 1024}. A non-cubical microstructure is used to permit considering a larger resolution in the principal fiber directions. Indeed, due to the anisotropic reinforcements, higher stresses are expected in-plane compared to out-of-plane. We apply a shear strain in the 1-2-plane, that is, in the fiber plane. The computation with N = 1024 serves as our reference.
The relative errors of the stress and the strain field (in L 2 ) as well as the convergence of the effective stress is shown in Figure 2B. We observe that the local fields converge as N −1∕2 in the L 2 -norm. Moreover, the errors in the effective stresses decay as 1∕N. In particular, the convergence asymptotics match the simpler spherical case which we considered previously.
We conclude this section by investigating a microstructure where an analytic solution for the local solution fields is available, at least for a specific macroscopic loading. More precisely, we consider the neutral inclusion introduced by Hashin. 105 The microstructure consists of three phases, a spherical core an annular region and the remaining "matrix" with the geometric center x c = ( , , ) of the microstructure and two radii r 1 < r 2 < 2 . Each of these three phases is occupied by an isotropic elastic medium. Hashin 105 showed that under macroscopic compression it is possible to select the elastic moduli of the phases Ω 1 and Ω 2 in such a way that the effective compression modulus of the microstructure coincides with the compression modulus of the matrix phase. The shear moduli 1 and 3 of the core and the matrix are actually irrelevant for the consideration. Typically, we fix the compression modulus K 3 ≡ K eff of the matrix and the ratio K 2 ∕ 2 of the intermediate region. There is a relationship 105 between the compression moduli K 1 and K 2 , that is, where 1 = (r 1 ∕r 2 ) 3 and 2 = 1 − 1 . For our investigation at hand, there are two cases to consider. For the first case, which we term "porous", the compression modulus K 1 is less than K 3 ≡ K eff . It follows that the compression modulus K 2 of the intermediate region exceeds K 3 ≡ K eff . With decreasing K 1 , K 2 increases. However, as K 1 goes the zero, K 2 remains bounded. Thus, we approach the case of a single pore with a special coating embedded in the matrix.
In the second case, which we call "rigid", the compression modulus K 1 exceeds K 3 ≡ K eff , whereas K 2 is smaller than K 3 ≡ K eff . For increasing K 1 , K 2 decreases. As K 1 → ∞, K 2 remains bounded away from zero. Thus, in the limit K 1 → ∞, we are faced with a rigid inclusion with special coating inside the matrix.
These two scenarios permit us to consider different material contrasts in the compression modulus.
For the study at hand, We use the radii where e = ∑ ∞ k=0 1 k! is Euler's number. We set K eff = 1∕3 and fix i = K i (i = 1, 2, 4), that is, we set Poisson's ratio to i = 1∕8 (i = 1, 2, 3). We compare the relative error in the effective compression modulus and the L 2 -error of the local strain and stress fields, see Figure 3. Inspecting the error in the effective properties for the "porous" case, see Figure 3A, we observe that the two coarsest discretizations, N = 16 and N = 32, do not conform to the expected scaling 1∕N. This contrasts with the geometries previously considered in this section, and is caused by the comparatively coarse resolution of the regions Ω 1 and Ω 2 . For the finer discretizations the theoretically predicted 1∕N-convergence rate is confirmed. We observe that the error increases with the contrast, with no clear limit as → ∞.
For the "rigid" case, see Figure 3B, the behavior is similar, and more or less the same conclusions may be drawn. However, the errors are slightly larger than for the "porous" case. Please note that resolving the latter is computationally much more expensive than dealing with the "porous" case due to the primal formulation. A remedy consists in working with the dual formulation. 62,106 Taking a look at the local L 2 -strain error of the "porous" case, see Figure 3C, we also observe no clear trend in the convergence behavior for low resolution due to the geometric error induced by the coarse voxelation. Only for the highest resolutions, the theoretically predicted N −1∕2 convergence behavior is clearly seen. We observe that the strain error increases with the contrast, even strongly so. Keeping in mind that the axes are scaled logarithmically, a divergence of the strain error is evident. This behavior is a well-known and understood characteristic of the Moulinec-Suquet discretization. 47,50,57 Taking a look at the stress error for the "porous" case, see Figure 3E, we observe a completely different behavior. As the contrast increases to infinity, there is a clear saturation behavior of the stress error, following the predicted N −1∕2 -rate.
For the "rigid" case, see Figure 3E,F, the roles of stress and strain are reversed-as expected-but the conclusions to be drawn remain the same.

CONCLUSION
The goal of this article was to provide a sharp error analysis for the Moulinec-Suquet discretization, 4,5 providing the theoretical basis for its practical success in computational micromechanics. Previous works were mostly restricted to qualitative convergence results. 66,72,73 The only quantitative results that we are aware of either require higher smoothness of the coefficients, 71 are restricted to specific microstructures, that is, laminates, 107 or deal with the Galerkin setup. 108 The principal difficulty involved was the low regularity characterizing the solution of the micromechanical problem and composite materials. Indeed, both the stress and the strain fields are, in general, discontinuous across material interfaces, a situation that is inherently linked to the physics of the problem.
We employed a strategy based on smoothing and elementary energy estimates to deal with the problem, losing half of the optimal 1∕N-rate in the process due to balancing the smoothing and the approximation error. Moreover, we could handle the higher derivatives appearing in the interpolation estimates by using Bernstein estimates for the gradient of the trigonometric polynomial. Usually it is better to work out estimates in terms of the non-discretized, that is, continuous, solution field. In our case, the opposite strategy turned out to be vital. Another key ingredient of the approach is the use of Féjer means to be able to preserve the L ∞ -bound upon trigonometric approximation. It would be interesting to see whether this theoretical tool is also useful for practical matters.
In addition to providing rates for the strain and the stress fields in L 2 , we also investigated the error of the effective stiffness (or the effective stresses). For this purpose, we took a close look at the error induced by trigonometric interpolation. On top of the Ye-Chung superconvergence argument, 73 a similar reasoning as for the strain/stress estimates was necessary, replacing the L 2 -by L 1 -estimates.
Let us take a critical look at the assumptions and possibilities for future work. This work was restricted to linear elasticity and smooth material interfaces as we relied critically on the Li-Nirenberg estimates 91 which were formulated in this setting. Suitable a-priori estimates 72 and a Ye-Chung type argument 76 are available for more general stress operators of strongly monotone type. Thus, generalizations of the Li-Nirenberg estimates to this kind of nonlinearity would probably enable extending the work at hand to this more general case, as well. Actually, computational experiments 76 suggest that the reported rates also hold for relevant classes of nonlinear composites and interesting classes of non-smooth material interfaces. To treat the latter analytically, finer tools are required, for example, based on weighted function spaces. 109,110 Another intriguing question concerns whether it is possible to design computational methods on Cartesian grids which converge at a higher rate than the Moulinec-Suquet discretization (and its relatives), but comes with a similar computational efficiency, for example, by using the fast Fourier transform.

Support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)-255730231-and the European
Research Council within the Horizon Europe program-project 101040238-is gratefully acknowledged. The author is grateful to D. Knees for shedding light on some results of regularity theory of elliptic systems, L. Risthaus for constructive feedback on an initial version of the manuscript and the anonymous reviewers for carefully reading the submitted manuscript. Open Access funding enabled and organized by Projekt DEAL.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.

APPENDIX A. DERIVATION OF THE REPRESENTATION OF THE QUADRATURE ERROR IN ONE SPATIAL DIMENSION
The purpose of this appendix is to show the validity of equation (48), that is, for any natural number N ∈ N there is a function w N ∶ [0, 2 ] → R, bounded uniformly in N, s.t., for any continuous and periodic function f ∶ [0, 2 ] → R with integrable weak derivative f ′ the representation holds, where − ∫ refers to the mean integral. For 0 ≤ a < b ≤ 2 one observes, by the fundamental theorem of calculus, On the other hand, the product rule gives that is, the compact formula For a = 2 j∕N and b = 2 (j + 1)∕N (j = 0, … , N − 1), this formula becomes 2 2N ( f Summing up from j = 0 to j = N − 1 gives equation (48) with the piece-wise defined function w N (x) = Nx − (2j + 1), 2 j∕N ≤ x ≤ 2 (j + 1)∕N.

APPENDIX B. DERIVATION OF THE HIGHER GRADIENT BOUNDS IN THE INTERPOLATION ESTIMATE
The goal of this appendix is to show the estimate (86) required in Section 3.3. The product rule implies For the first term, we observe by the L 2 -Bernstein estimate (35). We split the term For the first contribution, we utilize the Bernstein estimate (35), the boundedness of the Féjer means and the convergence-rate estimate (46) (B5) To treat the second term, we notice that taking the gradient commutes with Féjer approximation, where the last estimate is obtained as a part of argument (81). Thus, we are let to the bound which, in turn, leads to the expression For the second term in the estimate (B2), we observe as a consequence of estimate (45) for k = | | as well as p = 2, and using the Bernstein estimate (32) and the conservation of boundedness offered by the Féjer means (25). The latter two estimates combine to the bound Thus, in view of the contributions (B8) and (B11), the estimate (B2) implies the statement (B1), which was to be shown.