Grigor Nika
email: grigor.nika@kau.se

Adrian Muntean
email: adrian.muntean@kau.se

Hypertemperature effects in heterogeneous media and thermal flux at small-length scales

Keywords: MSC 2020 74Q05, 35B27, 35G05, 35G15, 35J58, 35Q79 correctors, scale-size thermal effects, generalized Fourier's law, microstructure

We employ an enriched microscopic heat conduction model that can account for size effects in heterogeneous media. Through, physically, relevant scaling arguments we improve the regularity of the corrector in the classical problem of periodic homogenization in the three-dimensional setting and, in doing so, we clarify the intimate role correctors play in measuring the difference between the heterogeneous solution (microscopic) and the homogenized solution (macroscopic). Moreover, if the data are of the form f = div F F F with F F F ∈ L 3 (Ω, R 3 ), then we can prove the classical corrector convergence theorem as well.

Introduction

Analysis of correctors in homogenization theory for second-order elliptic equations with highly oscillatory coefficients provides an important link between the microscopic and macroscopic problem since it measures the difference between the heterogeneous solution and the homogenized solution. Simultaneously, correctors point out some of the limitations of such second-order elliptic systems due to their lack in accounting for scale-size effects in e.g. microheterogeneous bodies. The prototypical and landmark example for the homogenization of second-order elliptic equations is Fourier's law of heat conduction. In its simplest form, Fourier's law relates the heat flux q q q as a linear function of the temperature gradient, that is:

q q q = -κ∇u, (1.1)
where u is the absolute temperature and κ > 0 the thermal conductivity that depends on properties of the material. In general, the coefficient κ may depend on temperature, space, and/or time, but often varies so little in cases of interest that it is reasonable to neglect this variation.

The theory of periodic homogenization led to a deeper exploration of Fourier's heat conduction law as it pertains to heterogeneous periodic material with different conductivities.

Effective heat fluxes were derived taking into account microstructure morphology and volume fraction. The classical problem in the periodic homogenization for stationary heat conduction states:

-div K( x x x ε )∇u ε = f in Ω, u ε = 0 on ∂Ω, (1.2)
where K(y y y) ∈ L ∞ (Y, R 3×3 ) is uniformly elliptic, symmetric, and Y -periodic with Y = [0, 1) 3 . If Ω ⊂ R 3 is uniformly Lipschitz open set, then there exists a unique solution u ε to (1.2) such that it converges weakly to a function u in H 1 0 (Ω), where the function u ∈ H 1 0 (Ω)

is the unique solution to,

-div K eff ∇u = f in Ω, u = 0 on ∂Ω, (1.3)
with K eff ij := 3 l=1 Y K il (y y y)(δ jl -∂ y l w j ) dy y y and w j ∈ {H 1 per (Y ) | Y w j dy y y = 0} solution to, -div y (K(y y y)(∇ y w j -e e e j )) = 0 in Y, w j is Y -periodic.

(1.4)

The convergence of u ε to u in H 1 0 (Ω) is only weak. If one seeks to improve the convergence then, usually, a corrector type term is introduced,

u ε -u -εû(•, • ε ) → 0 in H 1 (Ω), (1.5)
where û(x x x, y y y):= -N k=1 w k (y y y) ∂u ∂x k (x x x) is the correction term. In order for the expression to belong in H 1 (Ω) it is required that w k ∈ W 1,∞ (Y ) and u ∈ H 2 (Ω) (see [START_REF] Bensoussan | Asymptotic Analysis for Periodic Structures[END_REF]Sect. 5,pg. 33]). With the introduction of two-scale convergence [START_REF] Nguetseng | A general convergence result for a functional related to the theory of homogenization[END_REF] (see also [START_REF] Allaire | Homogenization and two-scale convergence[END_REF] and [START_REF] Lukkassen | Two-scale convergence[END_REF]) a rigorous justification was provided for the multiple-scale method and the corrector result (in N -dimensional space) in (1.5) was made rigorous through the following theorem:

Proposition 1.1. [1, Thm 2.6], [8, Prop. 9.12], [START_REF] Lukkassen | Two-scale convergence[END_REF]Thm. 22] Let û be given by,

û(x x x, y y y) = - N k=1 w k (y y y) ∂u ∂x k (x x x) + ũ1 (x x x), (1.6)
and suppose that ∇ y w k ∈ L r (Y, R N ), k = 1, . . . , N and ∇u ∈ L s (Ω, R N ) with 1 < r, s < ∞ and such that,

1 r + 1 s = 1 2 . (1.7)
Then,

∇u ε -∇u -∇ y û(•, • ε ) → 0 in L 2 (Ω, R N ). (1.8)
As one can immediately observe, the main unease with the above proposition is the higher integrability required on both local and homogenized solutions. However, with the introduction of periodic unfolding operators [START_REF] Ciorȃnescu | Éclatement périodique et homogénéisation[END_REF][START_REF] Ciorȃnescu | The periodic unfolding method in homogenization[END_REF], one can obtain a more general corrector result without requiring any regularity assumption on the cell function w k , stating

∇u ε -∇u - N k=1 Q ε ( ∂u ∂x k )∇ y w k → 0 in L 2 (Ω, R N ), (1.9)
where Q ε is the scale-splitting operator defined in [START_REF] Ciorȃnescu | Éclatement périodique et homogénéisation[END_REF]; see the elegant proof in [START_REF] Ciorȃnescu | The Periodic Unfolding Method. Theory and Applications to Partial Differential Problems[END_REF]. Moreover, based on the scale-splitting operator, upper bound estimates on the convergence rate in terms of ε were obtained in [START_REF] Griso | Error estimates and unfolding for periodic homogenization[END_REF] (with some additional regularity assumptions on the homogenized solution). Furthermore, the upper bound estimates on the convergence rate can be made tighter by using boundary layer correctors (see e.g. [START_REF] Moskow | First-order corrections to the homogenised eigenvalues of a periodic composite medium. a convergence proof[END_REF], [START_REF] Allaire | Boundary layer tails in periodic homogenization[END_REF], [START_REF] Onofrei | Error estimates for periodic homogenization with non-smooth coefficients[END_REF]). Therefore, it seems that in order to lift the (restrictive) regularity conditions in Proposition 1.1 one must have knowledge of the operator Q ε .

In this work, we commence with a higher-gradient heat equation model so as when passing to the limit the homogenized problem coincides with (1.3), however, the local solutions w k satisfy a higher-gradient local problem and, hence, possess better regularity properties than classical local solutions. This aforementioned higher regularity of the local solution compensated by the mild assumption that the data are of the form

f := div F F F , F F F ∈ L 3 (Ω, R 3 )
allows us to prove Proposition 1.1 with weaker than usual assumptions and without using the scale-splitting operator. The approach marshalled in this article is motivated from generalized continuum mechanic theories (see e.g. [START_REF] Forest | Milieux continus généralisés et matériaux hétérogènes[END_REF]). It is well understood that the effective properties of heterogeneous materials can depend not only on the volume fraction of the phases or their geometrical distribution but also on the absolute size of the constituents (see e.g. [START_REF] Forest | Asymptotic analysis of heterogeneous Cosserat media[END_REF]). By introducing physically relevant scaling arguments related to the absolute size of the constituents, we introduce a new length scale parameter that models scale-dependent thermal effects (see e.g. [START_REF] Fehér | Size effects and beyond-Fourier heat conduction in room-temperature experiments[END_REF][START_REF] Fehér | On the evaluation of non-Fourier effects in heat pulse experiments[END_REF]). Hence, by arguing as in the work of [START_REF] Nika | Effective medium theory for second-gradient nonlinear elasticity with chirality[END_REF], the solution of the enriched microscopic problem can be seen as a vanishingviscosity solution that coincides with the classical homogenized solution of (1.3). The higher regularity that the enriched microscopic model possesses is transferred to the local solution with the minimal assumptions of a uniformly Lipschitz open set Ω and non-smooth coefficients.

We have organized the paper as follows: In Section 2 we explain in detail the scaling argument we employ, present the enriched microscopic model, provide some motivation for its use, and prove some general qualitative results as they pertain to existence and uniqueness of solution as well as the variational nature of the problem. In section Section 3

we state the main results, we discuss their consequences, and demonstrate symmetry relations for the higher-gradient effective coefficients as well as explore their variational structure. Section 4 is dedicated to proving the main results in Section 3. Finally, we reserve Section 5 for some discussion and remarks.

2 Problem set-up

Generalized Fourier's law

Fourier's law of heat conduction is widely regarded as a limiting approximation of some general (potentially nonlinear) constitutive law for the heat flux that may depend on higher-gradients [START_REF] Coleman | Thermodynamics and departures from Fourier's law of heat conduction[END_REF] [7], [START_REF] Pan | A generalized heat conduction model in rarefied gas[END_REF]. For instance, the temperature of a rarefied gas at the slip regime, namely when 0.001 < Kn < 0.3 where Kn is the non-dimensional Knudsen number, deviates from Fourier's law of heat conduction [START_REF] Pan | A generalized heat conduction model in rarefied gas[END_REF]. Moreover, in the same article a generalized heat conduction model, from a phenomenological point of view, was postulated under the assumption that the gas is isotropic. The authors' speculation was that the heat flux in a rarefied gas in the slip regime, depends linearly on on the temperature gradient but also on higher-order temperature derivatives,

q q q = K∇u + L:∇∇u + M . . . ∇∇∇u, (2.1)
where K is a second order tensor, L is a third order tensor, and M is a fourth order tensor with : and . . . denoting second and third order contraction, respectively. In the case where we have isotropy and the medium is centro-symmetric, the higher-gradient heat flux can reduce further,

q q q = K∇u + M . . . ∇∇∇u, (2.2)
where K ij :=ηδ ij and M ijkl :=ℓ 1 δ ik δ jl + ℓ 2 δ il δ jk + ℓ 3 δ ij δ kl for scalars η and ℓ i , i = 1, 2, 3 that will be assumed to be constant or piece-wise constant. The heat conduction inequality asserts that, q q q.∇u ≤ 0, (

for all temperature fields. Furthermore, the tensor M has a physical meaning and is often referred to as the spatial retardation (see [START_REF] Christov | On a higher-gradient generalization of Fourier's law of heat conduction[END_REF]).

In addition to the above mentioned works, emerging primarily from the physics community, non-classical laws of Fourier's heat conduction have, for many years now, attracted considerable attention from the theoretical mechanics community as well (see e.g. [START_REF] Forest | Hypertemperature in thermoelastic solids[END_REF], [START_REF] Berezovski | Internal variables representation of generalized heat equations[END_REF]). In recent years, the motivation for deriving non-classical heat conduction models in the mechanics field stemmed from trying to understand the presence of thermal fluctuation fields in heterogeneous materials with a microstructure. Specifically, the authors in [START_REF] Forest | Hypertemperature in thermoelastic solids[END_REF] postulate the existence of a free energy function that has an added dependence on the gradient of the entropy density variable. Based on this enhanced free energy, an enhanced heat equation was derived containing a term with a characteristic length related to material parameters that can account for scale-size thermal effects in microheterogeneous bodies.

Finally, all of the above theoretical or computational non-classical approaches seem to have found some validation in recent experimental work where evidence of size-dependent thermal effects were reported in heterogeneous materials (see e.g. [START_REF] Fehér | Size effects and beyond-Fourier heat conduction in room-temperature experiments[END_REF][START_REF] Fehér | On the evaluation of non-Fourier effects in heat pulse experiments[END_REF]).

Scale-dependent thermal effects

In this paragraph, we postulate the modified heat flux in (2.2) and derive an additional length scale parameter that encapsulates the size-dependent thermal effects in the context of an idealized periodic microstructure. We assume that our working domain Ω is periodic with period ℓ and with overall characteristic length L. Define the dimensionless coordinates and temperature,

x x x * = x x x L , u * (x x x * ) = u(x x x) L . (2.4)
Moreover, we can define the following non-dimensional tensors:

KK * = K, MM * = M, (2.5)
where

K:= max z z z∈Y ℓ |K(z z z)|, M:= max z z z∈Y ℓ |M(z z z)| with Y ℓ :=(-ℓ/2, ℓ/2] 3 the periodic cell
characterizing the body Ω. We can now introduce an additional length scale relation between K and M as follows:

M = ℓ 2 TE K. (2.6)
Thus, the non-dimensional heat flux becomes,

q i = 3 j=1   KK * ij ∂u * ∂x * j + 3 k,l=1 K ℓ TE L 2 ∂ ∂x * j M * ijkl ∂ 2 u * ∂x * k ∂x * l   . (2.7)
If we use the notation, q q q * :=K -1 q q q then we have a non-dimensional form of the heat flux,

q q q * = K * ∇u * + ℓ TE L 2 div * (M * :∇ * ∇ * u * ) . (2.8)
We remark that since the coefficients K and M are Y ℓ periodic the corresponding normalized coefficients K * and M * are Y * periodic where

Y * := ℓ L Y with Y :=(-1/2, 1/2] 3 . Finally,
henceforth, if no confusion arises we will drop the * notation in order to expedite our presentation.

The microscopic problem

We consider a material with a periodic microstructure of period ε:=ℓ/L ≪ 1 occupying a region Ω ⊂ R 3 . The region Ω that the heterogeneous material occupies is assumed to be a uniformly Lipschitz open set (see [START_REF] Demengel | Functional Spaces for the Theory of Elliptic Partial Differential Equations[END_REF]Definition 2.65]). The exterior boundary component will be denoted by Σ:=∂Ω while the vector n n n will denote the unit normal on Σ pointing in the outward direction. The ε periodic problem, generated by defining the nondimensional number ε as the ratio of ℓ/L, will permit us to obtain an effective equation when ε → 0. However, unlike in classical homogenization problems, different cases ought to be considered depending on how the intrinsic length scale ℓ TE scales with ℓ (or L).

Here, since we are interested in recovering Fourier's classical law of heat conduction as an effective limit, we will only consider the scaling,

ℓ TE /ℓ ∼ 1.
(2.9)

The physical meaning of the above scaling, is that the intrinsic length ℓ TE is comparable with the length of the heterogeneities. Naturally, other type of scalings are possible, however, we will not address other cases here. We refer the reader to [START_REF] Nika | Effective medium theory for second-gradient nonlinear elasticity with chirality[END_REF] for different type of scalings in the context of generalized continuum mechanics.

Therefore, under the scaling in (2.9), the (generalized) heat flux becomes,

q q q ε = K( x x x ε )∇u ε + ε 2 div M( x x x ε ):∇∇u ε . (2.10)
The microscopic problem is then characterized by the following equation and boundary conditions,

-div K( x x x ε )∇u ε -ε 2 div M( x x x ε ):∇∇u ε = f in Ω, ε 2 M( x x x ε )∇∇u ε :n n n ⊗ n n n = 0 on Σ, u ε = 0 on Σ, (2.11)
where f is some given source that belongs in L 2 (Ω). We remark, that prescribing a homogeneous Dirichlet boundary condition, as is usually the case, is no longer sufficient.

We require, additionally, to prescribe a zero heat flux for what we refer to as a normal double heat flux that is directly related to the spatial retardation coefficient M (see next section on the weak formulation of the problem). In direct analogy with second-gradient elasticity, such a term is referred to as normal double traction (see [START_REF] Mindlin | On first strain-gradient theories in linear elasticity[END_REF][START_REF] Germain | La méthode des puissances virtuelles en mécanique des milieux continus, I: Théorie du second gradient[END_REF][START_REF] Germain | The method of virtual power in continuum mechanics. Part 2: Microstructure[END_REF]).

Notation and assumptions

-We employ the Einstein notation of repeated indices unless otherwise stated.

-Throughout the work we assume that the uniform strong ellipticity condition holds, i.e., there exist positive (generic) constants c 1 and c 2 such that:

c 1 |w w w| 2 ≤ w w w.K(x x x).w w w ≤ c 2 |w w w| 2 ,
c 1 |w w w| 2 |q q q| 2 ≤ w w w ⊗ q q q:M(x x x):w w w ⊗ q q q ≤ c 2 |w w w| 2 |q q q| 2 , (2.12) for all w w w, q q q ∈ R 3 -{0 0 0}.

Auxiliary formulas

For the readers' convenience and for the expediency of the our results, we introduce certain formulas that we will make use of in obtaining the variational formulation of (2.11). These formulas, among others, can also be found in [START_REF] Germain | La méthode des puissances virtuelles en mécanique des milieux continus, I: Théorie du second gradient[END_REF]Appendix].

For any smooth enough scalar function ξ defined on Σ or on a neighborhood of Σ, the Thus, we can write down a useful integration by parts formula on surfaces,

Σ ∇ s ξ ds = Σ ξ(div n n n)n n n ds + ∂Σ ξν ν ν dℓ, (2.15)
where

ν i = ϵ ϵ ϵ ijk t j n k , (2.16)
is the unit normal vector on ∂Σ and tangent to Σ, t t t is the unit tangent vector to ∂Σ, and ϵ ϵ ϵ is the Levi-Civita tensor. Lastly, we remark, the jump term on (2.15) is on a ridge, i.e.

the line on Σ where the tangent plane of Σ is discontinuous.

Remark 2.1. The above formulas are used with a high degree of frequency in emulsions and capillary fluids (see e.g. [START_REF] Nika | Rate of convergence for a multiscale model of dilute emulsions with non-uniform surface tension[END_REF]). We refer the reader to the appendix of reference [START_REF] Germain | La méthode des puissances virtuelles en mécanique des milieux continus, I: Théorie du second gradient[END_REF] or [START_REF] Germain | The method of virtual power in continuum mechanics. Part 2: Microstructure[END_REF] for an excellent exposition of the above formulae and related topics.

Weak formulation

The primary setting for for the variational formulation of (2.11) is the space H 2 (Ω)∩H 1 0 (Ω)

where the Sobolev space H 2 (Ω) is a Hilbert space with norm,

∥u∥ H 2 (Ω) = ∥u∥ 2 L 2 (Ω) + ∥∇u∥ 2 L 2 (Ω,R 3 ) + ∥∇∇u∥ 2 L 2 (Ω,R 3×3 ) 1/2 . ( 2
- Σ K ij ( x x x ε )∂ x j u ε -ε 2 ∂ x k M ikpq ( x x x ε )∂ 2 xpxq u ε n i v ds + Ω K ij ( x x x ε )∂ x j u ε ∂ x j v dx x x-ε 2 Σ M ikpq ( x x x ε )∂ 2 xpxq u ε n k n m n i ∂ xm v dσ + ∂Σ M ikpq ( x x x ε )∂ 2 xpxq u ε n k Π iℓ ν ℓ v dl- Σ Π ℓm ∂ xm (M ikpq ( x x x ε )∂ xpxq u ε n k Π iℓ )v dσ +ε 2 Ω M ikpq ( x x x ε )∂ xpxq u ε ∂ x i x k v dx x x= Ω f v dx x x. (2.18)
Using the fact that we have imposed a homogeneous Dirichlet boundary condition and a zero normal double heat flux for the spatial retardation on Σ, we can see that the variational formulation (in vectorial form) reduces to the following: Find

u ε ∈ H 2 (Ω) ∩ H 1 0 (Ω) such that, Ω K( x x x ε )∇u ε .∇v dx x x+ε 2 Ω M( x x x ε )∇∇u ε :∇∇v dx x x= Ω f v dx x x, (2.19
∥u ε ∥ 2 H 1 (Ω) + ε 2 ∥∇∇u ε ∥ 2 L 2 (Ω,R 3×3 ) 1/2 ≤ c(Ω) ∥f ∥ L 2 (Ω) .
(2.20)

Variational formulation

The weak solution to (2.19) can be classified as the unique minimum of the functional

J ε (θ), u ε = arg min θ∈H 2 (Ω)∩H 1 0 (Ω) J ε (θ), (2.21)
where

J ε (θ):= 1 2 Ω K( x x x ε )∇θ.∇θ dx x x+ 1 2 Ω ε 2 M( x x x ε )∇∇θ:∇∇θ dx x x- Ω f θ dx x x. (2.22)
A standard computation of the variational derivative of J ε will recover (2.19) and the Euler-Lagrange equations in succession.

3 Main results

Homogenization via unfolding Γ-convergence

We define the following domain decompositions (see [START_REF] Ciorȃnescu | Éclatement périodique et homogénéisation[END_REF][START_REF] Damlamian | An elementary introduction to periodic unfolding[END_REF][START_REF] Ciorȃnescu | The periodic unfolding method in homogenization[END_REF][START_REF] Ciorȃnescu | The Periodic Unfolding Method. Theory and Applications to Partial Differential Problems[END_REF]):

K - ε := ℓ ∈ Z d | ε(ℓ + Y ) ⊂ Ω , Ω - ε := int ∪ ℓ∈K - ε ε(ℓ + Y ) , Λ - ε := Ω\Ω - ε . (3.1) ε ε Ω Ω - ε 0 ∈ R 3 ε x ε ∈ Z 3 εy ∈ εY x ∈ Ω
F : L p (Ω × Y ) → R (p > 1)
. We say that F ε unfolding Γ-converges to

F if for all u ∈ L p (Ω × Y ):
1. For every sequence

u ε ∈ L p (Ω) such that T ε (u ε ) ⇀ u in L p (Ω × Y ) one has, F (u) ≤ lim inf ε→0 F ε (u ε ); (3.4)
2. there exists a sequence

u ε ∈ L p (Ω) such that T ε (u ε ) ⇀ u in L p (Ω × Y ) and F (u) = lim ε→0 F ε (u ε ). (3.5)
Theorem 3.1. The sequence of functionals J ε unfolding Γ-converge in the weak topology of H 2 (Ω) ∩ H 1 0 (Ω) to J eff as ε → 0,

J ε Γ → J eff , (3.6)
where Moreover, let f :=div Griso's result in (1.9), we see that (3.11) is more accessible to computations, hence, more practical.

J eff (θ):= 1 2 Ω K eff ∇θ.∇θ dx x x - Ω f θ dx x x, (3.7
F F F , F F F ∈ L 3 (Ω, R 3 ) then, ∇u ε -∇u -∇ y û(•, • ε ) → 0 in L 2 (Ω, R 3 ), ( 3

Remarks on the homogenized coefficients

Symmetry

Since we work within a variational framework, the homogenized coefficients inherit the symmetry that is imposed on them from the framework. If one were to obtain the effective tensor through a multiple scale expansion then the tensor would have the following, nonsymmetric, form:

K eff ik := Y K ij (y y y)(δ jk - ∂w k ∂y j
) dy y y.

(3.12)

Naturally, the two forms, (3.8) and (3.12), are equivalent (since the tensors K and M are assumed to be isotropic). The proof is virtually identical to the classical case (see [START_REF] Sanchez-Palencia | Non-homogeneous media and vibration theory[END_REF]), however, we feel it should be included here since it is a higher-gradient generalization of the classical case:

Multiply with a test function (3.9) (3.15)

Variational characterization of the higher-order effective tensor

By (3.8), we can see that K eff is symmetric as well, and is determined in its entirety from the knowledge of the quadratic form K eff ξ ξ ξ • ξ ξ ξ for any constant vector ξ ξ ξ ∈ R 3 . Using definition (3.8) one can check that,

K eff ξ ξ ξ • ξ ξ ξ = Y K(y y y)(-ξ ξ ξ+∇ y w ξ ).(-ξ ξ ξ+∇ y w ξ ) dy y y + Y M(y y y)∇ y ∇ y w ξ :∇ y ∇ y w ξ dy y y, (3.16)
where w ξ is the solution of the following cell problem,

-div y (K (∇ y w ξ -ξ ξ ξ) -div y (M:∇ y ∇ y w ξ )) =0 in Y, y y y → w ξ (y y y) is Y -periodic. (3.17)
Using completely standard techniques stemming from the calculus of variations we can write for our case,

K eff ξ ξ ξ • ξ ξ ξ= inf v∈H 2 per (Y) Y K(y y y)(∇ y v-ξ ξ ξ).(∇ y v-ξ ξ ξ) dy y y+ Y M(y y y)∇ y ∇ y v:∇ y ∇ y v dy y y , (3.18)
with

H 2 per (Y):={v ∈ H 2 (Y) | v is Y -periodic}. It
(Ω) ∩ H 1 0 (Ω) to J eff as ε → 0, J ε Γ → J eff , (3.6)
where

J eff (θ):= 1 2 Ω K eff ∇θ.∇θ dx x x - Ω f θ dx x x, (3.7)
K eff ik := Y K(y y y)(∇ y w k -e e e k
).(∇ y w i -e e e i ) + M(y y y)∇ y ∇ y w k :∇ y ∇ y w i dy y y,

with w k the Y -periodic local solution to, -div y (K(y y y)(∇ y w k -e e e k ) -div y (M(y y y)

∇ y ∇ y w k )) =0 in Y. (3.9)
Proof. Define,

E ε (θ):= 1 2 a(θ, θ) if θ ∈ H 2 (Ω) ∩ H 1 0 (Ω) +∞ if θ ∈ L 2 (Ω)\H 2 (Ω) ∩ H 1 0 (Ω), (4.1)
with a(θ, θ):

= Ω K( x x x ε )∇θ.∇θ dx x x+ Ω ε 2 M( x x x ε )∇∇θ:∇∇θ dx x x, (4.2)
and

E eff (θ):= 1 2 Ω K eff ∇θ.∇θ dx x x if θ ∈ H 2 (Ω) ∩ H 1 0 (Ω) +∞ if θ ∈ L 2 (Ω)\H 2 (Ω) ∩ H 1 0 (Ω). (4.3)
It suffices to show that E ε Γ → E eff where the energy E ε is finite.

Γ-limit inferior inequality.

Let u ε ∈ H 2 (Ω) ∩ H 1 0 (Ω) be a solution to (2.11) then based on the estimates in (2.20) and properties of the periodic unfolding operator T ε (see [START_REF] Ciorȃnescu | The periodic unfolding method in homogenization[END_REF] for a definition) we have,

• T ε (u ε ) ⇀ u in L 2 (Ω, H 2 (Y)) • T ε (∇u ε ) ⇀ ∇ x u + ∇ y û in L 2 (Ω, H 1 (Y, R 3 )) • T ε (ε∇∇u ε ) ⇀ ∇ y ∇ y û in L 2 (Ω × Y, R 3×3 ).
Then, for all φ φ φ ∈ L 2 (Ω, C 1 per (Y, R 3 )) we have,

0≤ 1 2 Ω K( x x x ε )(∇u ε -φ φ φ).(∇u ε -φ φ φ) dx x x+ 1 2 Ω ε 2 M( x x x ε )(∇∇u ε -∇φ φ φ):(∇∇u ε -∇φ φ φ) dx x x. (4.4)
Opening up the parentheses in the above expression we obtain, Since the space L 2 (Ω, C 1 per (Y, R 3 )) is dense in L 2 (Ω × Y, R 3 ), the above inequality holds for a sequence of regular functions of the form φ φ φ ε (x x x, y y y) Y -periodic in y y y with the following convergence properties,

E ε (u ε )≥ Ω K( x x x ε )∇u ε .φ φ φ dx x x- 1 2 Ω K( x x x ε )φ φ φ.φ φ φ dx x x + Ω ε 2 M( x x x ε )∇∇u ε :∇φ φ φ dx x x- 1 2 Ω ε 2 M( x x x ε )
• T ε (φ φ φ ε ) ⇀ ∇ x u + ∇ y û in L 2 (Ω × Y, R 3 ) • T ε (ε∇φ φ φ ε ) ⇀ ∇ y ∇ y û in L 2 (Ω × Y, R 3×3 ).
Upon extracting a (a non-relabeled) sub-sequence, we obtain in the limit inferior of (4.5),

lim inf ε→0 E ε (u ε )≥ 1 2 Ω×Y K(y y y)(∇ x u+∇ y û).(∇ x u+∇ y û)
+M(y y y)∇ y ∇ y û:∇ y ∇ y û dy y ydx x x

≥ 1 2 Ω inf v∈H 2 per (Y) Y K(y y y)(∇ x u+∇ y v).(∇ x u+∇ y v) +M(y y y)∇ y ∇ y v:∇ y ∇ y v dy y ydx x x = 1 2 Ω K eff ∇ x u.∇ x u dx x x
Γ-limit superior inequality.

We will construct the recovery sequence for smooth functions initially and then use a diagonalization argument to complete the proof. Define the sequence,

u ε (x x x) = u(x x x) + εû( x x x ε ). (4.8)
Then,

• T ε (u ε ) → u in L 2 (Ω, H 2 (Y)) • T ε (∇u ε ) → ∇ x u + ∇ y û in L 2 (Ω, H 1 (Y, R 3 )) • T ε (ε∇∇u ε ) → ∇ y ∇ y û in L 2 (Ω × Y).
Thus,

E ε (u ε ) = 1 2 Ω K( x x x ε )∇u ε .∇u ε +ε 2 M( x x x ε
)∇∇u ε :∇∇u ε dy y ydx x x = 1 2 Ω×Y K(y y y)T (∇u ε ).T (∇u ε )+ε 2 M(y y y)T (∇∇u ε ):T (∇∇u ε ) dy y ydx x x.

(4.9) Hence, passing to the limit as ε → 0 in the expression above we obtain,

lim ε→0 E ε (u ε ) = Ω×Y (K(y y y)(∇ x u + ∇ y û).(∇ x u + ∇ y û)+M(y y y)∇ y ∇ y û):∇ y ∇ y û) dy y ydx x x = 1 2 Ω inf v∈H 2 per (Y) Y K(y y y)(∇ x u+∇ y v).(∇ x u+∇ y v)+M(y y y)∇ y ∇ y v:∇ y ∇ y v dy y ydx x x = 1 2 Ω K eff ∇ x u.∇ x u dx x x (4.10)
We can conclude the proof using a density and diagonalization argument. for all ϕ ∈ H 1 0 (Ω). In the sense of distributions we can recover,

-div K eff ∇u = f in D(Ω), u = 0 on ∂Ω, (4.12)
with K eff given in (3.8).

The above equation is precisely the limit problem that has been obtained countless times in the homogenization literature. The difference in this work is that the local solution used to construct the effective tensor K eff belong satisfy a higher-gradient problem with w k ∈ H 2 per (Y ) (and not only in H 1 (Y )). This newly found local solution regularity is combined with the theorem below to prove the classical corrector convergence result.

The next theorem is included only for completion purposes and is a consequence of G.

Stampacchia's interpolation theorem [START_REF] Stampacchia | The spaces L (p,λ) , N (p,λ) and interpolation[END_REF]. Practically, it allows one to control the L p norm of the gradient of the weak solution of a second order elliptic system with constant coefficients, for p ∈ [2, ∞), by controlling the integrability of the data. We refer the interested reader to [3, Thm. 3.28, pg. 76] (and references therein) for a modern treatment on the subject. Our aim here is to apply the theorem below to the second order elliptic system with constant coefficients in (4.12) to obtain the necessary integrability for the homogenized solution to compensate for the remaining regularity that is needed in the classical corrector convergence theorem.

Theorem 4.1. Let u ∈ H 1 0 (Ω, R 3 ) be a weak solution of the Dirichlet problem

-div(A∇u) = -div F F F , (4.13)
where the constant coefficients A ijαβ satisfy the Legendre-Hadamard condition, and

F F F ∈ L p (Ω, R 3 ) for some 2 ≤ p < ∞. Then, ∇u ∈ L p (Ω, R 3 ) and ∥∇u∥ L p (Ω,R 3 ) ≤ c ∥F F F ∥ L p (Ω,R 3 ) (4.14)
A proof of this theorem can be found in [START_REF] Giaquinta | An Introduction to the Regularity Theory for Elliptic Systems, Harmonic Maps and Minimal Graphs[END_REF][Thm. (4.17)

:=div F F F , F F F ∈ L 3 (Ω, R 3 ) then, ∇u ε -∇u -∇ y û(•, • ε ) → 0 in L 2 (Ω, R 3 ), ( 3
On the last step we made use of equation (4.10).

Discussion

Lastly, remark that it has been observed in [START_REF] Nika | A gradient system for a higher-gradient generalization of Fourier's law of heat conduction[END_REF] that the above generalized heat conduction problem can be formulated as a classical gradient system with an induced gradient flow equation,

∂ t u = +∇ G DS (u) in Ω, (5.1)
with G the metric tensor and L (u) := Ω log u dx x x the entropy. Naturally, assumptions have to be made on the coefficients so that model is thermodynamically correct. Although, we will address the time-dependent case in future work, we would like to point out that this added gradient structure could be used to pass to the limit in the time dependent case through evolutionary Γ-convergence (see e.g. [START_REF] Mielke | On evolutionary Γ-convergence[END_REF]).

  tangential and normal components of ∇ξ are,(∇ξ) τ = -n n n × (n n n × ∇ξ)=∇ξ -(∇ξ) n n n n, (∇ξ) n :=∇ξ • n n n.(2.13) Moreover, we introduce the surface gradient of ξ using the projection operator Π:=I-n n n⊗n n n, ∇ s ξ=(I -n n n ⊗ n n n)∇ξ=Π∇ξ.(2.14)

Figure 1 : 3 . 2 )( 3 . 3 )

 13233 Figure 1: Unfolding operator on a periodic grid Let [z z z] Y = (⌊z 1 ⌋, ⌊z 2 ⌋, ⌊z 3 ⌋) denote the integer part of z z z ∈ R 3 and denote by {z z z} Y the difference z z z-[z z z] Y which belongs to Y . Regarding our multiscale problem that depends on a small length parameter ε > 0, we can decompose any x x x ∈ R 3 using the maps [•] Y : R 3 → Z 3 and {•} Y : R 3 → Y the following way (see Fig. 3.1 (right)), x x x = ε x x x ε Y + x x x ε Y . (3.2)

. 9 ) 3 . 1 .Theorem 3 . 2 .

 93132 y y)(∇ y w k -e e e k ).(∇ y w i -e e e i ) + M(y y y)∇ y ∇ y w k :∇ y ∇ y w i dy y y,(3.8) with w k the Y -periodic local solution to,-div y (K(y y y)(∇ y w k -e e e k ) -div y (M(y y y)∇ y ∇ y w k )) =0 in Y. (3RemarkThe first part of (3.8) is what one would expect as a result from the homogenization of second order linear elliptic equations while the second part is new, specific to the inclusion of the second-gradient thermal effects. Let u ε and u be the unique minimizers of (2.22) and (3.7), respectively.

. 10 ). 11 ) 3 . 2 .

 101132 as ε → 0 with û(x x x, y y y) = -w k (y y y) ∂u ∂x k (x x x) + ũ1 (x x x). (3Remark We note in the above theorem, the assumptions that ∇ y w k ∈ L r and ∇u ∈ L s with 1/r + 1/s = 1/2 are no longer needed. Furthermore, if we compare (3.11) with G.

1 Theorem 3 . 1 .

 131 follows immediately that the effective coefficients are elliptic (satisfy Legendre condition), i.e. K eff ξ ξ ξ • ξ ξ ξ ≥ K|ξ ξ ξ| 2 where K = Y K(y y y) dy y y. 4 Proof of the main results 4.1 Proof of Theorem 3.The sequence of functionals J ε unfolding Γ-converge in the weak topology of H 2

  ∇φ φ φ:∇φ φ φ dx x x = Ω×Y K(y y y)T ε (∇u ε ).T ε (φ φ φ) dy y ydx x x -1 2 Ω×Y K(y y y)T ε (φ φ φ).T ε (φ φ φ) dy y ydx x x + Ω×Y ε 2 M(y y y)T ε (∇∇u ε ):T ε (∇φ φ φ) dy y ydx x x-1 2 Ω×Y ε 2 M(y y y)T ε (∇φ φ φ):T ε (∇φ φ φ) dy y ydx x x ε→0 → Ω×Y K(y y y)(∇ x u + ∇ y û).φ φ φ dy y ydx x x-1 2 Ω×Y K(y y y)φ φ φ.φ φ φ dy y ydx x x + Ω×Y M(y y y)∇ y ∇ y û:∇ y φ φ φ dy y ydx x x-1 2 Ω×Y M(y y y)∇ y φ φ φ:∇ y φ φ φ dy y ydx x x(4.5)

Let u u u ∈ H 2 (

 2 Ω) ∩ H 1 0 (Ω). Using the density of C ∞ 0 (Ω) in H 2 (Ω) ∩ H 1 0 (Ω) we can suppose without loss of generality that u ∈ C ∞ 0 (Ω). Furthermore, let û ∈ C 2 per (Y) be a minimizer of, inf v∈H 2 per (Y) Y K(y y y)(∇ x u+∇ y v).(∇ x u+∇ y v) dy y y+M(y y y)∇ y ∇ y v:∇ y ∇ y v . (4.6) Existence of such a minimizer is shown using classical arguments of coercivity and lower semi-continuity of, y y)(∇ y v-ξ ξ ξ).(∇ y v-ξ ξ ξ) dy y y+ Y M(y y y)∇ y ∇ y v:∇ y ∇ y v dy y y . (4.7)

Remark 4 . 1 .

 41 The Euler-Lagrange equation of (3.7) is, Ω K eff ∇u.∇ϕ dx x x = Ω f ϕ dx x x (4.11)

. 10 )

 10 as ε → 0 with û(x x x, y y y) = -w k (y y y) ∂u ∂x k (x x x) + ũ1 (x x x).(3.11)Proof. The proof of this theorem amounts to essentially a compensated regularity argument, where the burden of regularity is split (un-equally) between the local solution and the homogenized solution. We argue as follows: by the assumption on the body force f we can apply Theorem 4.1 and obtain ∇u ∈ L 3 (Ω, R 3 ). Moreover, since the local solution w k ∈ H 2 per (Y ) we extend it by periodicity to the entire space. By applying standard Sobolev embedding theory we can obtain ∇ y w k ∈ L 6 (R 3 , R 3 ) (see e.g. [14, Theorem 2.31]). Therefore, the above compensated regularity argument makes the expression ∇u ε -∇u -∇ y û(•, • ε ) well defined in the L 2 norm and hence, we can prove Proposition 9.12 in [8, pg. 185] directly. In what follows, we provide the actual steps of the argument that is given in [8, Proposition 9.12, pg. 185] with the appropriate adjustments that the second-gradient thermal effects introduce.Staring from the ellipticity of the second order tensor K we have,c 1 ∇u ε (x x x)-∇u(x x x) -∇ y û(x x x, ε (x x x)-∇u(x x x)-∇ y û(x x x, x x x ε )).(∇u ε (x x x)-∇u(x x x)-∇ y û(x x x, ε (x x x) -∇u(x x x) -∇ y û(x x x, x x x ε )).(∇u ε (x x x) -∇u(x x x) -∇ y û(x x x, ε -∇ y ∇ y û):(ε∇∇u ε -∇ y ∇ y û) dx x x (4.15)where the last inequality is due to the ellipticity condition of the tensor M. Hence, opening up the expression above, using the structure of equation (2.11), and the symmetry of the tensors K and M we have, c 1 ∇u ε (x x x)-∇u(x x x)-∇ y û(x x x, y û).(∇u+∇ y û)dx x x ε :∇ y ∇ y ûdx x x+ Ω M( x x x ε )∇ y ∇ y û:∇ y ∇ y ûdx x x.

(4. 16 )c 1

 161 By unfolding, we can see that we are able to pass to the limit in each expression due to estimate (2.20). Thus,lim ε→0 ∥∇u ε -∇u-∇ y û∥ 2 L 2 (Ω,R 3 ) ≤ Ω×Y fudy y ydx x x -Ω×Y K(y y y)(∇u+∇ y û).(∇u+∇ y û)dy y ydx x x -Ω×Y M(y y y)∇ y ∇ y û:∇ y ∇ y ûdy y ydx x x=0.

  .17)Since the varational formulation of(2.11) is not a standard one, we write down the details for the readers convenience using the notation introduced in Section 2.3.2. Hence, if we multiply (2.11) by a test function v ∈ {ϕ ∈ C ∞ (Ω) | ϕ = 0 on Σ} and integrate by parts several times (including integration by parts on surfaces using formula (2.15)) we obtain,

  If we add and subtract coordinate y i on the w i term of the second integral we obtain,

										and integrate by parts to obtain,
		0 =	Y	M mjpq	∂ 2 w k ∂y p ∂y q	∂ 2 v ∂y m ∂y j	dy y y +	Y	K mj	∂(w k -y k ) ∂y j	∂v ∂y m	dy y y.	(3.13)
	Selecting v = w i we have,
		0 =	Y	M mjpq	∂ 2 w k ∂y p ∂y q	∂ 2 w i ∂y m ∂y j	dy y y +	Y	K mj	∂(w k -y k ) ∂y j	∂w i ∂y m	dy y y.	(3.14)
	Y	K ij (y y y)(δ jk -	∂w k ∂y j	) dy y y
		=	Y	M mjpq	∂ 2 w k ∂y p ∂y q	∂ 2 w i ∂y m ∂y j	dy y y +

Y K mj ∂(w k -y k ) ∂y j ∂(w i -y i )

∂y m dy y y.

Acknowledgements

The authors gratefully acknowledges the financial support by the Knowledge Foundation (project nr. KK 2020-0152).