Variational-asymptotic homogenization of thermoelastic periodic materials with
thermal relaxation

Deison Préve1,111Corresponding authors: Tel: +39 0103352556 – Fax: +39
0103352517
E-mail addresses<EMAIL_ADDRESS><EMAIL_ADDRESS>Andrea
Bacigalupo2,11footnotemark: 1, Marco Paggi1
1 IMT School for Advanced Studies Lucca, Piazza San Francesco 19, 55100 Lucca,
Italy
2 DICCA, Università degli Studi di Genova, Via Montallegro 1, 16145 Genova,
Italy

Abstract

A multiscale asymptotic homogenization method for periodic microstructured
materials in presence of thermoelasticity with periodic spatially dependent
one relaxation time is introduced. The asymptotic expansions of the micro-
displacement and the micro-temperature fields are rewritten on the transformed
Laplace space and expressed as power series of the microstructural length
scale, leading to a set of recursive differential problems over the periodic
unit cell. The solution of such cell problems leads to the perturbation
functions. Up-scaling and down-scaling relations are then defined, and the
latter allow expressing the microscopic fields in terms of the macroscopic
ones and their gradients. Average field equations of infinite order are also
derived. The efficiency of the proposed technique was tested in relation to a
bi dimensional orthotropic layered body with orthotropy axis parallel to the
direction of the layers, where the mechanical and temperature constitutive
properties were well stabilised. The dispersion curves of the homogenized
medium, truncated at the first order are compared with the dispersion curves
of the heterogeneous continuum obtained by the Floquet-Bloch theory. The
results obtained with the two different approaches show a very good agreement.

Keywords: Dynamic variational-asymptotic homogenization, Periodic materials,
Generalized thermoelasticity, One thermal relaxation time, Wave propagation.

## Introduction

The thermoelasticity phenomenon bonds together a set of theories as [1]
quoted, namely, the general theory of heat conduction, thermal stresses, and
strains set up by thermal flow in elastic bodies, and the reverse effect of
temperature distribution caused by the elastic deformation itself leading to
thermoelastic energy dissipation. The main theory of thermoelasticity based on
Fourier law of heat conduction, leaves thermal perturbations propagate at
infinite velocity in a diffusive manner, once the coupled displacement-
temperature governing equation is a parabolic-type partial differential
equation [1, 2, 3]. In practical terms, this means that when a temperature
gradient is suddenly produced in some region of the sample, this entails in an
instantaneous disturbance at each point of the material [4, 5]. However, there
are some experimental observations, as [6] pinpointed, where the temperature
acts like a wave propagating through the body with finite speed, commonly
referred to as second sound [7, 8]. The pattern of the heat wave propagation
has been observed in superfluids, inhomogeneous materials like sand and
processed meat, [9, 10], and in pure crystals, [11]. Moreover, this factor
contradicts the physical affirmation that for a finite time interval a
disturbance of bounded support may only generate the response of a bounded
support [12, 13, 6].

In order to overcome such a paradox, recent developments were made as [14]
mentions, and a couple of generalized thermoelastic theories have flourished.
[15] introduced the theory of generalized thermoelasticity with one relaxation
time, originally proposed by [16] in the context of theory of gases, and later
by [17] in the context of heat conduction in rigid bodies [6]. In this theory,
a modified law of heat conduction including both the heat flux and its time
derivative replaces the conventional Fourier law [18]. The heat equation
associated to this theory becomes hyperbolic and hence eliminates the paradox
of the propagation of thermal signals with an infinite speed. The equations of
motion and constitutive relations, remain the same for both theories of
thermoelasticity, the uncoupled theory and the coupled one [15].

Modelling multiphase materials containing a periodic microstructure, where
additionally elasticity and heat transfer are combined phenomena, is a topic
of vast importance in modern applications, such as aerospace, aircraft,
biomedical and electronic [19, 20, 21, 22, 23]. Nonetheless, a plethora of
problems are assumed on heterogeneous continuum with a small accuracy over the
microstructural length scale. Deriving the solution of the governing
thermoelastic partial differential equation with one relaxation time may be
cumbersome both analytically and numerically, due the periodicity of the media
[24, 25]. The multiscale asymptotic homogenization approaches, well stabilised
in earlier works by [26, 27, 28, 29], show themselves as effective tools for
determining the responses of the microscopic phases on the overall properties
of the composites, by replacing an heterogeneous continuum by an equivalent
homogenized model whose solutions are good approximations of the real ones,
but are characterized by constitutive tensors not affected by the fast
variable which gives rapid oscillations due to the underlying microstructure.
Therefore, such procedure might be quite effective computationally [30, 31].

At this stage, it must be highlighted that various homogenization methods have
being applied to study the overall properties of multiphase periodic
materials, which may be classified in asymptotic technique [32, 33, 34, 35,
36], variational asymptotic schemes [37, 38, 39, 31, 40], analytical paths
[41, 42, 43, 44, 45], and computational approaches [46, 47, 48, 49, 50, 34,
51, 52, 53, 54]. In what concerns to periodic materials in presence of
thermoelastic effects, as seen in [55], the homogenization method has been
investigate therein, and developed in the context of thermoelastic periodic
material with pertinent physical applications highlighted in [56].
Furthermore, a variational-asymptotic technique for thermoelastic periodic
materials was brought by [57].

As a matter of fact, cutting edge research regarding asymptotic homogenization
techniques, variational-asymptotic approaches and computational methods are
being made along several multiphase materials in presence of certain
phenomena. For instance, asymptotic homogenization technique over
piezoelectric composite materials may be found in [58, 59, 60, 61] and on
thermal-piezoelectric materials with a periodic microstructure via asymptotic
schemes can be checked in [62]. Also, multiscale homogenization schemes have
been applied to characterize the behaviour and the global constitutive
properties of viscoelastic heterogeneous materials with periodic
microstructure [63, 64, 65, 66]. Particularly, in the latter case, the
variational-asymptotic homogenization method was proposed in order to
characterize the propagation of dispersive waves in viscoelastic materials
with periodic microstructure. Likewise, the homogenization approach for
describing the elastic, thermal and diffusive properties of periodic materials
(also on periodic layered materials) in presence of thermal-diffusion has been
explored by [67, 68, 69]. Effective analysis on heterogeneous continuum in
presence of thermal-mechanical and thermal-magneto-electro-elastic
deformations have been studied by [70, 71, 72] and [73], respectively.

Inspired by those studies where the heterogeneous continuum is rather bartered
by some homogeneous model led by an homogenization procedure, the herein
proposed original theoretical framework combines the generalized theory of
thermoelasticity with periodic spatially dependent one time relaxation with
the homogenization technique in periodic composite materials. Having said
that, the goal of this work is to present a novel method to characterize the
overall properties, namely elastic and thermal, of multiphase periodic
materials governed by Lord-Shulman generalized thermoelasticity equation via
asymptotic homogenization. In the following, the Lord-Shulman generalized
thermoelastic governing equations at the micro-scale, which describe the non-
homogeneous medium, are determined in the time domain and subsequently
transformed by the bi-lateral Laplace transform once the relaxation time
varies within the material phase, in order to separate it from the partial
derivative in time contained in the differential operator induced by Maxwell-
Cattaneo law. Such a procedure hence takes the input to the Laplace domain
(complex frequency domain). With the aim of separating the fast from the slow
variable, the micro-displacement and micro-temperature fields are rewritten as
power series expansions of the microstructural length scale, such that a
cascade of recursive non-homogeneous differential problems is defined over the
periodic unit cell. From this point, the solvability conditions are imposed to
these recursive differential problems, arising from the down-scaling relations
which are written in terms of the perturbation functions. From the
substitution of the down-scaling relations into the microscopic field
equations, the governing equations at infinite order are given in closed form
in terms of the microscopic constitutive properties and the perturbation
functions. Finally, by truncating the governing equations, the first order
homogenized field equations are derived from which the overall constitutive
tensors are derived. The capability of the method is then assessed through the
study of thermal wave propagation over a two-dimensional bi-layered media with
a periodic microstructure with orthotropic phases and axis of orthotropy
parallel to the direction of layering, where the overall mechanical and
thermal properties are analytically determined. The outcome dispersion
functions from the homogenized model are therefore set to be compared with
those of the heterogeneous thermoelastic material obtained from the solutions
of the Floquet–Bloch theory [74, 75, 76].

The work is sorted out as follows: Sec. 1 recalls the fundamentals of
thermoelasticity with one relaxation time in a periodic media, describes the
microscopic governing equations in the time domain and in the Laplace space.
Sec. 2 is dedicated in developing the cascade of differential problems as well
as their solutions. In Sec. 2.1 is derived the cell problems and their
respective perturbation functions. The down-scaling and up-scaling relations
are defined in Sec. 2.2. The variational-asymptotic to establish an
equivalence between the equations at macro-scale and micro-scale procedure and
the zeroth order truncation is attacked in Sec. 3. Sec. 4 reports wave
propagation and the related dispersion functions. Throughout Sec. 5, the
herein method is applied for studying the overall properties of two-
dimensional bi-layered orthotropic composites. Specifically, the overall
elastic, thermal dilatation and thermal conduction tensors are determined in
theirs analytic form in Sec. 5.1, whereas benchmark and some analysis are
pursued in Sec. 5.2. Further mathematical details about (i) asymptotic
expansions of the field equations ( Sec. A.1.), (ii) recursive differential
equations (Sec. A.2.), (iii) the average field equation of infinite order
(A.3.), (iv) power-like functional (Secs. A.4., A.5.), (v) frequency spectrum
via Floquet-Bloch theory within the heterogeneous continuum (Secs. B.1., B.2.,
B.3.), are vastly exploited in the Supplementary Material.

##  1 Derivation of thermoelasticity with one relaxation time on the Laplace
domain

Throughout this Section, the generalized Lord-Shulman thermoelasticity theory
is recalled, in relation to an heterogeneous periodic composite material,
which will be taken to the Laplace domain.

Let consider an heterogeneous composite material $\mathfrak{L},$ Fig. 1, under
the assumption of small strains (for instance a stretching load), which leads
to a process of exchanging mechanical energy into thermal energy under the
action of externally applied thermal-mechanical loadings. Such procedure is
then followed by temperature variances and strains within the body, all of
which vanish upon the removal of the mentioned thermal-mechanical loadings.
The continuum $\mathfrak{L}$ is described as a linear thermoelastic Cauchy
medium [12, 13, 6] under stresses induced by body forces and temperature
changes due to heat source.

On each point of the material is endowed with a displacement field
$\boldsymbol{u}(\textbf{x},t)=u_{i}\boldsymbol{e}_{i}$ and a relative
temperature field $\theta(\textbf{x},t)=T(\textbf{x},t)-T_{0},$ where
$T(\textbf{x},t)$ is the absolute temperature, $T_{0}$ is a reference stress-
free temperature. The coupled constitutive relations link the stress tensor
$\boldsymbol{\sigma}(\textbf{x},t)=\sigma_{ij}\boldsymbol{e}_{i}\otimes\boldsymbol{e}_{j},$
the heat flux vector
${\mathbf{q}}(\textbf{x},t)=q_{ij}\boldsymbol{e}_{i}\otimes\boldsymbol{e}_{j},$
and the entropy per unit of volume
$\boldsymbol{\eta}(\textbf{x},t)=\eta_{ij}\boldsymbol{e}_{i}\otimes\boldsymbol{e}_{j},$
to the aforementioned relevant fields $\boldsymbol{u}(\textbf{x},t),$
$\theta(\textbf{x},t),$ that is

$\displaystyle{\mathbf{q}}(\textbf{x},t)+\tau(\textbf{x})\dfrac{\partial\mathbf{q}(\textbf{x},t)}{\partial
t}=-{\bar{\mathbf{K}}}(\textbf{x})\nabla\theta(\textbf{x},t)\,,$ (1a)
$\displaystyle\boldsymbol{\eta}(\textbf{x},t)=\boldsymbol{\alpha}(\textbf{x})\boldsymbol{\varepsilon}(\textbf{x},t)+\dfrac{C_{E}}{T_{0}}\theta(\textbf{x},t)\,,$
(1b)

where
$\boldsymbol{\varepsilon}(\textbf{x},t)=\varepsilon_{ij}\boldsymbol{e}_{i}\otimes\boldsymbol{e}_{j}$
is the strain tensor,
$\textswab{C}=C_{ijhk}\textbf{e}_{i}\otimes\textbf{e}_{j}\otimes\textbf{e}_{h}\otimes\textbf{e}_{k}$
is the fourth-order elastic tensor,
${\bar{\mathbf{K}}}=\bar{K}_{ij}\textbf{e}_{i}\otimes\textbf{e}_{j}$ is the
second-order thermal conductivity tensor,
${{\boldsymbol{\alpha}}}=\alpha_{ij}\textbf{e}_{i}\otimes\textbf{e}_{j}$ is
the second-order stress-temperature tensor, $C_{E}$ is the specific heat at
zero strain, and for convenience we define $p=C_{E}/T_{0}.$ Besides, as
aforementioned, to avoid the physical paradox of infinite speed for the
propagation of heat signal in the classical thermoelastic equations [13],
which are parabolic, a space dependent relaxation time $\tau$ is introduced,
known as the Maxwell-Cattaneo law replacing the Fourier law for thermal
conduction, in order to transform them into hyperbolic equations.

As the material is under the effect of small displacements, then the micro-
strain tensor is defined as

$\boldsymbol{\varepsilon}(\textbf{x},t)=\dfrac{1}{2}\left(\nabla\boldsymbol{u}(\textbf{x},t)+\nabla^{T}\boldsymbol{u}(\textbf{x},t)\right),$
(1b)

where $\nabla\boldsymbol{u}$ is the gradient of the micro-displacement
$\boldsymbol{u}(\textbf{x},t).$

The following balance equations hold

$\displaystyle\nabla\cdot(\boldsymbol{\sigma}(\textbf{x},t))+\textbf{b}(\textbf{x},t)=\rho(\textbf{x})\ddot{{\mathbf{u}}}(\textbf{x},t),$
(1ca)
$\displaystyle\nabla\cdot({\mathbf{q}}(\textbf{x},t))-\bar{r}(\textbf{x},t)=-T_{0}\dot{\boldsymbol{\eta}}(\textbf{x},t)\,,$
(1cb)

with $\textbf{b}=b_{i}\textbf{e}_{i}$ being the body forces, $\bar{r}$ the
heat source per unit time per unit volume and $\rho$ is the mass density.
Throughout, $t$ stands for the time coordinate, the superimposed dot denotes
the derivative in relation to time, and consider $\bar{r}=r/T_{0}.$

Figure 1: $3D$ continuum domain $\mathfrak{L}$ with periodic cell
$\mathcal{A}$ and its dimensionless unit cell $\mathcal{Q}$.

With the thermoelasticity constitutive equations and equilibrium equations
well established, we may thus introduce the periodic microstructure of an
heterogeneous composite material theory [77]. Let a point be identified by the
position vector
$\textbf{x}=x_{1}\boldsymbol{e}_{1}+x_{2}\boldsymbol{e}_{2}+x_{3}\boldsymbol{e}_{3}\in\mathbb{R}^{3}$
related to a system of coordinates with origin at point $O$ and written in
terms of the orthogonal base
$\\{\boldsymbol{e}_{1},\boldsymbol{e}_{2},\boldsymbol{e}_{3}\\},$ on the
three-dimensional heterogeneous material $\mathfrak{L}$ characterized by a
certain periodic microstructure. Let also
$\mathcal{A}=[0,\varepsilon]\times[0,\delta\varepsilon]\times[0,\gamma\varepsilon]$
be a periodic cell with characteristic size $\varepsilon$ defined by three
orthogonal periodic vectors $\boldsymbol{v}_{1}$, $\boldsymbol{v}_{2}$ and
$\boldsymbol{v}_{3},$ written as
$\boldsymbol{v}_{1}=d_{1}\boldsymbol{e}_{1}=\varepsilon\boldsymbol{e}_{1}$,
$\boldsymbol{v}_{2}=d_{2}\boldsymbol{e}_{2}=\delta\varepsilon\boldsymbol{e}_{2}$
and
$\boldsymbol{v}_{3}=d_{3}\boldsymbol{e}_{3}=\gamma\varepsilon\boldsymbol{e}_{3}$.
The material domain is set up by attaching three-dimensionally the cell
$\mathcal{A}$ in accordance to the directions $\boldsymbol{v}_{1}$,
$\boldsymbol{v}_{2}$ and $\boldsymbol{v}_{3}$, see Fig. 1.

Due to the periodicity of the material, throughout the work the superscript m
refers to the micro-scale. Therefore,
$\textswab{C}^{m}(\textbf{x},t)=C_{ijhk}^{m}\textbf{e}_{i}\otimes\textbf{e}_{j}\otimes\textbf{e}_{h}\otimes\textbf{e}_{k}$
is the fourth-order micro-elastic tensor,
$\bar{\mathbf{K}}^{m}(\textbf{x},t)=\bar{K}_{ij}^{m}\textbf{e}_{i}\otimes\textbf{e}_{j}$
is the second-order micro-thermal-conductivity tensor,
$\boldsymbol{\alpha}^{m}(\textbf{x},t)=\alpha_{ij}^{m}\textbf{e}_{i}\otimes\textbf{e}_{j}$
is the second-order micro-stress temperature tensor, $\rho^{m}(\textbf{x})$ is
the material density, $p^{m}(\textbf{x})$ is the specific heat at zero micro
strain, and $\tau^{m}(\textbf{x})$ is the micro relaxation time. Furthermore,
they obey the following conditions

$\displaystyle\textswab{C}^{m}(\textbf{x}+\boldsymbol{v}_{i})=\textswab{C}^{m}(\textbf{x}),$
(1da)
$\displaystyle\bar{\mathbf{K}}^{m}(\textbf{x}+\boldsymbol{v}_{i})=\bar{\mathbf{K}}^{m}(\textbf{x}),$
(1db)
$\displaystyle\boldsymbol{\alpha}^{m}(\textbf{x}+\boldsymbol{v}_{i})=\boldsymbol{\alpha}^{m}(\textbf{x}),$
(1dc)
$\displaystyle\rho^{m}(\textbf{x}+\boldsymbol{v}_{i})=\rho^{m}(\textbf{x}),$
(1dd) $\displaystyle p^{m}(\textbf{x}+\boldsymbol{v}_{i})=p^{m}(\textbf{x}),$
(1de)
$\displaystyle\tau^{m}(\textbf{x}+\boldsymbol{v}_{i})=\tau^{m}(\textbf{x}),$
(1df)

where $i=1,2,3,$ for all $\textbf{x}\in\mathcal{A}.$

In the following derivation, let us define $L$ and $\varepsilon$ as the
macroscopic length and the microstructural length, respectively, where
$L>$$>\varepsilon,$ which gives the condition where the scales might be
separated since the macroscopic length is taken to be much larger than the
microstructural one. Let us also suppose that the heterogeneous material
undergoes to $\mathfrak{L}$-periodic body forces
$\boldsymbol{b}(\textbf{x},t),$ having zero mean values over the continuum
$\mathfrak{L}=[0,L]\times[0,\delta L]\times[0,\gamma L],$ hence $\mathfrak{L}$
defined as above represents a portion of the continuum. Now, regarding the
dimensionless cell, let us rescale the periodic cell $\mathcal{A}$ by a factor
equals to the characteristic length $\varepsilon,$ this implies that the non-
dimensional cell model of the periodic microstructure is then
$\mathcal{Q}=[0,1]\times[0,\delta]\times[0,\gamma].$ Moreover, from such
rescaling process it arises two variables namely, the macroscopic (or slow)
variable $\textbf{x}\in\mathcal{A}$ in charge of measuring the slow
oscillations within the continuum, and the microscopic (or fast) variable
$\boldsymbol{\xi}=\textbf{x}/\varepsilon\in\mathcal{Q}$ responsible in
evaluating the fast heat propagation wave within the composite [28, 32, 29,
78]. Thanks to cell $\mathcal{Q},$ the properties presented in (1d) are now
redefined as $\mathcal{Q}$-periodic over $\mathcal{Q}$ so the terms can be
expressed by the microscopic variable $\boldsymbol{\xi}$ as

$\displaystyle\vskip 6.0pt plus 2.0pt minus
2.0pt\textswab{C}^{m}(\textbf{x})=\textswab{C}^{m}\left(\mathbf{x},\boldsymbol{\xi}=\dfrac{\mathbf{x}}{\varepsilon}\right),$
(1ea) $\displaystyle\vskip 6.0pt plus 2.0pt minus
2.0pt\bar{\mathbf{K}}^{m}(\textbf{x})=\bar{\mathbf{K}}^{m}\left(\mathbf{x},\boldsymbol{\xi}=\dfrac{\mathbf{x}}{\varepsilon}\right),$
(1eb) $\displaystyle\vskip 6.0pt plus 2.0pt minus
2.0pt\boldsymbol{\alpha}^{m}(\textbf{x})=\boldsymbol{\alpha}^{m}\left(\mathbf{x},\boldsymbol{\xi}=\dfrac{\mathbf{x}}{\varepsilon}\right),$
(1ec) $\displaystyle\vskip 6.0pt plus 2.0pt minus
2.0pt\rho^{m}(\textbf{x})=\rho^{m}\left(\mathbf{x},\boldsymbol{\xi}=\dfrac{\mathbf{x}}{\varepsilon}\right),$
(1ed) $\displaystyle\vskip 6.0pt plus 2.0pt minus
2.0ptp^{m}(\textbf{x})=p^{m}\left(\mathbf{x},\boldsymbol{\xi}=\dfrac{\mathbf{x}}{\varepsilon}\right),$
(1ee) $\displaystyle\vskip 6.0pt plus 2.0pt minus
2.0pt\tau^{m}(\textbf{x})=\tau^{m}\left(\mathbf{x},\boldsymbol{\xi}=\dfrac{\mathbf{x}}{\varepsilon}\right).$
(1ef)

Along with the $\mathcal{Q}$-periodicity assumptions made for the relations
(1d), regarding the micro-scale, the governing equations are finally obtained
by plugging the micro-scale constitutive equations $(\ref{sigma})$ to
$(\ref{microstrain})$ into the micro-scale balance equations $(\ref{bal1})$
and $(\ref{bal2}),$ also hiding the arguments for a cleaner notation, one
provides

$\displaystyle\nabla\cdot\left(\textswab{C}^{m}\nabla{\mathbf{u}}-{{\boldsymbol{\alpha}^{m}}\theta}\right)+{\mathbf{b}}=\rho^{m}{\mathbf{\ddot{u}}}\,,$
(1fa)
$\displaystyle\mathbf{q}+\tau^{m}\mathbf{\dot{q}}=-(T_{0})^{-1}\bar{\mathbf{K}}^{m}\nabla\theta\,,$
(1fb)
$\displaystyle\nabla\cdot{\mathbf{{q}}}+{\boldsymbol{\alpha}^{m}}\nabla{\mathbf{\dot{u}}}+p^{m}\dot{\theta}=r\,.$
(1fc)

In the particular scenario where the relaxation time
$\tau^{m}(\mathbf{x},\boldsymbol{\xi})=\tau>0$ is assumed as a constant, the
governing equations $(\ref{gov1})$ and $(\ref{gov2})$ result into:

$\displaystyle\nabla\cdot\left({\textswab{C}^{m}\nabla{\mathbf{u}}}-{\boldsymbol{\alpha}^{m}\theta}\right)+{\mathbf{b}}=\rho^{m}\mathbf{\ddot{u}}\,,$
(1ga)
$\displaystyle\nabla\cdot\left({{-(T_{0})^{-1}\bar{\mathbf{K}}^{m}}}\nabla\theta\right)+{{\boldsymbol{\alpha}^{m}}}\mathcal{D}(\nabla\mathbf{\dot{u}})+p^{m}\mathcal{D}(\dot{\theta})=\mathcal{D}(r)\,,$
(1gb)

where the differential operator $\mathcal{D}=1+\tau^{m}\partial/\partial t$
depending on the relaxation time $\tau^{m}$ is introduced.

At this point it must be highlighted that, since the periodic cell
$\mathcal{Q}$ is composed by two different phase materials, the relaxation
time $\tau^{m}$ must be considered as space dependent, which means
$\tau(\mathbf{x},\boldsymbol{\xi})$ is a $\mathcal{Q}$-periodic function and
hence varies within the composite material. For this reason, rather than
continuing with the equation $(\ref{gov2}),$ we apply the time bilateral
Laplace transform, or analogously the Fourier transform with complex argument
[79, 80] to the equations $(\ref{q}),$ $(\ref{eta}),$ $(\ref{gov1})$ and
$(\ref{gov3}),$ in order to separate the time relaxation $\tau^{m}$ from the
partial derivative in time contained in the differential operator
$\mathcal{D}.$

To clarify ideas, the time bilateral Laplace transform of a function
$f:\mathbb{R}\rightarrow\mathbb{R}$ is defined as [79],

$\mathcal{L}\left\\{f(t)\right\\}=\hat{f}(s)=\int_{-\infty}^{+\infty}f(t)\mathrm{e}^{-st}\mathrm{d}t,\quad
s\in\mathbb{C},$ (1h)

where the Laplace argument $s$ and the Laplace transform $\hat{f}$ are complex
valued (i.e. $\hat{f}:\mathbb{C}\to\mathbb{C}$ ). Whereas from an analytical
standing point in [79], the inverse Laplace transform of a function
$\hat{f}(s)$ is defined as

$\mathcal{L}^{-1}\left\\{\hat{f}(s)\right\\}=f(t)=\dfrac{1}{2\pi\mathrm{i}}\int_{x-\mathrm{i}\infty}^{x+\mathrm{i}\infty}\mathrm{e}^{ts}\hat{f}(s)\mathrm{d}s.$
(1i)

The Laplace transform of the $n-$th derivative of $f$ with respect to $t$is
given by the identity $\mathcal{L}\left\\{\partial^{n}f(t)/\partial
t^{n}\right\\}=s^{n}\hat{f}(s).$ The convolution between two functions $f_{1}$
and $f_{2}$ is defined as $\left(f_{1}(t)\ast
f_{2}(t)\right)=\int_{-\infty}^{\infty}f_{1}(u)f_{2}(t-\omega)\mathrm{d}\omega,$
and it follows that the Laplace transform convolution rule of two time
dependent functions $f_{1}$ and $f_{2}$ is given as
$\mathcal{L}\left\\{f_{1}(t)\ast
f_{2}(t)\right\\}=\mathcal{L}\left\\{f_{1}(t)\right\\}\mathcal{L}\left\\{f_{2}(t)\right\\}.$

Having established this, the governing equations on the transform space are

$\displaystyle\nabla\cdot\left(\textswab{C}^{m}\nabla\hat{{\mathbf{u}}}-{{\boldsymbol{\alpha}^{m}}\hat{\theta}}\right)+\hat{{\mathbf{b}}}=\rho
s^{2}\hat{{\mathbf{u}}}\,,$ (1ja)
$\displaystyle\hat{\mathbf{{q}}}+\tau^{m}s\hat{\mathbf{q}}=-(T_{0})^{-1}\bar{\mathbf{K}}^{m}\nabla\hat{\theta}\,,$
(1jb)
$\displaystyle\nabla\cdot\hat{\mathbf{{q}}}+{{\boldsymbol{\alpha}^{m}}}\nabla\hat{{\mathbf{u}}}+p^{m}s\hat{\theta}=\hat{r}\,,$
(1jc)

where $s$ is the unknown complex angular frequency $(s\in\mathbb{C}),$ and
$\hat{\mathbf{u}}$, $\hat{\theta},$ $\hat{{\mathbf{q}}},$
$\hat{{\mathbf{b}}},$ $\hat{r},$ are respectively, the bilateral Laplace
transformed of the micro-displacement field, the micro-temperature field, the
heat flux, the body forces and the heat source.

Proceeding with the substitutions, Eqs. $(\ref{gov1}),$ $(\ref{gov2})$ and
$(\ref{gov3})$ become on the Laplace transform space,

$\displaystyle\nabla\cdot\left({{\textswab{C}^{m}}\nabla\hat{{\mathbf{u}}}}-{{\boldsymbol{\alpha}^{m}}\hat{\theta}}\right)+\hat{{\mathbf{b}}}=\rho^{m}s^{2}\hat{{\mathbf{u}}}\,,$
(1ka)
$\displaystyle\nabla\cdot\left({{-\mathbf{K}^{m}}}\nabla\hat{\theta}\right)+{s{\boldsymbol{\alpha}^{m}}}\nabla\hat{{\mathbf{u}}}+p^{m}s\hat{{\theta}}=\hat{r}\,,$
(1kb)

where, due the plugging we assume

$\mathbf{K}^{m}\left({\mathbf{x}},\boldsymbol{\xi},s\right)=\dfrac{\bar{\mathbf{K}}^{m}(\mathbf{x},\boldsymbol{\xi})}{(1+\tau^{m}(\mathbf{x},\boldsymbol{\xi})s)T_{0}},$
(1l)

with a $\mathbf{K}^{m}$ well behaved over the defined domain.

The interface $\Sigma$ between two different phases $i$ and $j$ in the
periodic cell $\mathcal{A},$ holds the jump of the values of function $f$ at
it, written as $[[f]]=f^{i}(\Sigma)-f^{j}(\Sigma),$ follows that

$\displaystyle\left.\left[\left[\hat{\boldsymbol{u}}\right]\right]\right|_{\mathbf{x}\in\Sigma}=\mathbf{0},$
(1ma)
$\displaystyle\left.\left[\left[\left(\textswab{C}^{m}\nabla\hat{\boldsymbol{u}}-\boldsymbol{\alpha}^{m}\hat{\theta}\right)\cdot\mathbf{n}\right]\right]\right|_{\mathbf{x}\in\Sigma}=\mathbf{0},$
(1mb)
$\displaystyle\left.\left[\left[\hat{\theta}\right]\right]\right|_{\mathbf{x}\in\Sigma}=\mathbf{0},$
(1mc)
$\displaystyle\left.\left[\left[-\mathbf{K}^{m}\nabla\hat{\theta}\cdot\mathbf{n}\right]\right]\right|_{\mathbf{x}\in\Sigma}=\mathbf{0},$
(1md)

representing the conditions of continuity for a bounded interface, where the
outward normal to the interface $\Sigma$ is indicated by the vector
$\mathbf{n}.$ In case that $\tau^{m}$ is constant within the phases of the
material, and applying Laplace inverse transform we arrive in equations
$(\ref{eq:gov1})$ and $(\ref{eq:gov2})$ again.

On behalf of the time domain $t,$ let us apply the inverse Laplace transform
to the field Eqs. $(\ref{eq:Gov1T})$ and $(\ref{eq:Gov2T}),$ and due to the
dependency on the variable $s$ that $\mathbf{K}$ holds on equation
$(\ref{plugging}),$ which gives

$\mathcal{L}^{-1}\left\\{\mathbf{K}^{m}\nabla\hat{\theta}\right\\}=\mathcal{L}^{-1}\left\\{\mathbf{K}^{m}\right\\}\ast\nabla\theta,$
(1n)

the field equations $(\ref{eq:gov1})$ and $(\ref{eq:gov2}),$ might be
rewritten as an integral-differentiable form, therefore

$\displaystyle\nabla\cdot\left(\textswab{C}^{m}\nabla{\mathbf{u}}-{{\boldsymbol{\alpha}^{m}}\theta}\right)+{\mathbf{b}}=\rho^{m}{\mathbf{\ddot{u}}}\,,$
(1oa)
$\displaystyle\nabla\cdot\left(\mathcal{L}^{-1}\left\\{\mathbf{K}^{m}\right\\}\ast\nabla\theta-\boldsymbol{\alpha}^{m}\dot{\boldsymbol{u}}\right)+r=p^{m}\dot{\theta}\,$
(1ob)

known also as the field equations of a continuum with one time relaxation of
first order.

We must point out here that, from the $\mathcal{Q}$-periodicity of microscopic
constitutive tensors in Eqs. (1e), (1ka), (1kb), and bearing in mind the
$\mathcal{L}$-periodicity of the body forces, the micro-fields, displacement
and temperature (either both Laplace transformed or not) are dependent on the
fast variable $\boldsymbol{\xi}$ and the slow one $\mathbf{x},$ besides both
fields may be written as $\vskip 6.0pt plus 2.0pt minus
2.0pt\hat{\boldsymbol{u}}=\hat{\boldsymbol{u}}\left(\mathbf{x},\boldsymbol{\xi}=\frac{\mathbf{x}}{\varepsilon}\right)$
and
$\hat{\theta}=\hat{\theta}\left(\mathbf{x},\boldsymbol{\xi}=\frac{\mathbf{x}}{\varepsilon}\right).$
In addition, given the fact of the $\mathcal{Q}$-periodicity of their
coefficients, any attempt in deriving the solution of the set of PDE’s
$(\ref{eq:Gov1T})$ and $(\ref{eq:Gov2T})$ might be rather cumbersome both
analytically and numerically. Nevertheless, the homogenization technique
arrives in this context to barter the microstructured continuum by an
equivalent homogeneous one, where the solutions will be an approximation to
those from $(\ref{eq:Gov1T})$ and $(\ref{eq:Gov2T}),$ however on this scenario
the coefficients will not be affected by the fast variable $\boldsymbol{\xi},$
which gives rapid oscillations led by the underlying microstructure previously
held, and overall such procedure might be quite effective computationally.

Forward to the next Section, the homogenization technique on a bi-phase
periodic microstructure composite subjected to thermal sources and periodic
body forces will be applied over the equations of thermoelasticity with one
time relaxation on the Laplace transformed domain, namely equations
$(\ref{eq:Gov1T})$ and $(\ref{eq:Gov2T}).$

##  2 Asymptotic expansions of microscopic field equations on the Laplace
domain

Inspired by the asymptotic approach developed in [26, 28, 29, 39], the micro-
displacement $\boldsymbol{u}$ and micro-temperature $\theta$ are expressed as
power series in terms of $\varepsilon$, as well known as the asymptotic
expansions in terms of $\varepsilon$ that separates the fast variable
$\boldsymbol{\xi}=\textbf{x}/\varepsilon$ from the slow one $\textbf{x},$ i.e.
in the hypothesis of scale separation, diving into the components, follows
that

$\displaystyle
u_{h}\left(\textbf{x},\dfrac{\textbf{x}}{\varepsilon},t\right)=\sum_{l=0}^{+\infty}\varepsilon^{l}u^{(l)}_{h}=u^{(0)}_{h}\left(\textbf{x},\dfrac{\textbf{x}}{\varepsilon},t\right)+\varepsilon
u^{(1)}_{h}\left(\textbf{x},\dfrac{\textbf{x}}{\varepsilon},t\right)+\varepsilon^{2}u^{(2)}_{h}\left(\textbf{x},\dfrac{\textbf{x}}{\varepsilon},t\right)+\mathcal{O}(\varepsilon^{3}),$
(1pa)
$\displaystyle\theta\left(\textbf{x},\dfrac{\textbf{x}}{\varepsilon},t\right)=\sum_{l=0}^{+\infty}\varepsilon^{l}\theta^{(l)}=\theta^{(0)}\left(\textbf{x},\dfrac{\textbf{x}}{\varepsilon},t\right)+\varepsilon\theta^{(1)}\left(\textbf{x},\dfrac{\textbf{x}}{\varepsilon},t\right)+\varepsilon^{2}\theta^{(2)}\left(\textbf{x},\dfrac{\textbf{x}}{\varepsilon},t\right)+\mathcal{O}(\varepsilon^{3}).$
(1pb)

The Laplace transform $\eqref{eqn:lap}$ is applied to equations
$\eqref{eqn:u}$ and $\eqref{eqn:theta},$ which leads to

$\displaystyle\mathcal{L}\left(u_{h}\left(\textbf{x},\dfrac{\textbf{x}}{\varepsilon},t\right)\right)=\sum_{l=0}^{+\infty}\varepsilon^{l}\hat{u}^{(l)}_{h}=\hat{u}^{(0)}_{h}\left(\textbf{x},\dfrac{\textbf{x}}{\varepsilon},s\right)+\varepsilon\hat{u}^{(1)}_{h}\left(\textbf{x},\dfrac{\textbf{x}}{\varepsilon},s\right)+\varepsilon^{2}\hat{u}^{(2)}_{h}\left(\textbf{x},\dfrac{\textbf{x}}{\varepsilon},s\right)+\mathcal{O}(\varepsilon^{3}),$
(1qa)
$\displaystyle\mathcal{L}\left(\theta\left(\textbf{x},\dfrac{\textbf{x}}{\varepsilon},t\right)\right)=\sum_{l=0}^{+\infty}\varepsilon^{l}\hat{\theta}^{(l)}=\hat{\theta}^{(0)}\left(\textbf{x},\dfrac{\textbf{x}}{\varepsilon},s\right)+\varepsilon\hat{\theta}^{(1)}\left(\textbf{x},\dfrac{\textbf{x}}{\varepsilon},s\right)+\varepsilon^{2}\hat{\theta}^{(2)}\left(\textbf{x},\dfrac{\textbf{x}}{\varepsilon},s\right)+\mathcal{O}(\varepsilon^{3}).$
(1qb)

It is noteworthy that both asymptotic expansions above are equivalent to the
asymptotic expansions on the time domain $t,$ $\eqref{eqn:u}$ and
$\eqref{eqn:theta}.$

The homogenization procedure of the equations (1ka) and (1kb) searches for
solutions of the micro-displacement $\hat{u}_{h}$ and micro-temperature
$\hat{\theta},$ as decompositions of increasing powers of the microscopic
length $\varepsilon.$ In order to do so, the replacement of the asymptotic
expansions $\leavevmode\nobreak\ \eqref{eqr:ulaplace}$ and
$\leavevmode\nobreak\ \eqref{eqr:thetalaplace}$ must be done into the
microscopic field equations on the Laplace domain $\leavevmode\nobreak\
\eqref{eq:Gov1T}$ and $\eqref{eq:Gov2T},$ respectively. Further mathematical
details developing the asymptotic expansions in terms with equal power of
$\varepsilon$ of the field equations, among their interface conditions are
stressed in the Sec. A.1. of the Supplementary Material of the present work.

In fact, the arrived fields equations $(\ref{eqn:uhomogenized})$ and
$(\ref{eqn:thetahomogenized})$ in the Sec. A.1. of the Supplementary Material,
as power series obtained from the homogenization asymptotic procedure, have
the following structure,

$\displaystyle\dfrac{1}{\varepsilon^{2}}f_{i}^{(0)}(\textbf{x})+\dfrac{1}{\varepsilon}f_{i}^{(1)}(\textbf{x})+\varepsilon^{0}f_{i}^{(2)}(\textbf{x})+\varepsilon
f_{i}^{(3)}(\textbf{x})+\cdots+\varepsilon^{l}f_{i}^{(l+2)}(\textbf{x})+\hat{b}_{i}(\textbf{x})=0,$
(1ra)
$\displaystyle\dfrac{1}{\varepsilon^{2}}g^{(0)}(\textbf{x})+\dfrac{1}{\varepsilon}g^{(1)}(\textbf{x})+\varepsilon^{0}g^{(2)}(\textbf{x})+\varepsilon
g^{(3)}(\textbf{x})+\cdots+\varepsilon^{l}g^{(l+2)}(\textbf{x})+\hat{r}(\textbf{x})=0,$
(1rb)

where the set of functions $f_{i}^{(0)},\cdots,f_{i}^{(l+2)}$ and
$g^{(0)},\cdots,g^{(l+2)},$ with $l\in\mathbb{N},$ are such that the
dependency goes only over the slow variable $\textbf{x},$ and are determined
imposing the solvability conditions (seen in Sec. A.2.) on the class of the
$\mathcal{Q}$-periodic functions. Furthermore, the Sec. A.2. of the
Supplementary Material shows that the asymptotic field equations
$\eqref{eqn:uhomogenized}$ and $\eqref{eqn:thetahomogenized}$ produce a
cascade of recursive differential problems that determine sequentially the
structure of the solutions of the displacements $\boldsymbol{\hat{u}}^{(0)},$
$\boldsymbol{\hat{u}}^{(1)}$…, and also the structure of the solutions of the
temperature $\hat{\theta}^{(1)},$ $\hat{\theta}^{(1)}$…, respectively, in
terms of suitable $\mathcal{Q}$-periodic perturbation functions which are
determined by the solution of suitable cell problems detailed in the following
sub-session.

###  2.1 Cell problems and perturbation functions

The solutions of Eqs. $(\ref{eq:Gov1T})$ and $(\ref{eq:Gov2T})$ provide that
the inhomogeneous cell problems at the different orders of $\varepsilon$ can
be expanded in terms of the perturbation functions. In regards to them, such
perturbation functions depend exclusively on the microstructure features such
as, material geometry and mechanical properties, where the last one influences
the displacements and temperature due to the heterogeneity of the material.
Despite the existence of the homogenized solution $\mathcal{Q}$-periodic holds
without any other assumption [27, 28], it is noteworthy that imposing the
normalization conditions
$(\ref{eq:normalization1})-(\ref{eq:normalization3}),$ stated in the Sec. A.2.
of the Supplementary Material, to the cell problems, sooner stabilised on this
Section, the uniqueness of the perturbation functions is also held.

At this stage, the structure of these cell problems is exploited at the
different orders of $\varepsilon,$ for elasticity and thermal diffusion
problems. For this purpose, the recursive differential problems of the elastic
and thermal fields will be treated separately for each order of $\varepsilon.$

On behalf of the elastic problem at the order $\varepsilon^{-1},$ plugging the
Eq. $(\ref{eqn:solutiondisplac-1})$ , of the solution $u_{k}^{(1)},$ into the
differential problem $(\ref{eq:epsilon-1f}),$ both equation from the Sec. A.2.
of the Supplementary Material, lead to two cell problems:

$\displaystyle\vskip 6.0pt plus 2.0pt minus
2.0pt\left(C_{ijhk}^{m}\,N_{hpq_{1},k}^{(1,0)}\right)_{,j}+C_{ijpq_{1},j}^{m}=0,$
(1sa)
$\displaystyle\left(C_{ijhk}^{m}\tilde{N}_{h,k}^{(1,0)}\right)_{,j}-\alpha_{ij,j}^{m}=0.$
(1sb)

Similarly, from Eq. $(\ref{eq:epsilon-1g})$ and in consideration of the
solution $(\ref{eqn:solutiontemp-1})$ reported in the Sec. A.2. of the
Supplementary Material, after substitution, the cell problem at the order
$\varepsilon^{-1}$ takes the form

$\displaystyle\left(K_{ij}^{m}M_{q_{1},j}^{(1,0)}\right)_{,i}+K_{iq_{1},i}^{m}=0.$
(1t)

Therefore, for the terms at the order $\varepsilon^{-1},$ we obtain the
following set of cell problems

$\displaystyle\left(C_{ijhk}^{m}\,N_{hpq_{1},k}^{(1,0)}\right)_{,j}+C_{ijpq_{1},j}^{m}=0,$
(1ua)
$\displaystyle\left(C_{ijhk}^{m}\tilde{N}_{h,k}^{(1,0)}\right)_{,j}-\alpha_{ij,j}^{m}=0,$
(1ub)
$\displaystyle\left(K_{ij}^{m}M_{q_{1},j}^{(1,0)}\right)_{,i}+K_{iq_{1},i}^{m}=0,$
(1uc)

with interface conditions expressed in terms of perturbation functions
$N_{hpq_{1}}^{(1,0)},$ $\tilde{N}_{h}^{(1,0)},$ $M_{q_{1}}^{(1,0)}$ as

$\displaystyle\left.\left[\left[N_{hpq_{1}}^{(1,0)}\right]\right]\right|_{\boldsymbol{\xi}\in\Sigma_{1}}=0,$
(1va)
$\displaystyle\left.\left[\left[\left(C_{ijhk}^{m}\left(\delta_{hp}\,\delta_{kq_{1}}+N_{hpq_{1},k}^{(1,0)}\right)\right)n_{j}\right]\right]\right|_{\boldsymbol{\xi}\in\Sigma_{1}}=0,$
(1vb)
$\displaystyle\left.\left[\left[\tilde{N}_{h}^{(1,0)}\right]\right]\right|_{\boldsymbol{\xi}\in\Sigma_{1}}=0,$
(1vc)
$\displaystyle\left.\left[\left[\left(C_{ijhk}^{m}\left(\tilde{N}_{h,k}^{(1,0)}-\alpha^{m}_{ij}\right)\right)n_{j}\right]\right]\right|_{\boldsymbol{\xi}\in\Sigma_{1}}=0,$
(1vd)
$\displaystyle\left.\left[\left[M_{q_{1}}^{(1,0)}\right]\right]\right|_{\boldsymbol{\xi}\in\Sigma_{1}}=0,$
(1ve)
$\displaystyle\left.\left[\left[\left(K_{ij}^{m}\left(M_{q_{1},j}^{(1,0)}+\delta_{jq_{1}}\right)\right)n_{i}\right]\right]\right|_{\boldsymbol{\xi}\in\Sigma_{1}}=0,$
(1vf)

where the symmetries from the respective tensors are passed to the
perturbation functions, for instance
$N_{hpq_{1}}^{(1,0)}=N_{hq_{1}p}^{(1,0)},$ since
$C_{ijpq_{1},j}^{m}=C_{ijq_{1}p,j}^{m}.$

Once the perturbation functions $N_{hpq_{1}}^{(1,0)},$ $\tilde{N}_{h}^{(1,0)}$
and $M_{q_{1}}^{(1)}$ are determined, from the differential problem
$(\ref{eq:epsilon0f})$ and recalling the solution $u_{k}^{(2)}$ in Eq.
$(\ref{eqn:solutiondisplac0})$ stated in the Sec. A.2. of the Supplementary
Material, one derives the three following cell problems, at the order
$\varepsilon^{0},$ and symmetrizing with respect to indices $q_{1}$ and
$q_{2},$ leads to

$\displaystyle\begin{split}&\vskip 6.0pt plus 2.0pt minus
2.0pt\left(C_{ijhk}^{m}\,N_{hpq_{1}q_{2},k}^{(2,0)}\right)_{,j}+\dfrac{1}{2}\left[\left(C_{ijhq_{2}}^{m}\,N_{hpq_{1}}^{(1,0)}\right)_{,j}+C_{iq_{1}pq_{2}}^{m}+C_{iq_{2}hj}^{m}\,N_{hpq_{1},j}^{(1,0)}+\left(C_{ijhq_{1}}^{m}\,N_{hpq_{2}}^{(1,0)}\right)_{,j}+\right.\\\
&\vskip 6.0pt plus 2.0pt minus
2.0pt\left.+\,C_{iq_{2}pq_{1}}^{m}+C_{iq_{1}hj}^{m}\,N_{hpq_{2},j}^{(1,0)}\right]=\dfrac{1}{2}\left\langle
C_{iq_{1}pq_{2}}^{m}+C_{iq_{2}hj}^{m}\,N_{hpq_{1},j}^{(1,0)}+C_{iq_{2}pq_{1}}^{m}+C_{iq_{1}hj}^{m}\,N_{hpq_{2},j}^{(1,0)}\right\rangle,\end{split}$
(1wa) $\displaystyle\begin{split}&\vskip 6.0pt plus 2.0pt minus
2.0pt\left(C_{ijhk}^{m}\,\tilde{N}_{hq_{1},k}^{(2,1)}\right)_{,j}+\left[\left(C_{ijhq_{1}}^{m}\,\tilde{N}_{h}^{(1,0)}\right)_{,j}+C_{iq_{1}hj}^{m}\,\tilde{N}_{h,j}^{(1,0)}-\left(\alpha_{ij}^{m}\,M_{q_{1}}^{(1,0)}\right)_{,j}-\alpha_{iq_{1}}^{m}\right]=\left\langle
C_{iq_{1}hj}^{m}\,\tilde{N}_{h,j}^{(1,0)}-\alpha_{iq_{1}}^{m}\right\rangle\\\
&\left(C_{ijhk}^{m}\,N_{hp,k}^{(2,2)}\right)_{,j}-\rho^{m}\delta_{ip}=-\left\langle\rho^{m}\right\rangle\delta_{ip},\end{split}$
(1wb)

with interface conditions in terms of perturbation functions as

$\displaystyle\left.\left[\left[N_{hpq_{1}q_{2}}^{(2,0)}\right]\right]\right|_{\boldsymbol{\xi}\in\Sigma_{1}}=0,$
(1xa)
$\displaystyle\left.\left[\left[\left[C_{ijhk}^{m}N_{hpq_{1}q_{2},k}^{(2,0)}+\frac{1}{2}\left(C_{ijhq_{2}}^{m}N_{hpq_{1}}^{(1,0)}+C_{ijhq_{1}}^{m}N_{hpq_{2}}^{(1,0)}\right)\right]n_{j}\right]\right]\right|_{\boldsymbol{\xi}\in\Sigma_{1}}=0,$
(1xb)
$\displaystyle\left.\left[\left[\tilde{N}_{hq_{1}}^{(2,1)}\right]\right]\right|_{\boldsymbol{\xi}\in\Sigma_{1}}=0,$
(1xc)
$\displaystyle\left.\left[\left[\left(C_{ijhk}^{m}\tilde{N}_{hq_{1},k}^{(2,1)}+\left(C_{ijhq_{2}}^{m}\,\tilde{N}_{h}^{(1,0)}\delta_{kq_{1}}-\alpha_{ij}^{m}M_{q_{1}}^{(1,0)}\right)\right)n_{j}\right]\right]\right|_{\boldsymbol{\xi}\in\Sigma_{1}}=0,$
(1xd)
$\displaystyle\left.\left[\left[N^{(2,2)}_{hp}\right]\right]\right|_{\boldsymbol{\xi}\in\Sigma_{1}}=0,$
(1xe)
$\displaystyle\left.\left[\left[\left(C_{ijhk}^{m}N_{hp,k}^{(2,2)}\right)n_{j}\right]\right]\right|_{\boldsymbol{\xi}\in\Sigma_{1}}=0.$
(1xf)

Analogously, for the thermal homogenized field at the order $\varepsilon^{0},$
replacing solution $(\ref{eqn:solutiontemp0})$ into the differential problem
$(\ref{eq:epsilon0g})$ derived in the Sec. A.2. of the Supplementary Material,
one gives the next three cell problems

$\displaystyle\begin{split}&\vskip 6.0pt plus 2.0pt minus
2.0pt\left(K_{ij}^{m}\,M_{q_{1}q_{2},j}^{(2,0)}\right)_{,i}+\dfrac{1}{2}\left[\left(K_{iq_{2}}^{m}\,M_{q_{1}}^{(1,0)}\right)_{,i}+K_{q_{1}q_{2}}^{m}+K_{iq_{2}}^{m}\,M_{q_{1},i}^{(1,0)}+\left(K_{iq_{1}}^{m}\,M_{q_{2}}^{(1,0)}\right)_{,i}+K_{q_{2}q_{1}}^{m}+\right.\\\
&\vskip 6.0pt plus 2.0pt minus
2.0pt\left.K_{q_{1}i}^{m}\,M_{q_{2},i}^{(1,0)}\right]=\dfrac{1}{2}\left\langle
K_{q_{1}q_{2}}^{m}+K_{q_{2}i}^{m}\,M_{q_{1},i}^{(1,0)}+K_{q_{2}q_{1}}^{m}+K_{q_{1}i}^{m}\,M_{q_{2},i}^{(1,0)}\right\rangle,\end{split}$
(1ya) $\displaystyle\vskip 6.0pt plus 2.0pt minus
2.0pt\left(K_{ij}^{m}\,\tilde{M}_{pq_{1},j}^{(2,1)}\right)_{,i}-\left[\alpha_{ij}^{m}\,N_{ipq_{1},j}^{(1,0)}+\alpha_{pq_{1}}^{m}\,\right]=-\left\langle\alpha_{ij}^{m}\,N_{ipq_{1},j}^{(1,0)}+\alpha_{pq_{1}}^{m}\right\rangle$
(1yb) $\displaystyle\vskip 6.0pt plus 2.0pt minus
2.0pt\left(K_{ij}^{m}\,M_{,j}^{(2,1)}\right)_{,i}-\left[\alpha_{ij}^{m}\,\tilde{N}_{i,j}^{(1,0)}+p^{m}\right]=-\left\langle\alpha_{ij}^{m}\,\tilde{N}_{i,j}^{(1,0)}+p^{m}\right\rangle,$
(1yc)

and, in terms of perturbation functions, the interface conditions become

$\displaystyle\left.\left[\left[M_{q_{1}q_{2}}^{(2,0)}\right]\right]\right|_{\boldsymbol{\xi}\in\Sigma_{1}}=0,$
(1za)
$\displaystyle\left.\left[\left[\left[K_{ij}^{m}M_{q_{1}q_{2},j}^{(2,0)}+\frac{1}{2}\left(K_{iq_{2}}^{m}\,M_{q_{1}}^{(1,0)}+K_{iq_{1}}^{m}\,M_{q_{2}}^{(1,0)}\right)\right]n_{i}\right]\right]\right|_{\boldsymbol{\xi}\in\Sigma_{1}}=0.$
(1zb)
$\displaystyle\left.\left[\left[\tilde{M}_{pq_{1}}^{(2,1)}\right]\right]\right|_{\boldsymbol{\xi}\in\Sigma_{1}}=0,$
(1zc)
$\displaystyle\left.\left[\left[\left(K_{ij}^{m}(s)\tilde{M}_{pq_{1},j}^{(2,1)}\right)n_{i}\right]\right]\right|_{\boldsymbol{\xi}\in\Sigma_{1}}=0.$
(1zd)
$\displaystyle\left.\left[\left[M^{(2,1)}\right]\right]\right|_{\boldsymbol{\xi}\in\Sigma_{1}}=0,$
(1ze)
$\displaystyle\left.\left[\left[\left(K_{ij}^{m}(s)M_{,j}^{(2,1)}\right)n_{i}\right]\right]\right|_{\boldsymbol{\xi}\in\Sigma_{1}}=0.$
(1zf)

###  2.2 Down-scaling and up-scaling relations

On one hand, the feasibility of expressing the microscopic fields
$\hat{u}_{h}(\mathbf{x},\boldsymbol{\xi},s)$ and
$\hat{\theta}(\mathbf{x},\boldsymbol{\xi},s)$ as asymptotic expansions of
powers of the microscopic length $\varepsilon$ in terms of the macroscopic
fields $\hat{U}^{M}_{h}(\mathbf{x},s)$, $\hat{\Theta}^{M}(\mathbf{x},s),$
their gradients and in terms of the $\mathcal{Q}$-periodic perturbation
functions arises from the solution of the cell problems at different orders of
$\varepsilon$ developed in Sec. 2.1. From the expansions $(\ref{eqn:theta})$
and $(\ref{eqn:u}),$ and considering the solutions brought in the Sec. A.2. of
the Supplementary Material $\leavevmode\nobreak\
\eqref{eqn:solutiondisplac-2},$ $\leavevmode\nobreak\
\eqref{eqn:solutiontemp-2},$ $\leavevmode\nobreak\
\eqref{eqn:solutiondisplac-1},$ $\leavevmode\nobreak\
\eqref{eqn:solutiontemp-1},$ $\leavevmode\nobreak\
\eqref{eqn:solutiondisplac0}$ and $\leavevmode\nobreak\
\eqref{eqn:solutiontemp0}$ of cell problems at the different orders of
$\varepsilon,$ one derives the down-scaling relations of the transformed
micro-displacement field and the transformed micro-temperature field,
respectively,

$\displaystyle\begin{split}&\hat{u}_{h}\left(\mathbf{x},\frac{\mathbf{x}}{\varepsilon},s\right)=\left[\hat{U}^{M}_{h}(\mathbf{x},s)+\varepsilon\left(N^{(1,0)}_{hpq_{1}}(\boldsymbol{\xi})\dfrac{\partial\hat{U}^{M}_{h}}{\partial
x_{q_{1}}}+\tilde{N}^{(1,0)}_{h}(\boldsymbol{\xi})\hat{\Theta}^{M}\right)\right.+\\\
&+\left.\left.\varepsilon^{2}\left(N^{(2,0)}_{hpq_{1}q_{2}}(\boldsymbol{\xi})\dfrac{\partial^{2}\hat{U}^{M}_{h}}{\partial
x_{q1}\partial
x_{q2}}+\tilde{N}^{(2,1)}_{hq_{1}}(\boldsymbol{\xi})\dfrac{\partial\hat{\Theta}^{M}}{\partial
x_{q_{1}}}+s^{2}N^{(2,2)}_{hp}(\boldsymbol{\xi})\hat{U}^{M}_{h}\right)+\mathcal{O}(\varepsilon^{3})\right]\right|_{\boldsymbol{\xi}=\frac{\mathbf{x}}{\varepsilon}},\end{split}$
(1aaa)
$\displaystyle\begin{split}&\hat{\theta}\left(\mathbf{x},\frac{\mathbf{x}}{\varepsilon},s\right)=\left[\hat{\Theta}^{M}(\mathbf{x},s)+\varepsilon\left(M^{(1,0)}_{q_{1}}(\boldsymbol{\xi})\dfrac{\partial\hat{\Theta}^{M}}{\partial
x_{q_{1}}}\right)+\right.\\\
&+\left.\left.\varepsilon^{2}\left(M^{(2,0)}_{pq_{1}q_{2}}(\boldsymbol{\xi})\dfrac{\partial^{2}\hat{\Theta}^{M}}{\partial
x_{q1}\partial
x_{q2}}+s\tilde{M}^{(2,1)}_{pq_{1}}(\boldsymbol{\xi})\dfrac{\partial\hat{\Theta}^{M}}{\partial
x_{q_{1}}}+s^{2}M^{(2,1)}_{hp}(\boldsymbol{\xi})\hat{\Theta}^{M}\right)+\mathcal{O}(\varepsilon^{3})\right]\right|_{\boldsymbol{\xi}=\frac{\mathbf{x}}{\varepsilon}}.\end{split}$
(1aab)

Noteworthy, the $\mathcal{Q}$-periodic perturbation functions
$N_{hpq_{1}}^{(1,0)},\tilde{N}_{h}^{(1,0)},N_{hpq_{1}q_{2}}^{(2,0)},\tilde{N}^{(2,0)}_{hq_{1}},N_{hp}^{(2,2)},M^{(1,0)}_{q_{1}},$
$M^{(2,0)}_{q_{1}q_{2}},$ $\tilde{M}^{(2,1)}_{pq_{1}},$ $M^{(2,1)},$ cast the
effect of microstructural inhomogeneities of the material through their
dependency on the fast variable $\boldsymbol{\xi}=\mathbf{x}/\varepsilon$,
whereas the macro-fields $\hat{\boldsymbol{u}}(\mathbf{x},s),$ and
$\hat{\theta}(\mathbf{x},s)$ are $\mathfrak{L}$-periodic functions and
therefore rely only on the slow variable $\mathbf{x}$ [62].

On the other hand, also known as the up-scaling relations, they allow the
definition of the macroscopic fields in terms of the microscopic ones. This
means that the macroscopic fields may be defined as the mean values of
microscopic equations $(\ref{downscalingu})$ and $(\ref{downscalingtheta})$
over the unit cell $\mathcal{Q},$ thus it follows that

$\displaystyle\begin{split}\hat{U}_{h}(\mathbf{x},s)\doteq\left\langle\hat{u}_{h}\left(\mathbf{x},\dfrac{\mathbf{x}}{\varepsilon}+{\boldsymbol{\zeta}},s\right)\right\rangle=&\dfrac{1}{|\mathcal{Q}|}\int_{\mathcal{Q}}\hat{u}_{h}\left(\mathbf{x},\dfrac{\mathbf{x}}{\varepsilon}+{\boldsymbol{\zeta}},s\right)d\boldsymbol{\zeta}=\\\
=&\dfrac{1}{|\mathcal{Q}|}\int_{\mathcal{Q}}\hat{u}_{h}\left(\mathbf{x},\boldsymbol{\xi},s\right)d\boldsymbol{\xi}=\left\langle\hat{u}_{h}\left(\mathbf{x},{\boldsymbol{\xi}},s\right)\right\rangle,\end{split}$
(1aba)
$\displaystyle\begin{split}\hat{\Theta}(\mathbf{x},s)\doteq\left\langle\hat{\theta}\left(\mathbf{x},\dfrac{\mathbf{x}}{\varepsilon}+\boldsymbol{\zeta},s\right)\right\rangle=&\dfrac{1}{|\mathcal{Q}|}\int_{\mathcal{Q}}\hat{\theta}\left(\mathbf{x},\dfrac{\mathbf{x}}{\varepsilon}+{\boldsymbol{\zeta}},s\right)d\boldsymbol{\zeta}=\\\
=&\dfrac{1}{|\mathcal{Q}|}\int_{\mathcal{Q}}\hat{\theta}\left(\mathbf{x},\boldsymbol{\xi},s\right)d\boldsymbol{\xi}=\left\langle\hat{\theta}\left(\mathbf{x},\boldsymbol{\xi},s\right)\right\rangle,\quad\end{split}$
(1abb)

where the variable $\boldsymbol{\zeta}\in\mathcal{Q}$ has been introduced and
the vector $\varepsilon\boldsymbol{\zeta}\in\mathcal{A}$ defines all the
possible translations of the heterogeneous medium compared to a grid of cells
having characteristic size $\varepsilon,$ with respect to the
$\mathfrak{L}$-periodic body force $\hat{\mathbf{b}}(\mathbf{x})$ [29, 78].
Besides, the mean value operator of the left hand side is taken over the
variable $\boldsymbol{\zeta},$ whereas the mean value operator of the right
hand side has become over $\boldsymbol{\xi}.$

##  3 Approximation of the power-like functional through truncation of its
asymptotic expansion

A variational-asymptotic procedure [29, 39] is exploited to furnish a finite
order governing equation to establish an equivalence among the equations (1ka)
and (1kb) at the macro-scale. For this matter, on the periodic domain
$\mathfrak{L},$ let the power-like functional ${\Lambda}$ be expressed in
terms of the energy-like density $\lambda_{m}$ at the micro-scale

$\displaystyle{\Lambda}(\boldsymbol{u},\theta)$
$\displaystyle=\dfrac{\partial}{\partial
t}\int_{\mathfrak{L}}\lambda_{m}\left(\boldsymbol{x},\frac{\boldsymbol{x}}{\varepsilon}\right)d\mathbf{x}=$
$\displaystyle=\dfrac{\partial}{\partial
t}\int_{\mathfrak{L}}\left(\dfrac{1}{2}\rho^{m}\dot{\boldsymbol{u}}\ast\dot{\boldsymbol{u}}+\dfrac{1}{2}\nabla\boldsymbol{u}\ast(\mathfrak{C}^{m}\nabla\boldsymbol{u})-\dfrac{1}{2}\nabla\boldsymbol{u}\ast(\boldsymbol{\alpha}^{m}\theta)-\boldsymbol{u}\ast\boldsymbol{b}\right)d\mathbf{x}\quad+$
(1ac)
$\displaystyle-\int_{\mathfrak{L}}\left(\dfrac{1}{2}\nabla\theta\ast\left(\mathcal{L}^{-1}\left\\{\mathbf{K}^{m}\right\\}\ast\nabla\theta\right)+\dfrac{1}{2}\theta\ast\left(\boldsymbol{\alpha}^{m}\nabla\dot{\boldsymbol{u}}\right)+\dfrac{1}{2}\theta\ast(p^{m}\dot{\theta})-\theta\ast
r\right)d\mathbf{x}\,.$

Now, applying the Laplace transform on $\mathcal{L}({\Lambda}),$ the
functional above is taken to the Laplace domain and expressed in terms of the
power-like density $\hat{\lambda}_{m}$ as follows

$\displaystyle\mathcal{L}\\{{\Lambda}\\}$
$\displaystyle=\hat{\Lambda}(\hat{\boldsymbol{u}},\hat{\theta})=\int_{\mathfrak{L}}\hat{\lambda}_{m}\left(\boldsymbol{x},\frac{\boldsymbol{x}}{\varepsilon}\right)d\mathbf{x}=$
$\displaystyle=\int_{\mathfrak{L}}s\left(\dfrac{1}{2}\rho^{m}s^{2}\hat{\boldsymbol{u}}\cdot\hat{\boldsymbol{u}}+\dfrac{1}{2}\nabla\hat{\boldsymbol{u}}:(\mathfrak{C}^{m}\nabla\hat{\boldsymbol{u}})-\dfrac{1}{2}\nabla\hat{\boldsymbol{u}}:(\boldsymbol{\alpha}^{m}\hat{\theta})-\hat{\boldsymbol{u}}\cdot\hat{\boldsymbol{b}}\right)d\mathbf{x}\quad+$
(1ad)
$\displaystyle-\int_{\mathfrak{L}}\left(\dfrac{1}{2}\nabla\hat{\theta}\cdot\left(\mathbf{K}^{m}\nabla\hat{\theta}\right)+\dfrac{1}{2}s\hat{\theta}\left(\boldsymbol{\alpha}^{m}:\nabla\hat{\boldsymbol{u}}\right)+\dfrac{1}{2}s\hat{\theta}(p^{m}\hat{\theta})-\hat{\theta}\hat{r}\right)d\mathbf{x}\,,$

where the symbol $:$ indicates the tensor double inner product. Note that the
Euler-Lagrange equation from the functional (3) agrees with the demonstration
presented in the Sec. A.5. which leads to the field equations at the micro-
scale $(\ref{eq:Gov1T})$ and $(\ref{eq:Gov2T}).$

As it was introduced at the equations (1aba) and (1abb), the transformed
power-like functional $\hat{\Lambda}$ and its corresponding power-like density
$\hat{\lambda}_{m}$ are dependent by the translation parameter
$\boldsymbol{\zeta}\in\mathcal{Q}.$ Which means, the perturbation functions
$N_{hpq_{1}}^{(1,0)},\tilde{N}_{h}^{(1,0)},N_{hpq_{1}q_{2}}^{(2,0)},\tilde{N}^{(2,0)}_{hq_{1}},N_{hp}^{(2,2)},M^{(1,0)}_{q_{1}},$
$M^{(2,0)}_{q_{1}q_{2}},$ $\tilde{M}^{(2,1)}_{pq_{1}},$ $M^{(2,1)},$
determined from the cell problems in Sec. 2.1, also depend on the translation
variable $\boldsymbol{\zeta}.$ Besides, the power-like density
$\hat{\lambda}_{m}$ in the Laplace domain complies with the property

$\hat{\lambda}^{\zeta}_{m}\left(\textbf{x},\dfrac{\textbf{x}}{\varepsilon}\right)=\hat{\lambda}_{m}\left(\textbf{x},\frac{\textbf{x}}{\varepsilon}+\boldsymbol{\zeta}\right),$
(1ae)

implying that the Laplace transformed of the power-like functional
$\hat{\Lambda},$ depends on the parameter $\boldsymbol{\zeta},$ one has

$\displaystyle\hat{\Lambda}^{\zeta}=\hat{\Lambda}(\boldsymbol{\zeta})=\int_{\mathfrak{L}}\hat{\lambda}^{\zeta}_{m}\left(\textbf{x},\dfrac{\textbf{x}}{\varepsilon}\right)d\textbf{x}=\int_{\mathfrak{L}}\hat{\lambda}_{m}\left(\textbf{x},\frac{\textbf{x}}{\varepsilon}+\boldsymbol{\zeta}\right)d\textbf{x}.$
(1af)

Let $\hat{\Lambda}_{m}$ be the average transformed power-like functional at
the micro-scale

$\displaystyle\begin{split}&\hat{\Lambda}_{m}\dot{=}\langle\hat{\Lambda}^{\zeta}\rangle=\dfrac{1}{|\mathcal{Q}|}\int_{\mathcal{Q}}\left[\int_{\mathfrak{L}}\hat{\lambda}_{m}\left(\textbf{x},\frac{\textbf{x}}{\varepsilon}+\boldsymbol{\zeta}\right)d\textbf{x}\right]d\boldsymbol{\zeta}=\\\
&=\int_{\mathfrak{L}}\left[\dfrac{1}{|\mathcal{Q}|}\int_{\mathcal{Q}}\hat{\Lambda}(\boldsymbol{\zeta})d\boldsymbol{\zeta}\right]d\textbf{x}=\int_{\mathfrak{L}}\left\langle\hat{\lambda}_{m}\left(\textbf{x},\dfrac{\textbf{x}}{\varepsilon}\right)\right\rangle
d\textbf{x},\end{split}$ (1ag)

where the Fubini’s theorem was applied in the third equality.

Since the power-like functional $\hat{\Lambda}^{\zeta}$ is averaged by the
translated realizations of the microstructure, which entails that the
transformed power-like density at the micro-scale satisfies

$\displaystyle\Big{\langle}\hat{\lambda}_{m}\Big{(}\boldsymbol{x},\frac{\boldsymbol{x}}{\varepsilon}+\boldsymbol{\zeta}\Big{)}\Big{\rangle}=\frac{1}{|\mathcal{Q}|}\int_{\mathcal{Q}}\hat{\lambda}_{m}\Big{(}\boldsymbol{x},\frac{\boldsymbol{x}}{\varepsilon}+\boldsymbol{\zeta}\Big{)}d\boldsymbol{\zeta}=\frac{1}{|\mathcal{Q}|}\int_{\mathcal{Q}}\hat{\lambda}_{m}\Big{(}\boldsymbol{x},\boldsymbol{\xi}\Big{)}d\boldsymbol{\xi}=\langle\hat{\lambda}_{m}(\boldsymbol{x},\boldsymbol{\xi})\rangle,$
(1ah)

consequently, this ensures that the average transformed power-like functional
$\langle\hat{\Lambda}^{\zeta}\rangle$ at the micro-scale does not rely on the
translation variable $\boldsymbol{\zeta}.$

Through a similar variational approach made in the Sec. A.5., we will
determine the governing field equations at the macro-scale and the overall
constitutive tensors, finding the Euler-Lagrange equation of the power-like
function from the down-scaling relations (1aaa) and (1aab). So, let us once
again introduce the down-scaling relations related to the transformed micro-
displacement $\boldsymbol{\hat{u}}(\textbf{x},\boldsymbol{\xi},s),$ and to the
transformed micro-temperature $\hat{\theta}(\textbf{x},\boldsymbol{\xi},s)$,
i.e.

$\displaystyle\hat{u}_{h}\left(\mathbf{x},\frac{\mathbf{x}}{\varepsilon},s\right)=\hat{U}^{M}_{h}(\mathbf{x},s)+\varepsilon\left(N^{(1,0)}_{hpq_{1}}(\boldsymbol{\xi})\dfrac{\partial\hat{U}^{M}_{h}}{\partial
x_{q_{1}}}+\tilde{N}^{(1,0)}_{h}(\boldsymbol{\xi})\hat{\Theta}^{M}\right)+\mathcal{O}(\varepsilon^{2}),$
(1aia)
$\displaystyle\hat{\theta}\left(\mathbf{x},\frac{\mathbf{x}}{\varepsilon},s\right)=\hat{\Theta}^{M}(\mathbf{x},s)+\varepsilon\left(M^{(1,0)}_{q_{1}}(\boldsymbol{\xi})\dfrac{\partial\hat{\Theta}^{M}}{\partial
x_{q_{1}}}\right)+\mathcal{O}(\varepsilon^{2}),$ (1aib)

and applying the gradients (1bza) and (1bzb) to the down-scaling relations
above, one gives

$\displaystyle\dfrac{D\hat{u}_{h}}{Dx_{k}}$
$\displaystyle=\dfrac{\partial\hat{U}^{M}_{h}}{\partial
x_{k}}+N^{(1,0)}_{hpq_{1},k}\dfrac{\partial\hat{U}^{M}_{p}}{\partial
x_{q_{1}}}+\tilde{N}^{(1,0)}_{h,k}\hat{\Theta}^{M}+\varepsilon\left(N^{(1,0)}_{hpq_{1}}\dfrac{\partial^{2}\hat{U}^{M}_{p}}{\partial
x_{q_{1}}\partial
x_{k}}+\tilde{N}^{(1,0)}_{h}\dfrac{\partial\hat{\Theta}^{M}}{\partial
x_{k}}\right)+\mathcal{O}(\varepsilon^{2})=$
$\displaystyle=B^{(1,0)}_{hkpq_{1}}\dfrac{\partial\hat{U}^{M}_{h}}{\partial
x_{q_{1}}}+\tilde{B}^{(1,0)}_{hk}\hat{\Theta}^{M}+\mathcal{O}(\varepsilon),$
(1aj)

where $B^{(1,0)}_{hkpq_{1}}=\delta_{hp}\delta_{kq_{1}}+N^{(1,0)}_{hpq_{1},k}$
and $\tilde{B}^{(1,0)}_{hk}=\tilde{N}^{(1,0)}_{h,k},$ and

$\displaystyle\dfrac{D\hat{\theta}_{h}}{Dx_{j}}$
$\displaystyle=\dfrac{\partial\hat{\Theta}^{M}}{\partial
x_{j}}+M^{(1,0)}_{q_{1},j}\dfrac{\partial\hat{\Theta}^{M}}{\partial
x_{q_{1}}}+\varepsilon
M^{(1,0)}_{q_{1}}\dfrac{\partial^{2}\hat{\Theta}^{M}}{\partial
x_{q_{1}}\partial x_{j}}+\mathcal{O}(\varepsilon^{2})=$
$\displaystyle=A^{(1,0)}_{jq_{1}}\dfrac{\partial\hat{\Theta}^{M}}{\partial
x_{q_{1}}}+\mathcal{O}(\varepsilon),$ (1ak)

where $A^{(1,0)}_{jq_{1}}=\delta_{jq_{1}}+M^{(1,0)}_{q_{1},j}.$ The tensors
$B^{(1,0)}_{hkpq_{1}},$ $\tilde{B}^{(1,0)}_{hk}$ and $A^{(1,0)}_{jq_{1}}$ are
called localization tensors and are also $\mathcal{Q}$-periodic functions in
relation to the fast variable $\boldsymbol{\xi},$ once the perturbation
functions and their gradients are $\mathcal{Q}$-periodic functions.

As previously mentioned, such first variation must vanish for all admissible
$\delta\hat{U}^{M}_{h},$ $\delta\hat{\Theta}^{M}$ i.e,
$\delta\hat{\Lambda}_{m}(\hat{U}^{M}_{h},\delta\hat{U}^{M}_{h})=0,$ and hence
the Euler-Lagrange differential equation associated to the variational problem
(1dd) leads to the following two equations

$\displaystyle s^{2}\langle\rho^{m}\rangle\hat{U}^{M}_{l}-\left\langle
B^{(1,0)}_{ijlr_{1}}C^{m}_{ijhk}B^{(1,0)}_{hkpq_{1}}\right\rangle\dfrac{\partial^{2}\hat{U}^{M}_{p}}{\partial
x_{q_{1}}\partial x_{r_{1}}}+\left\langle
B^{(1,0)}_{ijlq_{1}}\alpha^{m}_{ij}-B^{(1,0)}_{hklq_{1}}C^{m}_{ijhk}\tilde{B}^{(1,0)}_{ij}\right\rangle\dfrac{\partial\hat{\Theta}^{M}}{\partial
x_{q_{1}}}-\hat{b}_{l}=0,$ (1ala) $\displaystyle\begin{split}&\left\langle
A^{(1,0)}_{ir_{1}}K^{m}_{ij}A^{(1,0)}_{jq_{1}}\right\rangle\dfrac{\partial^{2}\hat{\Theta}^{M}}{\partial
x_{r_{1}}x_{q_{1}}}+s\left(\left\langle
B^{(1,0)}_{hkpq_{1}}C^{m}_{ijhk}\tilde{B}^{(1,0)}_{ij}\right\rangle\dfrac{\partial\hat{U}^{M}_{l}}{\partial
x_{q_{1}}}-\left\langle
B^{(1,0)}_{ijpq_{1}}\alpha^{m}_{ij}\right\rangle\dfrac{\partial\hat{U}^{M}_{l}}{\partial
x_{q_{1}}}\right)\quad+\\\ &-s\langle
p^{m}+2\tilde{B}^{(1,0)}_{ij}\alpha^{m}_{ij}+\tilde{B}^{(1,0)}_{ij}C^{m}_{ijhk}\tilde{B}^{(1,0)}_{hk}\rangle\hat{\Theta}^{M}+\hat{r}=0,\end{split}$
(1alb)

which are clearly developed in terms of the Laplace transformed macro-
temperature and macro-displacement among with their gradients. Since the
equations (1ala) and (LABEL:eulerlagrange2) are set by positive definite
constitutive tensors, one implies that the existence and the uniqueness of
their solutions are endorsed by the Legendre-Hadamard condition.

Rearranging and rewriting the equations (1ala) and (LABEL:eulerlagrange2) in
terms of the constitutive tensors, one gives

$\displaystyle C_{lr_{1}pq_{1}}\dfrac{\partial^{2}\hat{U}^{M}_{p}}{\partial
x_{q_{1}}\partial
x_{r_{1}}}-\alpha_{lq_{1}}\dfrac{\partial\hat{\Theta}^{M}}{\partial
x_{q_{1}}}+\hat{b}_{l}=\rho s^{2}\hat{U}^{M}_{l},$ (1ama) $\displaystyle
K_{r_{1}q_{1}}\dfrac{\partial^{2}\hat{\Theta}^{M}}{\partial
x_{r_{1}}x_{q_{1}}}-s\alpha_{lq_{1}}\dfrac{\partial\hat{U}^{M}_{l}}{\partial
x_{q_{1}}}+\hat{r}=ps\hat{\Theta}^{M},$ (1amb)

where the overall constitutive tensors are defined as

$\displaystyle C_{lr_{1}pq_{1}}=\left\langle
B^{(1,0)}_{ijlr_{1}}C^{m}_{ijhk}B^{(1,0)}_{hkpq_{1}}\right\rangle,$ (1ana)
$\displaystyle K_{r_{1}q_{1}}=\left\langle
A^{(1,0)}_{ir_{1}}K^{m}_{ij}A^{(1,0)}_{jq_{1}}\right\rangle,$ (1anb)
$\displaystyle\alpha_{lq_{1}}=\left\langle
B^{(1,0)}_{ijlq_{1}}\alpha^{m}_{ij}-B^{(1,0)}_{hklq_{1}}C^{m}_{ijhk}\tilde{B}^{(1,0)}_{ij}\right\rangle,$
(1anc) $\displaystyle p=\langle
p^{m}+2\tilde{B}^{(1,0)}_{ij}\alpha^{m}_{ij}+\tilde{B}^{(1,0)}_{ij}C^{m}_{ijhk}\tilde{B}^{(1,0)}_{hk}\rangle,$
(1and) $\displaystyle\rho=\langle\rho^{m}\rangle.$ (1ane)

Note that, at this stage, the equations (1ama) and (1amb) are the macro-scale
equivalent to the field equations (1ka) and (1kb) on the transformed Laplace
space. Moreover, by applying the inverse Laplace transform (1i) to them, and
recalling the identity (1n), the field equations at the macro-scale
corresponding to the equations (1ga) and (1gb) are recast on the time domain
as

$\displaystyle C_{lr_{1}pq_{1}}\dfrac{\partial^{2}\ddot{U}^{M}_{p}}{\partial
x_{q_{1}}\partial
x_{r_{1}}}-\alpha_{lq_{1}}\dfrac{\partial\dot{\Theta}^{M}}{\partial
x_{q_{1}}}+\hat{b}_{l}=\rho\ddot{U}^{M}_{l},$ (1aoa)
$\displaystyle\mathcal{L}^{-1}\left\\{K_{r_{1}q_{1}}\right\\}\ast\dfrac{\partial^{2}{\Theta}^{M}}{\partial
x_{r_{1}}x_{q_{1}}}-\alpha_{lq_{1}}\dfrac{\partial\dot{U}^{M}_{l}}{\partial
x_{q_{1}}}+r=p\dot{\Theta}^{M}.$ (1aob)

##  4 Wave propagation in homogenized continuum

In this Section, we exploit the wave propagation along the homogenized
continuum, which is equivalent to the started thermoelastic periodic material.
Such procedure consists in applying the Fourier transform to equations (1ama)
and (1amb), with respect to the macroscopic variable
$\textbf{x}\in\mathbb{R}^{3},$ in order to derive these field equations at the
macro-scale defined over the complex frequency $s,$ and also defined over the
wave vector $\boldsymbol{k}=(k_{1},k_{2},k_{3})^{T}\in\mathbb{C}^{3}$
($k_{1},k_{2},k_{3}$ are the wave numbers and $T$ identifies transpose
vector), in other terms, the thermoelastic field equations (1ga) and (1gb),
initially defined on the space-time domain
$(\textbf{x},t)\in\mathbb{R}^{3}\times\mathbb{R},$ will be taken to the wave
vector-frequency domain $(\boldsymbol{k},s)\in\mathbb{C}^{3}\times\mathbb{C}$
[81].

Proceeding the same way as made in Sec. 1, we begin defining the space Fourier
transform and its properties that will be used throughout the work [79]. Let
$f:\mathbb{R}^{3}\rightarrow\mathbb{R}$ be an arbitrary function defined over
the domain $\textbf{x}\in\mathbb{R}^{3},$ the complex space Fourier transform
is defined as,

$\mathcal{F}(f(\textbf{x}))=\check{f}(\boldsymbol{k})=\int_{\mathbb{R}^{3}}f(\textbf{x})\mathrm{e}^{-\mathrm{i}\boldsymbol{k}\cdot\textbf{x}}\mathrm{d}\textbf{x},\quad\boldsymbol{k}\in\mathbb{C}^{3},$
(1ap)

where $\check{f}:\mathbb{C}^{3}\to\mathbb{C}.$

The Fourier transform of the derivative of order $n+m$ of $f$ with respect to
$x^{n}_{j}$ and $x^{m}_{r}$ is given by the
identity$\mathcal{F}\left\\{\partial^{n+m}f(\textbf{x})\partial
x^{n}_{j}\partial
x^{m}_{r}\right\\}=(\mathrm{i}{k}_{j})^{n}(\mathrm{i}{k}_{r})^{m}\check{f}(\boldsymbol{k}).$

Transforming the macro-scale field equations (1ama) and (1amb) into the
Fourier transformed space, one has

$\displaystyle
C_{lr_{1}pq_{1}}(\mathrm{i}{k}_{q_{1}})(\mathrm{i}{k}_{r_{1}})\check{\hat{U}}^{M}_{p}-\alpha_{lq_{1}}(\mathrm{i}{k}_{q_{1}})\check{\hat{\Theta}}^{M}+\check{\hat{b}}_{l}=\rho
s^{2}\check{\hat{U}}^{M}_{l},$ (1aqa) $\displaystyle
K_{r_{1}q_{1}}(\mathrm{i}{k}_{q_{1}})(\mathrm{i}{k}_{r_{1}})\check{\hat{\Theta}}^{M}-s\alpha_{lq_{1}}(\mathrm{i}{k}_{q_{1}})\check{\hat{U}}^{M}_{l}+\check{\hat{r}}=ps\check{\hat{\Theta}}^{M}.$
(1aqb)

In case of free wave propagation, we take the external influences (source
terms) being zero, such as the heat source equivalent $\check{\hat{r}},$ and
the body force equivalent $\check{\hat{\boldsymbol{b}}},$ so the equations
above might be rephrased to the compact form as

$\displaystyle\left[\bar{\mathfrak{C}}(\boldsymbol{k}\otimes\boldsymbol{k})+\rho
s^{2}\boldsymbol{I}\right]\cdot\check{\hat{\boldsymbol{U}}}^{M}+\mathrm{i}(\boldsymbol{\alpha}\boldsymbol{k})\check{\hat{\Theta}}^{M}=\boldsymbol{0},$
(1ara)
$\displaystyle\left[\mathbf{K}:(\boldsymbol{k}\otimes\boldsymbol{k})+ps\right]\check{\hat{\Theta}}^{M}+\mathrm{i}s(\boldsymbol{\alpha}\boldsymbol{k})\cdot\check{\hat{\boldsymbol{U}}}^{M}=0,$
(1arb)

where the tensor $\bar{\mathfrak{C}}$ corresponds to the shift
$\bar{\mathfrak{C}}=\bar{C}_{lpr_{1}q_{1}}\textbf{e}_{l}\otimes\textbf{e}_{p}\otimes\textbf{e}_{r_{1}}\otimes\textbf{e}_{q_{1}},$
with $\bar{C}_{lpr_{1}q_{1}}=C_{lpr_{1}q_{1}}.$ Note that, the system of
linear equations formed by (1ara) and (1arb) may be rewritten in a matricial
form as

$\left[\begin{array}[]{cc}\bar{\mathfrak{C}}(\boldsymbol{k}\otimes\boldsymbol{k})+\rho
s^{2}\boldsymbol{I}&\mathrm{i}\boldsymbol{k}\boldsymbol{\alpha}\\\
\mathrm{i}s\boldsymbol{k}\boldsymbol{\alpha}&\mathbf{K}:(\boldsymbol{k}\otimes\boldsymbol{k})+ps\end{array}\right]\left[\begin{array}[]{c}\check{\hat{\boldsymbol{U}}}^{M}\\\
\check{\hat{\Theta}}^{M}\end{array}\right]=\left[\begin{array}[]{c}\boldsymbol{0}\\\
0\end{array}\right].$ (1ay)

The matricial system (1ay), is an eigenvector-eigenvalue problem that provides
the frequency spectrum as an implicit function of the wave vector
$\boldsymbol{k}$ and the complex frequency $s,$ which is given by the
frequency equation
$F(\boldsymbol{k},s)=\textrm{det}([\boldsymbol{H}(\boldsymbol{k},s)])=0,$ also
known as implicit dispersion relation, where the matrices in the system (1ay)
are named as
$[\boldsymbol{H}(\boldsymbol{k},s)][\boldsymbol{V}]=[\boldsymbol{0}].$ Since
the implicit function $F(\boldsymbol{k},s)=0$ is a complex function, it can be
rephrased in terms of its real part and imaginary part as
$F(\boldsymbol{k},s)=\Re e(F(\boldsymbol{k},s))+\mathrm{i}\Im
m(F(\boldsymbol{k},s))=0,$ or simply

$\left\\{\begin{array}[]{lcr}\Re e(F(\boldsymbol{k},s))&=&0\\\ \Im
m(F(\boldsymbol{k},s))&=&0\end{array}\right.,$ (1bb)

which represents the collection of branches (frequency spectrum) by the hyper-
surfaces $\Re e(F(\boldsymbol{k},s))$ and $\Im m(F(\boldsymbol{k},s))$ into
the space $(\boldsymbol{k},s)\in\mathbb{C}^{4},$ such hyper-surface
$F(\boldsymbol{k},s)$ is also seen in mathematics as a level set of the
complexed-valued function $F$. As [81] quotes, from a mathematical standpoint,
there is no reason to choose real over complex values for either the frequency
or the wave vector. What matters is actually the physical meaning that might
be associated to them. Having described the implicit dispersion function, from
this stage on we emphasize two possible lanes that might provide us the
dispersion spectrum, and finally analyse the wave propagation either with
spatial damping or with the damping in time, [82].

In the first case, let us consider the complex angular frequency as
$s=\mathrm{i}\omega$ with $\omega\in\mathbb{R}$ in the implicit dispersion
relation (1bb). Thus the dispersion relation dependence goes over to the
complex wave vector $\boldsymbol{k}$ and the angular frequency $\omega.$
Therefore, the dispersion relation associated to generic inhomogeneous
harmonic waves holds

$\left\\{\begin{array}[]{lcr}\Re e(F(\Re e(\boldsymbol{k}),\Im
m(\boldsymbol{k}),\omega))&=&0\\\ \Im m(F(\Re e(\boldsymbol{k}),\Im
m(\boldsymbol{k}),\omega))&=&0\end{array}\right.,$ (1be)

where the intersection between these two hyper-surfaces defines the frequency
spectrum of the material as an hyper-curve immersed in $\mathbb{R}^{7}.$ More
specifically, such hyper-curve describes the relation between the angular
frequency $\omega$ and the complex vector $\boldsymbol{k},$ which can be
expressed as $\boldsymbol{k}=\Re e(\boldsymbol{k})+\mathrm{i}\Im
m(\boldsymbol{k})=\|\Re e(\boldsymbol{k})\|\boldsymbol{n}_{r}+\mathrm{i}\|\Im
m(\boldsymbol{k})\|\boldsymbol{n}_{i},$ where $\boldsymbol{n}_{r},$
$\boldsymbol{n}_{i}$ are versors (i.e
$\|\boldsymbol{n}_{r}\|=\|\boldsymbol{n}_{i}\|=1$ and
$\boldsymbol{n}_{r},\boldsymbol{n}_{i}\in\mathbb{R}^{3}$), representing the
direction of the normal to planes of constant phase and planes of constant
amplitude of the propagating wave, respectively. Moreover, in the particular
case of the homogeneous harmonic waves characterized by
$\boldsymbol{n}_{r}=\boldsymbol{n}_{i}=\boldsymbol{n}$ and for the complex
wave expressed in the form $\boldsymbol{k}=(\|\Re e(\boldsymbol{k})\|+i\|\Im
m(\boldsymbol{k})\|)\boldsymbol{n}=\kappa\boldsymbol{n},$ the dispersion
relation is obtained by specializing (1be) in the form

$\left\\{\begin{array}[]{lcr}\Re e(F(\Re e(\kappa),\Im
m(\kappa),\omega))&=&0\\\ \Im m(F(\Re e(\kappa),\Im
m(\kappa),\omega))&=&0\end{array}\right.,$ (1bh)

a three dimensional hyper-surface whose intersection is a curve immersed in
$\mathbb{R}^{3}.$

Nonetheless, in the second case of wave propagation with attenuation in time,
let us assume in the set of equations (1bb) the complex frequency being as
$s=\Re e(s)+\mathrm{i}\Im m(s),$ the wave vector such that
$\boldsymbol{k}\in\mathbb{R}^{3},$ and let us also fix a direction of the wave
vector $\boldsymbol{k}=\kappa\boldsymbol{n},$ where $\boldsymbol{n}$ is its
versor and $\kappa=||\boldsymbol{k}||\in\mathbb{R},$ thus (1bb) becomes

$\left\\{\begin{array}[]{lcr}\Re e(F(\kappa,\Re e(s),\Im m(s)))&=&0\\\ \Im
m(F(\kappa,\Re e(s),\Im m(s)))&=&0\end{array}\right.,$ (1bk)

and hence taking the intersection between the two hyper-surfaces from (1bk),
one yields a hyper-curve $s(\kappa)$-curves immersed in $\mathbb{R}^{3},$
given a certain propagation direction $\boldsymbol{n}.$

##  5 Illustrative example: homogenization procedure on a bi-phase orthotropic
layered material

In order to contrast the results found by the analysis of the heterogeneous
approach stressed in the Sec. B.1. of the Supplementary Material, the
homogenized formulation obtained in the Sec. 5 is now applied for a layered
two-dimensional infinite thermoelastic body with orthotropic phases having the
orthotropy axes parallel to the layering direction $\boldsymbol{e_{1}}.$
Analytical precise formulae for the overall elastic, thermal dilatation and
thermal conduction tensors are determined. Subsequently, we compare the
results of the heterogeneous model developed in the Sec. B.1., against those
from the model obtained proceeding with the first order homogenization
technique, both under the same hypothesis, by finding their dispersions
spectrum seen in the Sec. 4.

###  5.1 Perturbation functions and overall constitutive tensors

Let us assume a layered body composed by two phases where the first order
homogenization process will be applied on the system. As may seen in Fig. 2,
let the domain $\mathcal{A}$ made by two layered materials be defined having
thickness $\mathfrak{s}_{1}$ and $\mathfrak{s}_{2},$ where
$d_{2}=\varepsilon=\mathfrak{s}_{1}+\mathfrak{s}_{2},$ and
$\eta=\mathfrak{s}_{1}/\mathfrak{s}_{2}.$ Also the domain displays orthotropic
phases, which it coincides with the axis that determines the layering
direction $\boldsymbol{e}_{1}.$ The perturbation functions at the order
$\varepsilon^{-1},$ $N_{hpq_{1}}^{(1,0)}$ $\tilde{N}_{h}^{(1,0)},$
$M_{q_{1}}^{(1,0)}$ are analytically obtained by the solutions of the three
cell problems presented in (1u) formulated in Sec. 5, along with the interface
conditions in (1v). Due to the microstructure symmetry, these functions rely
exclusively on the microscopic (fast) component $\xi_{2},$ which is
perpendicular to the layering direction $\boldsymbol{e}_{1}.$ Considering only
non-zero perturbations functions, they are given by

$\displaystyle N_{112}^{(1,0)_{1}}$
$\displaystyle=N_{121}^{(1,0)_{1}}=\dfrac{\left(C_{1212}^{{}_{(2)}}-C_{1212}^{{}_{(1)}}\right)\xi_{2}^{{}_{(1)}}}{C_{1212}^{{}_{(2)}}\eta+C_{1212}^{{}_{(1)}}},$
$\displaystyle N_{112}^{(1,0)_{2}}$
$\displaystyle=N_{121}^{(1,0)_{2}}=\dfrac{\eta\left(C_{1212}^{{}_{(2)}}-C_{1212}^{{}_{(1)}}\right)\xi_{2}^{{}_{(2)}}}{C_{1212}^{{}_{(2)}}\eta+C_{1212}^{{}_{(1)}}},$
(1bla) $\displaystyle N_{211}^{(1,0)_{1}}$
$\displaystyle=\dfrac{\left(C_{1122}^{{}_{(2)}}-C_{1122}^{{}_{(1)}}\right)\xi_{2}^{{}_{(1)}}}{C_{2222}^{{}_{(2)}}\eta+C_{2222}^{{}_{(1)}}},$
$\displaystyle N_{211}^{(1,0)_{2}}$
$\displaystyle=\dfrac{\eta\left(C_{1122}^{{}_{(2)}}-C_{1122}^{{}_{(1)}}\right)\xi_{2}^{{}_{(2)}}}{C_{2222}^{{}_{(2)}}\eta+C_{2222}^{{}_{(1)}}},$
(1blb) $\displaystyle N_{222}^{(1,0)_{1}}$
$\displaystyle=\dfrac{\left(C_{2222}^{{}_{(2)}}-C_{2222}^{{}_{(1)}}\right)\xi_{2}^{{}_{(1)}}}{C_{2222}^{{}_{(2)}}\eta+C_{2222}^{{}_{(1)}}},$
$\displaystyle N_{222}^{(1,0)_{2}}$
$\displaystyle=\dfrac{\eta\left(C_{2222}^{{}_{(2)}}-C_{2222}^{{}_{(1)}}\right)\xi_{2}^{{}_{(2)}}}{C_{2222}^{{}_{(2)}}\eta+C_{2222}^{{}_{(1)}}},$
(1blc) $\displaystyle\tilde{N}_{2}^{(1,0)_{1}}$
$\displaystyle=\dfrac{\left(\alpha_{22}^{{}_{(1)}}-\alpha_{22}^{{}_{(2)}}\right)\xi_{2}^{{}_{(1)}}}{C_{2222}^{{}_{(2)}}\eta+C_{2222}^{{}_{(1)}}},$
$\displaystyle\tilde{N}_{2}^{(1,0)_{2}}$
$\displaystyle=\dfrac{\eta\left(\alpha_{22}^{{}_{(2)}}-\alpha_{22}^{{}_{(1)}}\right)\xi_{2}^{{}_{(2)}}}{C_{2222}^{{}_{(2)}}\eta+C_{2222}^{{}_{(1)}}},$
(1bld) $\displaystyle M_{2}^{(1,0)_{1}}$
$\displaystyle=\dfrac{\left(K_{22}^{{}_{(2)}}-K_{22}^{{}_{(1)}}\right)\xi_{2}^{{}_{(1)}}}{K_{22}^{{}_{(2)}}\eta+K_{22}^{{}_{(1)}}},$
$\displaystyle M_{2}^{(1,0)_{2}}$
$\displaystyle=\dfrac{\eta\left(K_{22}^{{}_{(1)}}-K_{22}^{{}_{(2)}}\right)\xi_{2}^{{}_{(2)}}}{K_{22}^{{}_{(2)}}\eta+K_{22}^{{}_{(1)}}},$
(1ble)

where the superscript $i=\\{1,2\\}$ denotes for the phase 1 and the phase 2,
respectively. Moreover, the dimensionless vertical variables
$\xi_{2}^{{}_{(1)}}$ and $\xi_{2}^{{}_{(2)}}$ centred in each layer are such
that, $\xi_{2}^{{}_{(1)}}\in\left[-\eta/2(\eta+1),\eta/2(\eta+1)\right]$ and
$\xi_{2}^{{}_{(2)}}\in\left[-1/2(\eta+1),1/2(\eta+1)\right],$ agreeing with
the Fig. 2.

Figure 2: Heterogeneous bi-dimensional domain bi-phase layered periodic cell.

In particular, once the perturbation functions $M_{q1}^{(1,0)}$ depend on the
complex angular frequency $s$ (see Eqs. (1l) and (1cob)), such dependency is
evaluated, but dimensionless. As may be seen in Fig. 3, the complex absolute
value of the dimensionless perturbation function, which is
$\|M_{2}^{(1,0)}\|,$ is analytically computed by the equations in (1ble), with
respect to the phase $1$ and the phase $2,$ taking the following values for
dimensionless parameters $\eta=1,$ $\tau^{{}_{(2)}}/\tau^{{}_{(1)}}=3$ and
$\bar{K}_{22}^{{}_{(2)}}/\bar{K}_{22}^{{}_{(1)}}=10.$ Such a function depends
on the fast variable $\xi_{2},$ which is perpendicular to the transversal
direction $\boldsymbol{e}_{1}$ and to the dimensionless complex parameter
$s\tau^{{}_{(1)}}=\Re e(s\tau^{{}_{(1)}})+\mathrm{i}\Im m(s\tau^{{}_{(1)}}).$
It must be noted in Fig. 3$(a)$ that the function $\|M_{2}^{(1,0)}(\xi_{2},\Re
e(s\tau^{{}_{(1)}}))\|$ presents a singularity along the variable $\xi_{2},$
whereas the function $\|M_{2}^{(1,0)}(\xi_{2},\Im m(s\tau^{{}_{(1)}}))\|$
displays the absence of any sort of singularities regardless the variable.

As previously mentioned, we also derive the non-vanishing overall elastic,
thermal dilatation and thermal conduction tensors corresponding to a first
order equivalent continuum. Thus, from the perturbation functions (1bla),
(1blb) and (1blc) are used into equation (1ana), hence the overall elastic
tensors $C_{lr_{1}pq_{1}}$ are expressed as

Figure 3: Dimensionless complex absolute value of the perturbation function
$M_{2}^{(1,0)}$ for the dimensionless parameters assumed $\eta=1,$
$\tau^{{}_{(2)}}/\tau^{{}_{(1)}}=3$ and
$\bar{K}_{22}^{{}_{(2)}}/\bar{K}_{22}^{{}_{(1)}}=10.$ $(a)$
$\|M_{2}^{(1,0)}\|$ vs. $\Re e(s\tau^{{}_{(1)}})$ and $\xi_{2};$ $(b)$
$\|M_{2}^{(1,0)}\|$ vs. $\Im m(s\tau^{{}_{(1)}})$ and $\xi_{2}.$

$\displaystyle
C_{1111}=\dfrac{\eta^{2}C_{1111}^{{}_{(1)}}C^{{}_{(2)}}_{2222}+\eta\left(C^{{}_{(2)}}_{1111}C^{{}_{(2)}}_{2222}+C^{{}_{(2)}}_{1111}C^{{}_{(1)}}_{2222}-\left(C_{1122}^{{}_{(1)}}-C_{1122}^{{}_{(2)}}\right)^{2}\right)+C_{1111}^{{}_{(2)}}C^{{}_{(1)}}_{2222}}{(\eta+1)\left(C^{{}_{(1)}}_{2222}+\eta
C^{{}_{(2)}}_{2222}\right)},$ (1bma) $\displaystyle C_{1122}=\dfrac{\eta
C^{{}_{(1)}}_{1122}C^{{}_{(2)}}_{2222}+C^{{}_{(2)}}_{1122}C^{{}_{(1)}}_{2222}}{C^{{}_{(1)}}_{2222}+\eta
C^{{}_{(2)}}_{2222}},$ (1bmb) $\displaystyle
C_{1212}=\dfrac{(\eta+1)C^{{}_{(1)}}_{1212}C^{{}_{(2)}}_{1212}}{C^{{}_{(1)}}_{1212}+\eta
C^{{}_{(2)}}_{1212}},$ (1bmc) $\displaystyle
C_{2222}=\dfrac{(\eta+1)C^{{}_{(1)}}_{2222}C^{{}_{(2)}}_{2222}}{C^{{}_{(1)}}_{2222}+\eta
C^{{}_{(2)}}_{2222}}.$ (1bmd)

Similarly, for the thermal dilatation components follows that from the
perturbation functions (1bla) to (1bld) into equations (1ana) and (1anc), one
provides

$\displaystyle\alpha_{11}=\dfrac{\eta^{2}C^{{}_{(2)}}_{2222}\alpha^{{}_{(1)}}_{11}+\eta\left(C^{{}_{(1)}}_{2222}\alpha^{{}_{(1)}}_{11}+C^{{}_{(2)}}_{2222}\alpha^{{}_{(2)}}_{11}-\left(\alpha^{{}_{(1)}}_{22}-\alpha^{{}_{(2)}}_{22}\right)\left(C^{{}_{(1)}}_{1122}-C^{{}_{(2)}}_{1122}\right)\right)+C^{{}_{(1)}}_{2222}\alpha^{{}_{(2)}}_{11}}{(\eta+1)(C^{{}_{(1)}}_{2222}+\eta
C^{{}_{(2)}}_{2222})},$ (1bna) $\displaystyle\alpha_{22}=\dfrac{\eta
C^{{}_{(2)}}_{2222}\alpha^{{}_{(1)}}_{22}+\alpha^{{}_{(2)}}_{22}C^{{}_{(1)}}_{2222}}{\eta
C^{{}_{(2)}}_{2222}+C^{{}_{(1)}}_{2222}}.$ (1bnb)

The specific heat $p$ and the mass density $\rho$ in the equations (1and) and
(1ane) respectively, are obtained following the same steps made latter for the
elastic, and thermal dilatations tensors, giving

$\displaystyle
p=\dfrac{p^{{}_{(1)}}\eta^{2}C^{{}_{(2)}}_{2222}+\left(p^{{}_{(1)}}C_{2222}^{{}_{(1)}}+p^{{}_{(2)}}C_{2222}^{{}_{(2)}}+2\left(\alpha_{22}^{{}_{(1)}}-\alpha_{22}^{{}_{(2)}}\right)^{2}\right)\eta+p^{{}_{(2)}}C_{2222}^{{}_{(1)}}}{\left(\eta
C^{{}_{(2)}}_{2222}+C^{{}_{(1)}}_{2222}\right)\left(\eta+1\right)},$ (1boa)
$\displaystyle\rho=\dfrac{\eta\rho^{{}_{(1)}}+\rho^{{}_{(2)}}}{\eta+1}.$
(1bob)

Still in this matter, from the perturbation functions (1ble) and equation
(1anb), the components of the overall thermal conduction tensor are

$\displaystyle K_{11}(s)=\dfrac{\eta
K^{{}_{(1)}}_{11}(s)+K^{{}_{(2)}}_{11}(s)}{\eta+1},$ (1bpa) $\displaystyle
K_{22}(s)=\dfrac{(\eta+1)K^{{}_{(1)}}_{22}(s)K^{{}_{(2)}}_{22}(s)}{\eta
K^{{}_{(2)}}_{22}(s)+K^{{}_{(1)}}_{22}(s)}.$ (1bpb)

Consequently, from the Eq. (1l), the corresponding parts of the overall
thermal conduction tensor namely $K_{11}(s)$ and $K_{22}(s),$ are also
dependent on the complex angular frequency $s,$ which ends up being rewritten
in terms of the relaxation times $\tau^{{}_{(}1)},$ $\tau^{{}_{(}2)},$ the
specific heat $T_{0}$ and also the thermal conductivity tensors defined at the
begining of the theory. After the respective substitutions of (1l) to both
equations in (1bp), the equivalent components for overall thermal conduction
tensor are, respectively:

Figure 4: Dimensionless complex absolute value of the non-vanishing overall
thermal conduction tensor components for the dimensionless parameters
$\eta=1,$ $\tau^{{}_{(2)}}/\tau^{{}_{(1)}}=3$ and
$\bar{K}_{22}^{{}_{(2)}}/\bar{K}_{22}^{{}_{(1)}}=10.$ $(a)$
$\|{K}_{11}T_{0}/\bar{K}_{11}^{{}_{(1)}}\|$ vs. $\Re
e(s\tau^{{}_{(1)}})\times\Im m(s\tau^{{}_{(1)}});$ $(b)$
$\|\bar{K}_{22}T_{0}/{K}_{22}^{{}_{(1)}}\|$ vs. $\Re
e(s\tau^{{}_{(1)}})\times\Im m(s\tau^{{}_{(1)}}).$

$\displaystyle{K}_{11}(s)=\dfrac{\eta\bar{K}^{{}_{(1)}}_{11}(s\tau^{{}_{(2)}}+1)+\bar{K}^{{}_{(2)}}_{11}(s\tau^{{}_{(1)}}+1)}{T_{0}(\eta+1)(s\tau^{{}_{(1)}}+1)(s\tau^{{}_{(2)}}+1)},$
(1bqa)
$\displaystyle{K}_{22}(s)=\dfrac{(\eta+1)\bar{K}^{{}_{(1)}}_{22}\bar{K}^{{}_{(2)}}_{22}}{T_{0}\left(\eta\bar{K}^{{}_{(2)}}_{22}(s\tau^{{}_{(1)}}+1)+\bar{K}^{{}_{(1)}}_{22}(s\tau^{{}_{(2)}}+1)\right)}.$
(1bqb)

The Fig. 4 displays the complex absolute value of the dimensionless overall
thermal conduction tensor components
$\|{K}_{11}T_{0}/\bar{K}_{11}^{{}_{(1)}}\|$ and
$\|{K}_{22}T_{0}/\bar{K}_{22}^{{}_{(1)}}\|$ in terms of the dimensionless
complex frequency $(\Re e(s\tau^{{}_{(1)}}),\Im m(s\tau^{{}_{(1)}})),$ for the
same dimensionless parameters $\eta=1,$ $\tau^{{}_{(2)}}/\tau^{{}_{(1)}}=3$
and $\bar{K}_{22}^{{}_{(2)}}/\bar{K}_{22}^{{}_{(1)}}=10$ already assumed when
it was analysed the perturbation function $M_{2}^{(1,0)}.$ Still in the regard
of the Fig. 4, note that the Eq. (1bqa) presents a quadratic behaviour in its
denominator and hence, as expected, it does display two poles, whereas the
equation for $K_{22}$, namely Eq. (1bqb), has a linear nature in its
denominator, giving only one pole.

###  5.2 Comparative analysis: heterogeneous material vs. homogenized solid

Under the same hypothesis assumed in the beginning of this Section, herein we
proceed with a comparative analysis between the results obtained from the
heterogeneous modelling procedure developed in the Sec. B.1. via Floquet-Bloch
theory, and the homogenized first order model approached in the previous
sections of the generalized thermoelastic problem with the derived overall
elastic, thermal dilatation and thermal conduction constants.

For the considered two-dimensional body with the orthotropy axes perpendicular
to the layering direction $\boldsymbol{e_{2}}$ and the wave vector taken such
as $k_{1}=0,$ i.e the propagation axes is exclusively along the
$\boldsymbol{e_{2}}.$ For this scenario we consider an uncoupled problem which
means that the dilatation tensor is assumed zero. Having said so, on the one
hand by plugging overall constitutive tensors displayed in the set of Eqs.
(1bm), (1bo) and (1bq) into the homogenized specialized governing equations
from (1aoa), (1aob) and following the procedure reported in Sec. 4, the
dispersive relation (1bk) (for $\kappa=k_{2}$) provides the explicit functions

$k_{2}(s)=\pm\dfrac{\mathrm{i}s\sqrt{C_{1212}^{{}_{(1)}}C_{1212}^{{}_{(2)}}\left(\eta\rho^{{}_{(1)}}+\rho^{{}_{(2)}}\right)\left(\eta
C_{1212}^{{}_{(2)}}+C_{1212}^{{}_{(1)}}\right)}}{C_{1212}^{{}_{(1)}}C_{1212}^{{}_{(2)}}(\eta+1)},$
(1br)

or in its inverse form

$s(k_{2})=\pm\dfrac{\mathrm{i}k_{2}(\eta+1)\sqrt{C_{1212}^{{}_{(1)}}C_{1212}^{{}_{(2)}}\left(\eta\rho^{{}_{(1)}}+\rho^{{}_{(2)}}\right)\left(\eta
C_{1212}^{{}_{(2)}}+C_{1212}^{{}_{(1)}}\right)}}{\left(\eta\rho^{{}_{(1)}}+\rho^{{}_{(2)}}\right)\left(\eta
C_{1212}^{{}_{(2)}}+C_{1212}^{{}_{(1)}}\right)},$ (1bs)

corresponding to the dispersion relation associated to the shear waves, giving
rise to the standard non-dispersive behaviour. Similarly

$k_{2}(s)=\pm\dfrac{\mathrm{i}s\sqrt{C_{2222}^{{}_{(1)}}C_{2222}^{{}_{(2)}}\left(\eta\rho^{{}_{(1)}}+\rho^{{}_{(2)}}\right)\left(\eta
C_{2222}^{{}_{(2)}}+C_{2222}^{{}_{(1)}}\right)}}{C_{2222}^{{}_{(1)}}C_{2222}^{{}_{(2)}}(\eta+1)},$
(1bt)

with the inversion function as

$s(k_{2})=\pm\dfrac{\mathrm{i}k_{2}(\eta+1)\sqrt{C_{2222}^{{}_{(1)}}C_{2222}^{{}_{(2)}}\left(\eta\rho^{{}_{(1)}}+\rho^{{}_{(2)}}\right)\left(\eta
C_{2222}^{{}_{(2)}}+C_{2222}^{{}_{(1)}}\right)}}{\left(\eta\rho^{{}_{(1)}}+\rho^{{}_{(2)}}\right)\left(\eta
C_{2222}^{{}_{(2)}}+C_{2222}^{{}_{(1)}}\right)},$ (1bu)

corresponding to the dispersion relation associated to the compressional
waves, providing once again the standard non-dispersive behaviour. Finally,
for dispersion relation associated to the thermal waves,

$k_{2}(s)=\pm\dfrac{\mathrm{i}\sqrt{\mathscr{Z}_{0}\left[s^{2}\left(\eta\bar{K}_{22}^{{}_{(2)}}\tau^{{}_{(1)}}+\bar{K}_{22}^{{}_{(1)}}\tau^{{}_{(2)}}\right)+s\left(\eta\bar{K}_{22}^{{}_{(2)}}+\bar{K}_{22}^{{}_{(1)}}\right)\right]}}{\bar{K}^{{}_{(1)}}_{22}\bar{K}^{{}_{(2)}}_{22}(\eta+1)\left(\eta
C_{2222}^{{}_{(2)}}+C_{2222}^{{}_{(1)}}\right)},$ (1bv)

where
$\mathscr{Z}_{0}=T_{0}\left(\eta^{2}p^{{}_{(1)}}C_{2222}^{{}_{(2)}}+\eta\left(p^{{}_{(1)}}C_{2222}^{{}_{(1)}}+2\left(\alpha_{22}^{{}_{(1)}}-\alpha_{22}^{{}_{(2)}}\right)^{2}\right)+p^{{}_{(2)}}C_{2222}^{{}_{(1)}}\right),$
or in its inverse form

$s(k_{2})=\dfrac{\mathscr{Z}_{0}\left(\eta\bar{K}_{22}^{{}_{(2)}}+\bar{K}_{22}^{{}_{(1)}}\right)\pm\sqrt{\mathscr{Z}_{1}(k_{2})}}{2\mathscr{Z}_{0}\left(\eta\bar{K}_{22}^{{}_{(2)}}\tau^{{}_{(1)}}+\bar{K}_{22}^{{}_{(1)}}\tau^{{}_{(2)}}\right)},$
(1bw)

with the function $\mathscr{Z}_{1}(k_{2})$ defined as

$\mathscr{Z}_{1}(k_{2})=-4k^{2}_{2}\mathscr{Z}_{0}\left(\eta\bar{K}_{22}^{{}_{(2)}}\tau^{{}_{(1)}}+\bar{K}_{22}^{{}_{(1)}}\tau^{{}_{(2)}}\right)(\eta+1)^{2}\left(\eta
C_{2222}^{{}_{(2)}}+C_{2222}^{{}_{(1)}}\right)^{2}\left(\bar{K}^{{}_{(1)}}_{22}\bar{K}^{{}_{(2)}}_{22}\right)^{2}+\mathscr{Z}_{0}^{2}\left(\eta\bar{K}_{22}^{{}_{(2)}}+\bar{K}_{22}^{{}_{(1)}}\right)^{2}.$
(1bx)

It is important to notice here that, the dispersion relations (1bv) and (1bw)
are expressed, as expected, in terms of both phases of the micro-thermal-
conductivity tensor components $\bar{K}^{{}_{(1)}}_{22}$ and
$\bar{K}^{{}_{(2)}}_{22},$ as well as the relaxation times $\tau^{{}_{(1)}}$
and $\tau^{{}_{(2)}}.$

Figure 5: Dimensionless dispersion functions associated to thermal waves
characterized when $k_{1}=0.$ Comparison between a first order homogenized
model (light blue curve and light red curve) with its respective heterogeneous
one (dark blue curve and dark red curve), assuming the numerical values for
the dimensionless parameters $\eta=1,$ $\tau^{{}_{(2)}}/\tau^{{}_{(1)}}=3,$
$\bar{K}_{22}^{{}_{(2)}}/\bar{K}_{22}^{{}_{(1)}}=3,$
$p^{{}_{(2)}}/p^{{}_{(1)}}=3$ and
$p^{{}_{(1)}}T_{0}/\left(\tau^{{}_{(1)}}\bar{K}_{22}^{{}_{(1)}}\right)=1.$
$(a)$ $\omega\tau^{{}_{(1)}}$ vs. $\Re e(k_{2}\varepsilon)\times\Im
m(k_{2}\varepsilon);$ $(b)$ zoomed view of the angular frequency spectrum
$\omega\tau^{{}_{(1)}}$ vs. $\Re e(k_{2}\varepsilon)\times\Im
m(k_{2}\varepsilon);$ $(c)$ view of the plane $\omega\tau^{{}_{(1)}}\times\Re
e(k_{2}\varepsilon);$ $(d)$ view of the plane $\omega\tau^{{}_{(1)}}\times\Im
m(k_{2}\varepsilon).$

In respect of characterizing harmonic waves in the homogenized continuum, it
is worth noting that for these above equations and their respective inverse
forms allow us to either describe them by spatial damping or time damping as
previously seen in Sec. 4.

On the other hand, followed from the Secs. B.1. and B.2. of the Supplementary
Material, the dispersion relations related to the heterogeneous bi-dimensional
body with the orthotropy axes perpendicular to the layering direction
$\boldsymbol{e_{2}}$ and the wave vector taken such as $k_{1}=0,$ with
$d_{2}=d^{{}_{(1)}}+d^{{}_{(2)}}$ are given by

Figure 6: Dimensionless dispersion functions associated to thermal waves
characterized when $k_{1}=0.$ Comparison between homogenized models (light
curves) and their respective heterogeneous ones (dark curves) given different
values of the parameters. $(a)$ setting $\eta=1,$
$p^{{}_{(1)}}T_{0}/\left(\tau^{{}_{(1)}}\bar{K}_{22}^{{}_{(1)}}\right)=1$ and
varying the parameters as $\tau^{{}_{(2)}}/\tau^{{}_{(1)}}=3,$
$\bar{K}_{22}^{{}_{(2)}}/\bar{K}_{22}^{{}_{(1)}}=3,$
$p^{{}_{(2)}}/p^{{}_{(1)}}=3$ (red curves),
$\tau^{{}_{(2)}}/\tau^{{}_{(1)}}=5,$
$\bar{K}_{22}^{{}_{(2)}}/\bar{K}_{22}^{{}_{(1)}}=5,$
$p^{{}_{(2)}}/p^{{}_{(1)}}=5$ (green curves),
$\tau^{{}_{(2)}}/\tau^{{}_{(1)}}=10,$
$\bar{K}_{22}^{{}_{(2)}}/\bar{K}_{22}^{{}_{(1)}}=10,$
$p^{{}_{(2)}}/p^{{}_{(1)}}=10$ (blue curves); $(b)$ fixing the parameters
$\tau^{{}_{(2)}}/\tau^{{}_{(1)}}=3,$
$\bar{K}_{22}^{{}_{(2)}}/\bar{K}_{22}^{{}_{(1)}}=3,$
$p^{{}_{(2)}}/p^{{}_{(1)}}=3,$
$p^{{}_{(1)}}T_{0}/\left(\tau^{{}_{(1)}}\bar{K}_{22}^{{}_{(1)}}\right)=1,$ and
varying the thickness as $\eta=1$ (red curves), $\eta=10$ (blue curves).

$k_{2}(s)=\dfrac{1}{\mathrm{i}d_{2}}\left[\ln\left(\dfrac{-\tilde{I}_{1}^{{}_{j}}(s)\pm\sqrt{\left(\tilde{I}_{1}^{{}_{j}}(s)\right)^{2}-4}}{2}\right)-2\pi\mathrm{i}n\right],\quad\textrm{with
$n\in\mathbb{Z},\quad j={1,2,3},$}$ (1by)

where the details about the invariants $\tilde{I}_{1}^{{}_{j}}(s)$ are brought
in the Sec. B.3. of the Supplementary Material.

For the case in which the wave number $k_{1}$ is assumed zero and when
$s=\mathrm{i}\omega,$ the Fig. 5 exhibits the behaviour of the thermal wave
propagation given by Eqs. (1bw) and (1bx) for the homogenized material
represented in light blue and light red curves, while the thermal wave
function (1by) for the heterogeneous solid via Floquet-Bloch theory is
presented in dark blue and light red curves, where the parameters have been
set as $\eta=1,$ $\tau^{{}_{(2)}}/\tau^{{}_{(1)}}=3,$
$\bar{K}_{22}^{{}_{(2)}}/\bar{K}_{22}^{{}_{(1)}}=3,$
$p^{{}_{(2)}}/p^{{}_{(1)}}=3$ and
$p^{{}_{(1)}}T_{0}/\left(\tau^{{}_{(1)}}\bar{K}_{22}^{{}_{(1)}}\right)=1.$
Both light and dark blue lines stand for the translations of their respective
dispersion curves due to the periodicity of the material along
$\boldsymbol{e}_{2}$ around $\Re e(k_{2}\varepsilon)=2\pi n,$ which they go
all along the real wave number axis $\Re e(k_{2}\varepsilon).$ The dispersion
relations of the homogenized model can also be determined seeking solutions
from the governing equations (1ama) and (1amb) (with source terms assumed
zero) of the forms
${\hat{\boldsymbol{U}}}^{M}=\check{\hat{\boldsymbol{U}}}^{M}\exp[\mathrm{i}(\boldsymbol{k}\cdot\textbf{x}+2\pi
n)]$ and
${\hat{\Theta}}^{M}=\check{\hat{\Theta}}^{M}\exp[\mathrm{i}(\boldsymbol{k}\cdot\textbf{x}+2\pi
n)]$ (with $n\in\mathbb{Z}$) for the displacement and temperature,
respectively. When $n=0$ one arrives to the functions (1bs), (1bu), (1bv),
while when $n=-1,$ $n=1$ one produces the translated curves drawn in light
blue displayed in the Figs. 5$(a)$ and 5$(c).$ Besides, since a first order
homogenization was applied one may notice an accurate approximation between
the light curves and dark ones along the thermal wave propagation, in other
words the dispersion function derived from the homogenization process shows a
very good agreement with the dispersion function obtained taking advantage of
the Floquet–Bloch theory, seen in the Secs. B.1., B.2. and B.3. of the
Supplementary Material, for the interval $2\pi(1-n)/3<\Re
e(k_{2}\varepsilon)<2\pi(1+n)/3,\forall n\in\mathbb{Z}.$ The Figs. 5$(a)$ and
5$(c)$ show the dispersion functions associated to $n=-1,$ $n=0$ and $n=1,$
where may be observed that for a low frequency the curves of the heterogeneous
continuum (dark red) are quite precise in respect to those obtained from the
homogenization process of the material (light red). For a high frequency can
be noticed that the curves of the heterogeneous continuum (dark red) are well
approximated with a tiny deviation when compared to the translated curves
($n=-1$ and $n=1$) gotten from the homogenization process of the material
(light blue). Fig. 5$(b)$ shows a zoom of the thermal wave propagation by
attenuation in space seen in Fig. 5$(a).$ The dimensionless angular frequency
$\omega\tau^{{}_{(1)}}$ by the dimensionless real wave number $\Re
e(k_{2}\varepsilon)$ plane and the dimensionless angular frequency
$\omega\tau^{{}_{(1)}}$ by the dimensionless attenuation factor $\Im
m(k_{2}\varepsilon)$ plane are represented in 5$(c)$ and 5$(d),$ respectively.

Figure 7: Dimensionless dispersion functions associated to compressional-
thermal waves characterized when $k_{1}=0.$ Comparison between a first order
homogenized model (light blue curve and light red curve) with its respective
heterogeneous one (dark blue curve and dark red curve), assuming the numerical
values for the dimensionless parameters $\eta=1,$
$\tau^{{}_{(2)}}/\tau^{{}_{(1)}}=3,$
${C}_{2222}^{{}_{(2)}}/{C}_{2222}^{{}_{(1)}}=3,$
$\bar{K}_{22}^{{}_{(2)}}/\bar{K}_{22}^{{}_{(1)}}=3,$
$p^{{}_{(2)}}/p^{{}_{(1)}}=3,$ $\rho^{{}_{(2)}}/\rho^{{}_{(1)}}=2,$
$(\alpha_{22}^{{}_{(1)}}T_{0})/{C}_{2222}^{{}_{(1)}}=1/100,$
$\alpha_{22}^{{}_{(2)}}T_{0}/{C}_{2222}^{{}_{(2)}}=1/10,$
$\alpha_{22}^{{}_{(1)}}\eta\sqrt{C_{2222}^{{}_{(1)}}/\rho^{{}_{(1)}}}/\bar{K}_{22}^{{}_{(1)}}=1/100,$
$\alpha_{22}^{{}_{(2)}}\eta\sqrt{C_{2222}^{{}_{(1)}}/\rho^{{}_{(1)}}}/\bar{K}_{22}^{{}_{(2)}}=1/10,$
$p^{{}_{(1)}}T_{0}\eta\sqrt{C_{2222}^{{}_{(1)}}/\rho^{{}_{(1)}}}/\bar{K}_{22}^{{}_{(1)}}=1$
and
$\tau^{{}_{(1)}}\sqrt{C_{2222}^{{}_{(1)}}/\rho^{{}_{(1)}}}/\varepsilon=1/10.$
$(a)$ $\omega\varepsilon\sqrt{\rho^{{}_{(1)}}/C_{2222}^{{}_{(1)}}}$ vs. $\Re
e(k_{2}\varepsilon)\times\Im m(k_{2}\varepsilon);$ $(b)$ zoomed view of the
angular frequency spectrum
$\omega\varepsilon\sqrt{\rho^{{}_{(1)}}/C_{2222}^{{}_{(1)}}}$ vs. $\Re
e(k_{2}\varepsilon)\times\Im m(k_{2}\varepsilon);$ $(c)$ view of the plane
$\omega\varepsilon\sqrt{\rho^{{}_{(1)}}/C_{2222}^{{}_{(1)}}}\times\Re
e(k_{2}\varepsilon);$ $(d)$ view of the plane
$\omega\varepsilon\sqrt{\rho^{{}_{(1)}}/C_{2222}^{{}_{(1)}}}\times\Im
m(k_{2}\varepsilon).$

For a better understanding of the first order homogenization facing its
heterogeneous version, the Fig. 6 has been made varying some parameters, where
dark curves represent the heterogeneous material obtained via Floquet-Bloch
theory, whereas light curves are for the first order homogenized model.
Essentially, Fig. 6$(a)$ displays three zoomed different situations for the
thermal wave propagation when $\eta=1$ and
$p^{{}_{(1)}}T_{0}/\left(\tau^{{}_{(1)}}\bar{K}_{22}^{{}_{(1)}}\right)=1,$
namely the curves in red were plotted setting
$\tau^{{}_{(2)}}/\tau^{{}_{(1)}}=3,$
$\bar{K}_{22}^{{}_{(2)}}/\bar{K}_{22}^{{}_{(1)}}=3,$
$p^{{}_{(2)}}/p^{{}_{(1)}}=3,$ for the second scenario the curves in green
were generated choosing $\tau^{{}_{(2)}}/\tau^{{}_{(1)}}=5,$
$\bar{K}_{22}^{{}_{(2)}}/\bar{K}_{22}^{{}_{(1)}}=5,$
$p^{{}_{(2)}}/p^{{}_{(1)}}=5,$ and for the third situation the curves in blue
were plotted choosing the parameters $\tau^{{}_{(2)}}/\tau^{{}_{(1)}}=10,$
$\bar{K}_{22}^{{}_{(2)}}/\bar{K}_{22}^{{}_{(1)}}=10,$
$p^{{}_{(2)}}/p^{{}_{(1)}}=10.$ Furthermore, regarding the behaviour led by
the different values of the dimensionless parameters taken, it can be noted in
the Fig. 6$(a)$ that an increase in the parameters
$\tau^{{}_{(2)}}/\tau^{{}_{(1)}},$
$\bar{K}_{22}^{{}_{(2)}}/\bar{K}_{22}^{{}_{(1)}}$ and
$p^{{}_{(2)}}/p^{{}_{(1)}},$ produces curves that bend less sharply, and hence
have smaller curvatures. With an analogous idea, Fig. 6$(b)$ displays two
scenarios varying the thickness that are $\eta=1$ represented by the red
curves and $\eta=10$ represented by the blue curves, while the parameters
$\tau^{{}_{(2)}}/\tau^{{}_{(1)}}=3,$
$\bar{K}_{22}^{{}_{(2)}}/\bar{K}_{22}^{{}_{(1)}}=3,$
$p^{{}_{(2)}}/p^{{}_{(1)}}=3$ and
$p^{{}_{(1)}}T_{0}/\left(\tau^{{}_{(1)}}\bar{K}_{22}^{{}_{(1)}}\right)=1,$ are
fixed. Analysing the Fig. 6$(b),$ it is notable that a decrease in the
dimensionless thickness induce to bigger curvatures along the wave
propagation. From these comparisons between the wave for the heterogeneous
solid and their respectives for the homogenized model illustrated in Fig. 6,
it may be verified a good matching for the wave propagation between the two
models exemplified considering the interval $2\pi(1-n)/3<\Re
e(k_{2}\varepsilon)<2\pi(1+n)/3,\forall n\in\mathbb{Z}.$

Figure 8: Dimensionless dispersion functions associated to compressional-
thermal waves characterized when $k_{1}=0.$ Comparison between homogenized
models (light curves) and their respective heterogeneous ones (dark curves)
given different values of the parameters. In both the images the dimensionless
parameters are fixed as
$(\alpha_{22}^{{}_{(1)}}T_{0})/{C}_{2222}^{{}_{(1)}}=1/100,$
$\alpha_{22}^{{}_{(2)}}T_{0}/{C}_{2222}^{{}_{(2)}}=1/10,$
$\alpha_{22}^{{}_{(1)}}\eta\sqrt{C_{2222}^{{}_{(1)}}/\rho^{{}_{(1)}}}/\bar{K}_{22}^{{}_{(1)}}=1/100,$
$\alpha_{22}^{{}_{(2)}}\eta\sqrt{C_{2222}^{{}_{(1)}}/\rho^{{}_{(1)}}}/\bar{K}_{22}^{{}_{(2)}}=1/10,$
$p^{{}_{(1)}}T_{0}\eta\sqrt{C_{2222}^{{}_{(1)}}/\rho^{{}_{(1)}}}/\bar{K}_{22}^{{}_{(1)}}=1$
and
$\tau^{{}_{(1)}}\sqrt{C_{2222}^{{}_{(1)}}/\rho^{{}_{(1)}}}/\varepsilon=1/10.$
$(a)$ setting $\eta=1,$ and varying the parameters as
$\tau^{{}_{(2)}}/\tau^{{}_{(1)}}=3,$
${C}_{2222}^{{}_{(2)}}/{C}_{2222}^{{}_{(1)}}=3,$
$\bar{K}_{22}^{{}_{(2)}}/\bar{K}_{22}^{{}_{(1)}}=3,$
$p^{{}_{(2)}}/p^{{}_{(1)}}=3,$ $\rho^{{}_{(2)}}/\rho^{{}_{(1)}}=2,$ (red
curves), $\tau^{{}_{(2)}}/\tau^{{}_{(1)}}=5,$
${C}_{2222}^{{}_{(2)}}/{C}_{2222}^{{}_{(1)}}=5,$
$\bar{K}_{22}^{{}_{(2)}}/\bar{K}_{22}^{{}_{(1)}}=5,$
$p^{{}_{(2)}}/p^{{}_{(1)}}=5,$ $\rho^{{}_{(2)}}/\rho^{{}_{(1)}}=4,$ (green
curves), $\tau^{{}_{(2)}}/\tau^{{}_{(1)}}=10,$
${C}_{2222}^{{}_{(2)}}/{C}_{2222}^{{}_{(1)}}=10,$
$\bar{K}_{22}^{{}_{(2)}}/\bar{K}_{22}^{{}_{(1)}}=10,$
$p^{{}_{(2)}}/p^{{}_{(1)}}=10,$ $\rho^{{}_{(2)}}/\rho^{{}_{(1)}}=6,$ (blue
curves); $(b)$ fixing the parameters $\tau^{{}_{(2)}}/\tau^{{}_{(1)}}=3,$
${C}_{2222}^{{}_{(2)}}/{C}_{2222}^{{}_{(1)}}=3,$
$\bar{K}_{22}^{{}_{(2)}}/\bar{K}_{22}^{{}_{(1)}}=3,$
$p^{{}_{(2)}}/p^{{}_{(1)}}=3,$ $\rho^{{}_{(2)}}/\rho^{{}_{(1)}}=2,$ and
varying the thickness as $\eta=1$ (red curves), $\eta=20$ (blue curves).
Figure 9: Dimensionless dispersion functions associated to shear waves
characterized when $k_{1}=0.$ Comparison between a first order homogenized
model (light blue curve and light red curve) with its respective heterogeneous
one (dark blue curve and dark red curve), assuming the numerical values for
the dimensionless parameters $\eta=1,$
${C}_{1212}^{{}_{(2)}}/{C}_{1212}^{{}_{(1)}}=1$ and
$\rho^{{}_{(2)}}/\rho^{{}_{(1)}}=2.$ $(a)$
$\omega\varepsilon\sqrt{\rho^{{}_{(1)}}/C_{1212}^{{}_{(1)}}}$ vs. $\Re
e(k_{2}\varepsilon)\times\Im m(k_{2}\varepsilon);$ $(b)$ zoomed view of the
angular frequency spectrum
$\omega\varepsilon\sqrt{\rho^{{}_{(1)}}/C_{1212}^{{}_{(1)}}}$ vs. $\Re
e(k_{2}\varepsilon)\times\Im m(k_{2}\varepsilon);$ $(c)$ view of the plane
$\omega\varepsilon\sqrt{\rho^{{}_{(1)}}/C_{1212}^{{}_{(1)}}}\times\Re
e(k_{2}\varepsilon);$ $(d)$ view of the plane
$\omega\varepsilon\sqrt{\rho^{{}_{(1)}}/C_{1212}^{{}_{(1)}}}\times\Im
m(k_{2}\varepsilon).$ Figure 10: Dimensionless dispersion functions associated
to shear waves characterized when $k_{1}=0.$ Comparison between homogenized
models (light curves) and their respective heterogeneous ones (dark curves)
given different values of the parameters. $(a)$ setting $\eta=1,$ and varying
the parameters as ${C}_{1212}^{{}_{(2)}}/{C}_{1212}^{{}_{(1)}}=1$ and
$\rho^{{}_{(2)}}/\rho^{{}_{(1)}}=2,$ (red curves);
${C}_{1212}^{{}_{(2)}}/{C}_{1212}^{{}_{(1)}}=5$ and
$\rho^{{}_{(2)}}/\rho^{{}_{(1)}}=5,$ (green curves);
${C}_{1212}^{{}_{(2)}}/{C}_{1212}^{{}_{(1)}}=10$ and
$\rho^{{}_{(2)}}/\rho^{{}_{(1)}}=10,$ (blue curves); $(b)$ the non-dimensional
parameters are adopted as ${C}_{1212}^{{}_{(2)}}/{C}_{1212}^{{}_{(1)}}=1,$
$\rho^{{}_{(2)}}/\rho^{{}_{(1)}}=2,$ and varying the thickness as $\eta=1$
(red curves), $\eta=20$ (blue curves).

Likewise the previous comparative test, under the hypothesis when the wave
number $k_{1}$ is assumed zero and the complex frequency is taken as
$s=\mathrm{i}\omega,$ the compressional-thermal wave function of the
heterogeneous continuum via Floquet-Bloch theory and the compressional-thermal
waves of the first order homogenized material are illustrated in the Fig. 7 by
the dark blue and red lines, and the light blue and red lines, respectively.
The dimensionless parameters were chosen as $\eta=1,$
$\tau^{{}_{(2)}}/\tau^{{}_{(1)}}=3,$
${C}_{2222}^{{}_{(2)}}/{C}_{2222}^{{}_{(1)}}=3,$
$\bar{K}_{22}^{{}_{(2)}}/\bar{K}_{22}^{{}_{(1)}}=3,$
$p^{{}_{(2)}}/p^{{}_{(1)}}=3,$ $\rho^{{}_{(2)}}/\rho^{{}_{(1)}}=2,$
$(\alpha_{22}^{{}_{(1)}}T_{0})/{C}_{2222}^{{}_{(1)}}=1/100,$
$\alpha_{22}^{{}_{(2)}}T_{0}/{C}_{2222}^{{}_{(2)}}=1/10,$
$\alpha_{22}^{{}_{(1)}}\eta\sqrt{C_{2222}^{{}_{(1)}}/\rho^{{}_{(1)}}}/\bar{K}_{22}^{{}_{(1)}}=1/100,$
$\alpha_{22}^{{}_{(2)}}\eta\sqrt{C_{2222}^{{}_{(1)}}/\rho^{{}_{(1)}}}/\bar{K}_{22}^{{}_{(2)}}=1/10,$
$p^{{}_{(1)}}T_{0}\eta\sqrt{C_{2222}^{{}_{(1)}}/\rho^{{}_{(1)}}}/\bar{K}_{22}^{{}_{(1)}}=1$
and
$\tau^{{}_{(1)}}\sqrt{C_{2222}^{{}_{(1)}}/\rho^{{}_{(1)}}}/\varepsilon=1/10.$
In the Figs. 7$(a)$ and 7$(c),$ both light and dark blue lines stand for the
translations of their respective dispersion curves due to the periodicity of
the material along $\boldsymbol{e}_{2}$ around $\Re e(k_{2}\varepsilon)=2\pi
n,$ which they go all along the real wave number axis $\Re
e(k_{2}\varepsilon).$ Moreover, one may observe a precise estimation between
the light lines and dark ones along the compressional-thermal wave
propagation. This means that the dispersion function derived from the first
order homogenization process yield to a quite good matching when contrasted
against the dispersion function determined with the Floquet–Bloch theory for
the heterogeneous solid elucidated in the Secs. B.1., B.2. and B.3. of the
Supplementary Material, for the interval $2\pi(1-n)/3<\Re
e(k_{2}\varepsilon)<2\pi(1+n)/3,\forall n\in\mathbb{Z}.$ The Figs. 7$(a)$ and
7$(c),$ exhibit the dispersion functions when $n=-1,$ $n=0$ and $n=1,$ where
may be noted that for a low frequency the curves of the heterogeneous
continuum (dark red) are matched very well for those obtained from the
homogenization process of the material (light red), in addition the same
observation may be seen for the translated blue lines. Fig. 7$(b)$ shows a
zoom of the compressional-thermal wave propagation by spatial damping seen in
Fig. 7$(a).$ A bi-dimensional view of the dimensionless angular frequency
$\omega\varepsilon\sqrt{\rho^{{}_{(1)}}/C_{2222}^{{}_{(1)}}}$ by the
dimensionless real wave number $\Re e(k_{2}\varepsilon)$ plane and the
dimensionless angular frequency
$\omega\varepsilon\sqrt{\rho^{{}_{(1)}}/C_{2222}^{{}_{(1)}}}$ by the
dimensionless attenuation factor $\Im m(k_{2}\varepsilon)$ plane are
represented in 7$(c)$ and 7$(d),$ respectively. It is worthy of note that, the
branches of the spectrum characterized when the real part $\Re
e(k_{2}\varepsilon)$ is set fixed as $\Im m(k_{2}\varepsilon)$ varies, as
expected, which characterize the first frequency stop-band for quasi-
compressional waves and which propagate perpendicularly to the layering
direction, they are unable to be approximated by the frequency spectrum
obtained from the homogenized first order model, as showed in the Fig. 7$(d).$
In order to overcome such encumbrance, it is feasible to use a perturbative
technique, which will lead to a local asymptotic approximation of the
compressional-thermal wave through local explicit and closed-form parametric
expressions in the space of complex wave number and the angular frequency [83,
69].

The Fig. 8 compares some cases of the compressional-thermal wave assuming
different values of the dimensionless parameters, where the dark lines
represent the waves obtained for the heterogeneous material via Floquet-Bloch
theory, while the light lines are for the first order homogenized model. The
Figs. 8$(a)$ and 8$(b)$ were generated taking the following values for
dimensionless parameters
$(\alpha_{22}^{{}_{(1)}}T_{0})/{C}_{2222}^{{}_{(1)}}=1/100,$
$\alpha_{22}^{{}_{(2)}}T_{0}/{C}_{2222}^{{}_{(2)}}=1/10,$
$\alpha_{22}^{{}_{(1)}}\eta\sqrt{C_{2222}^{{}_{(1)}}/\rho^{{}_{(1)}}}/\bar{K}_{22}^{{}_{(1)}}=1/100,$
$\alpha_{22}^{{}_{(2)}}\eta\sqrt{C_{2222}^{{}_{(1)}}/\rho^{{}_{(1)}}}/\bar{K}_{22}^{{}_{(2)}}=1/10,$
$p^{{}_{(1)}}T_{0}\eta\sqrt{C_{2222}^{{}_{(1)}}/\rho^{{}_{(1)}}}/\bar{K}_{22}^{{}_{(1)}}=1$
and
$\tau^{{}_{(1)}}\sqrt{C_{2222}^{{}_{(1)}}/\rho^{{}_{(1)}}}/\varepsilon=1/10.$
Assuming $\eta=1,$ Fig. 8$(a)$ shows three zoomed different cases for the
compressional-thermal wave propagation. Therefore, the red lines were plotted
setting $\tau^{{}_{(2)}}/\tau^{{}_{(1)}}=3,$
${C}_{2222}^{{}_{(2)}}/{C}_{2222}^{{}_{(1)}}=3,$
$\bar{K}_{22}^{{}_{(2)}}/\bar{K}_{22}^{{}_{(1)}}=3,$
$p^{{}_{(2)}}/p^{{}_{(1)}}=3,$ $\rho^{{}_{(2)}}/\rho^{{}_{(1)}}=2,$ the green
lines were drawn choosing $\tau^{{}_{(2)}}/\tau^{{}_{(1)}}=5,$
${C}_{2222}^{{}_{(2)}}/{C}_{2222}^{{}_{(1)}}=5,$
$\bar{K}_{22}^{{}_{(2)}}/\bar{K}_{22}^{{}_{(1)}}=5,$
$p^{{}_{(2)}}/p^{{}_{(1)}}=5,$ $\rho^{{}_{(2)}}/\rho^{{}_{(1)}}=4,$ and
finally the lines in blue represent the situation in which
$\tau^{{}_{(2)}}/\tau^{{}_{(1)}}=10,$
${C}_{2222}^{{}_{(2)}}/{C}_{2222}^{{}_{(1)}}=10,$
$\bar{K}_{22}^{{}_{(2)}}/\bar{K}_{22}^{{}_{(1)}}=10,$
$p^{{}_{(2)}}/p^{{}_{(1)}}=10,$ $\rho^{{}_{(2)}}/\rho^{{}_{(1)}}=6.$ Taking a
look closely at the curvatures of the compressional-thermal waves plotted in
the Fig. 8$(a),$ it may be affirmed that increasing the numerical values of
$\tau^{{}_{(2)}}/\tau^{{}_{(1)}},$
${C}_{2222}^{{}_{(2)}}/{C}_{2222}^{{}_{(1)}},$
$\bar{K}_{22}^{{}_{(2)}}/\bar{K}_{22}^{{}_{(1)}},$ $p^{{}_{(2)}}/p^{{}_{(1)}}$
and $\rho^{{}_{(2)}}/\rho^{{}_{(1)}},$ the curves show a tendency of having
smaller curvatures. The compressional-thermal waves behaviour is foretold in
matters of different values of the thickness $\eta,$ the Fig. 8$(b)$ displays
two scenarios where is fixed $\tau^{{}_{(2)}}/\tau^{{}_{(1)}}=3,$
${C}_{2222}^{{}_{(2)}}/{C}_{2222}^{{}_{(1)}}=3,$
$\bar{K}_{22}^{{}_{(2)}}/\bar{K}_{22}^{{}_{(1)}}=3,$
$p^{{}_{(2)}}/p^{{}_{(1)}}=3,$ $\rho^{{}_{(2)}}/\rho^{{}_{(1)}}=2,$ and
varying the thickness that are $\eta=1$ and $\eta=20$ represented by the red
and blue lines, respectively. Such comparisons provided by the Fig. 8 shows a
good agreement of the compressional-thermal wave propagation between the two
models exemplified considering the interval $2\pi(1-n)/3<\Re
e(k_{2}\varepsilon)<2\pi(1+n)/3,\forall n\in\mathbb{Z}.$

Some plots of the shear waves may be seen in Fig. 9, where red lines stand for
the central shear wave, blue lines for the translated shear waves, and the
nuances in the color distinguish the shear curve of the heterogeneous
continuum (dark lines) from the shear wave first order homogenized (light
lines). In regards to de capability of the first order homogenization process
for the shear wave, it also shows itself in a good agreement with its
heterogeneous continuum from the Floquet-Bloch theory. The Fig. 9 was drawn
adopting the dimensionless parameters as $\eta=1,$
${C}_{1212}^{{}_{(2)}}/{C}_{1212}^{{}_{(1)}}=1$ and
$\rho^{{}_{(2)}}/\rho^{{}_{(1)}}=2.$ Note that, as verified for the
compressional-thermal wave, and with a cleaner comprehension looking at the
Fig. 9$(d),$ also the frequency stop-bands of the shear wave propagating
perpendicularly to the layering direction is not well described by the first
order homogenization technique developed. As aforementioned, to remedy this
issue the shear waves may be approached by a perturbative method through
explicit and closed-form parametric expressions, in which will approximate
locally the real spectrum of the non-homogeneous material [83, 69].

The Fig. 10$(a)$ exhibits three zoomed different cases for the shear wave
propagation for a fixed thickness $\eta=1.$ The red lines represent the
choices for the parameters as ${C}_{1212}^{{}_{(2)}}/{C}_{1212}^{{}_{(1)}}=1,$
$\rho^{{}_{(2)}}/\rho^{{}_{(1)}}=2,$ the green lines stand for
${C}_{1212}^{{}_{(2)}}/{C}_{1212}^{{}_{(1)}}=5,$
$\rho^{{}_{(2)}}/\rho^{{}_{(1)}}=5,$ and lastly, the lines that approximate
locally the real spectrum of the heterogeneous material in blue correspond to
the case when ${C}_{1212}^{{}_{(2)}}/{C}_{1212}^{{}_{(1)}}=10,$
$\rho^{{}_{(2)}}/\rho^{{}_{(1)}}=10.$ From a quick analysis over the Fig.
10$(a),$ it can be assured that the curvatures of the shear waves are reduced
numerically due to a rise in the values of the non-dimensional parameters
${C}_{1212}^{{}_{(2)}}/{C}_{1212}^{{}_{(1)}}$ and
$\rho^{{}_{(2)}}/\rho^{{}_{(1)}}.$ Varying the thickness $\eta$ namely
$\eta=1$ (red curves) and $\eta=20$ (blue curves), and adopting the
dimensionless parameters as ${C}_{1212}^{{}_{(2)}}/{C}_{1212}^{{}_{(1)}}=1,$
$\rho^{{}_{(2)}}/\rho^{{}_{(1)}}=2,$ both shear waves for these assumptions
are contrasted in the Fig. 10$(b).$ Essentially, the Fig. 10 illustrates the
good estimation of the shear wave propagation between the two models developed
considering the interval $2\pi(1-n)/3<\Re
e(k_{2}\varepsilon)<2\pi(1+n)/3,\forall n\in\mathbb{Z}.$

## Conclusions

The present work has formulated an asymptotic homogenization approach for
describing composites that have a periodic microstructure characteristic in
presence of the generalized thermoelasticity theory with a single periodic
spatially dependent relaxation time [15]. The down-scaling relations herein
developed permit such technique to express the micro-fields namely,
displacement and temperature, in terms of the macroscopic fields and their
derivatives, through the $\mathcal{Q}$-periodic perturbation functions. The
cell problems i.e., a cascade of inhomogeneous recursive differential problems
defined on the unit cell $\mathcal{Q}$ and having zero mean values over this
unit cell, give the availability to seek for those perturbation functions.

The cell problems were derived substituting series expansions for both of the
micro-fields in powers of the microstructural characteristic size
$\varepsilon$ into the micro governing equations and rearranging the outcomes
as asymptotic expansions in terms of $\varepsilon$ exponents. Inserting the
down-scaling relations into the microscopic thermoelastic field equations has
led to the average field equations of infinite order, which are asymptotically
equivalent to the initial governing equations of the heterogeneous continuum,
in which the formal solutions might be derived by plugging an asymptotic
infinite series of the macro fields in powers of the microstructural size
$\varepsilon.$ Particularly, truncating such series of the macroscopic
displacement and temperature up to the zeroth order yield to the governing
field equations of the homogeneous first order (Cauchy) thermoelastic medium
equivalent to the heterogeneous one, which the elastic, thermal dilatation and
heat conduction overall constitutive tensors have been analytically stated. In
particular, the overall heat conduction tensor depends on the thermal
conductivity tensor and the relaxation time of the phases. Consequently, the
higher the order of the truncated asymptotic expansions of the macro fields
the better the estimation for solutions of the heterogeneous problem.

The reliability of such analysis has been evaluated through a benchmark test
between the presented first order homogenized model and the heterogeneous
continuum over an illustrative example of a bi-dimensional two-phase periodic
layered material endowed with an orthotropy axis parallel to the layering
direction. In what regards the investigation of the dispersion properties of
the heterogeneous thermoelastic medium the Floquet-Bloch theory has been
implemented. Such path has led to a sixth order eigenvector-eigenvalue problem
solved over the periodic cell subjected to Floquet-Bloch periodic boundary
conditions, which throughout the comparative analysis were considered the
uncoupled and coupled scenarios. The eigenproblem provided an imaginary
implicit function and a real implicit function of the complex frequency and
wave vector, where the intersection between them defines the frequency
spectrum. In both eigenvector-eigenvalue problems from the Floquet-Bloch
approach, characterizing the heterogeneous material, the thermal waves from
the uncoupled problem as well as the compressional-thermal waves from the
coupled problem, and the shear waves were compared against their respectives
when applied the first order homogenized approximation.

Over the comparisons carried out in this work, it has been achieved a very
good agreement between the dispersion curves computed from the two different
methods i.e., thermal, compressional-thermal and shear waves along and normal
direction to the layering. In essence, such strong connection observed between
the methods confirms that the multi-field homogenization technique of the
first order herein proposed has been shown itself to be a quite good tool to
estimate the macroscopical elastic and thermal overall properties of the
equivalent heterogeneous body having periodic microstructure under
thermoelastic phenomena with sufficient accuracy.

## Acknowledgments

The authors gratefully acknowledge financial support from National Group of
Mathematical Physics (GNFM-INdAM). AB would like to acknowledge financial
support from Compagnia San Paolo, project MINIERA no. I34I20000380007 and from
University of Trento, project UNMASKED 2020. MP would like to acknowledge
financial support from the Italian Ministry of Education, University and
Research (MIUR) to the research project of relevant national interest (PRIN
2017) “XFAST-SIMS: Extra-fast and accurate simulation of complex structural
systems” (CUP: D68D19001260001).

## References

  * Biot [1956] Maurice Anthony Biot.  Thermoelasticity and irreversible thermodynamics.  _Journal of applied physics_ , 27(3):240–253, 1956.
  * Kupradze et al. [1980] VD Kupradze, TG Gegelia, Mo Oo Basheleishvili, To V Burchuladze, and E Sternberg.  Three-dimensional problems of the mathematical theory of elasticity and thermoelasticity.  1980\.
  * Carcione et al. [2019] José M Carcione, Zhi-Wei Wang, Wenchang Ling, Ettore Salusti, Jing Ba, and Li-Yun Fu.  Simulation of wave propagation in linear thermoelastic media.  _Geophysics_ , 84(1):T1–T11, 2019.
  * Joseph and Preziosi [1989] Daniel D Joseph and Luigi Preziosi.  Heat waves.  _Reviews of Modern Physics_ , 61(1):41, 1989.
  * Liu et al. [1999] Jing Liu, Xu Chen, and Lisa X Xu.  New thermal wave aspects on burn evaluation of skin subjected to instantaneous heating.  _IEEE transactions on biomedical engineering_ , 46(4):420–428, 1999.
  * Ignaczak and Ostoja-Starzewski [2010] Józef Ignaczak and Martin Ostoja-Starzewski.  _Thermoelasticity with finite wave speeds_.  Oxford University Press, 2010.
  * Chandrasekharaiah [1998] DS Chandrasekharaiah.  Hyperbolic thermoelasticity: a review of recent literature.  1998\.
  * Zamani et al. [2011] Arash Zamani, Richard B Hetnarski, and M Reza Eslami.  Second sound in a cracked layer based on lord–shulman theory.  _Journal of Thermal Stresses_ , 34(3):181–200, 2011.
  * Kaminski [1990] W Kaminski.  Hyperbolic heat conduction equation for materials with a nonhomogeneous inner structure.  _Journal of Heat Transfer_ , 112(3):555–560, 1990\.
  * Mitra et al. [1995] K Mitra, S Kumar, A Vedevarz, and MK Moallemi.  Experimental evidence of hyperbolic heat conduction in processed meat.  _Journal of Heat Transfer_ , 117(3):568–573, 1995\.
  * Narayanamurti and Dynes [1972] V Narayanamurti and RC Dynes.  Observation of second sound in bismuth.  _Physical Review Letters_ , 28(22):1461, 1972\.
  * Hetnarski and Ignaczak [1999] Richard B Hetnarski and Jozef Ignaczak.  Generalized thermoelasticity.  _Journal of Thermal Stresses_ , 22(4-5):451–476, 1999.
  * Hetnarski et al. [2009] Richard B Hetnarski, M Reza Eslami, and GML Gladwell.  _Thermal stresses: advanced theory and applications_ , volume 41.  Springer, 2009.
  * Sherief et al. [2010] Hany H Sherief, AMA El-Sayed, and AM Abd El-Latief.  Fractional order theory of thermoelasticity.  _International Journal of Solids and structures_ , 47(2):269–275, 2010.
  * Lord and Shulman [1967] Harold Wesley Lord and Y Shulman.  A generalized dynamical theory of thermoelasticity.  _Journal of the Mechanics and Physics of Solids_ , 15(5):299–309, 1967.
  * Maxwell [1867] James Clerk Maxwell.  On the dynamical theory of gases.  _Philosophical transactions of the Royal Society of London, 157_ , pages 49–88, 1867.
  * Cattaneo [1948] Carlo Cattaneo.  Sulla conduzione del calore.  _Atti Sem. Mat. Fis. Univ. Modena_ , 3:83–101, 1948.
  * Dhaliwal and Sherief [1980] Ranjit S Dhaliwal and Hani H Sherief.  Generalized thermoelasticity for anisotropic media.  _Quarterly of Applied Mathematics_ , 38(1):1–8, 1980.
  * Prunty [1978] J Prunty.  Dimensionally stable graphite composites for spacecraft structures.  _SAMPE Quarterly_ , 9(2):41–51, 1978.
  * Dodson and Rule [1989] Kelly J Dodson and John E Rule.  Thermal stability considerations for space flight optical benches.  _sampe_ , pages 1578–1589, 1989.
  * Farmer et al. [1992] Jeffery T Farmer, Deborah M Wahls, Robert L Wright, and A Louis Tahernia.  Thermal distortion analysis of an antenna-support truss in geosynchronous orbit.  _Journal of Spacecraft and Rockets_ , 29(3):386–393, 1992.
  * Kim et al. [2000] Ran Y Kim, Allan S Crasto, and Gregory A Schoeppner.  Dimensional stability of composite in a space thermal environment.  _Composites science and technology_ , 60(12-13):2601–2608, 2000.
  * Zhu and Sun [2003] RP Zhu and CT Sun.  Effects of fiber orientation and elastic constants on coefficients of thermal expansion in laminates.  _Mechanics of Advanced Materials and Structures_ , 10(2):99–107, 2003.
  * Nayfeh and Nemat-Nasser [1971] Adnan Nayfeh and S Nemat-Nasser.  Thermoelastic waves in solids with thermal relaxation.  _Acta Mechanica_ , 12(1-2):53–69, 1971.
  * Hawwa and Nayfeh [1995] Muhammad A Hawwa and Adnan H Nayfeh.  The general problem of thermoelastic waves in anisotropic periodically laminated composites.  _Composites Engineering_ , 5(12):1499–1517, 1995\.
  * Sanchez Palencia [1974] E Sanchez Palencia.  Comportements local et macroscopique d’un type de milieux physiques heterogenes.  _International Journal of Engineering Science_ , 12(4):331–351, 1974.
  * Papanicolau et al. [1978] G Papanicolau, A Bensoussan, and J-L Lions.  _Asymptotic analysis for periodic structures_ , volume 5.  Elsevier, 1978.
  * Bakhvalov and Panasenko [1984] NS Bakhvalov and GP Panasenko.  _Homogenization: Averaging Processes in Periodic Media_.  Kluwer Academic Publishers, Dordrecht-Boston-London, 1984.
  * Smyshlyaev and Cherednichenko [2000] Valery P Smyshlyaev and KD Cherednichenko.  On rigorous derivation of strain gradient effects in the overall behaviour of periodic heterogeneous media.  _Journal of the Mechanics and Physics of Solids_ , 48(6):1325–1357, 2000.
  * Bacigalupo and Gambarotta [2014a] Andrea Bacigalupo and Luigi Gambarotta.  Computational dynamic homogenization for the analysis of dispersive waves in layered rock masses with periodic ures.  _Computers and Geotechnics_ , 56:61–68, 2014a.
  * Nassar et al. [2016] Hussein Nassar, Q-C He, and Nicolas Auffray.  On asymptotic elastodynamic homogenization approaches for periodic media.  _Journal of the Mechanics and Physics of Solids_ , 88:274–290, 2016.
  * Gambin and Kroner [1989] B. Gambin and E. Kroner.  Higher-order terms in the homogenized stress-strain relation of periodic elastic media.  _physica status solidi (b)_ , 151(2):513–519, 1989.
  * Andrianov et al. [2008] Igor V Andrianov, Vladimir I Bolshakov, Vladyslav V Danishevs’ kyy, and Dieter Weichert.  Higher order asymptotic homogenization and wave propagation in periodic composite materials.  _Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences_ , 464(2093):1181–1201, 2008.
  * Bacigalupo and Gambarotta [2010] Andrea Bacigalupo and Luigi Gambarotta.  Second-order computational homogenization of heterogeneous materials with periodic microstructure.  _ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik_ , 90(10-11):796–811, 2010.
  * Yang et al. [2016] Zhiqiang Yang, Junzhi Cui, Yi Sun, Jun Liang, and Zihao Yang.  Multiscale analysis method for thermo-mechanical performance of periodic porous materials with interior surface radiation.  _International Journal for Numerical Methods in Engineering_ , 105(5):323–350, 2016.
  * Fish et al. [2019] Jacob Fish, Zhiqiang Yang, and Zifeng Yuan.  A second-order reduced asymptotic homogenization approach for nonlinear periodic heterogeneous materials.  _International Journal for Numerical Methods in Engineering_ , 119(6):469–489, 2019.
  * Willis [1981] John R Willis.  Variational and related methods for the overall properties of composites.  In _Advances in applied mechanics_ , volume 21, pages 1–78. Elsevier, 1981.
  * Smyshlyaev [2009] Valery P Smyshlyaev.  Propagation and localization of elastic waves in highly anisotropic periodic composites via two-scale homogenization.  _Mechanics of Materials_ , R59:434–447, 2009.
  * Bacigalupo and Gambarotta [2014b] Andrea Bacigalupo and Luigi Gambarotta.  Second-gradient homogenized model for wave propagation in heterogeneous periodic media.  _International Journal of Solids and Structures_ , 51(5):1052–1065, 2014b.
  * Kamotski and Smyshlyaev [2019] IIia V Kamotski and Valery P Smyshlyaev.  Bandgaps in two-dimensional high-contrast periodic elastic beam lattice materials.  _Journal of the Mechanics and Physics of Solids_ , 123:292–304, 2019.
  * Bigoni and Drugan [2007] D Bigoni and WJ Drugan.  Analytical derivation of cosserat moduli via homogenization of heterogeneous elastic materials.  2007\.
  * Milton and Willis [2007] Graeme W Milton and John R Willis.  On modifications of newton’s second law and linear continuum elastodynamics.  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 463(2079):855–880, 2007.
  * Bacca et al. [2013a] Mattia Bacca, Francesco Dal Corso, Daniele Veber, and Davide Bigoni.  Mindlin second-gradient elastic properties from dilute two-phase cauchy-elastic composites part i: closed form expression for the effective higher-order constitutive tensor.  _International Journal of Solids and Structures_ , 50(24):4010–4019, 2013a.
  * Bacca et al. [2013b] Mattia Bacca, Francesco Dal Corso, Daniele Veber, and Davide Bigoni.  Mindlin second-gradient elastic properties from dilute two-phase cauchy-elastic composites part ii: Higher-order constitutive properties and application cases.  _International Journal of Solids and Structures_ , 50(24):4020–4029, 2013b.
  * Bacigalupo et al. [2018] Andrea Bacigalupo, Marco Paggi, F Dal Corso, and D Bigoni.  Identification of higher-order continua equivalent to a cauchy elastic composite.  _Mechanics Research Communications_ , 93:11–22, 2018.
  * Ostoja-Starzewski et al. [1999] Martin Ostoja-Starzewski, Stéphane D Boccara, and Iwona Jasiuk.  Couple-stress moduli and characteristic length of a two-phase composite.  _Mechanics Research Communications_ , 26(4):387–396, 1999.
  * Kouznetsova et al. [2002] Varvara Kouznetsova, Marc GD Geers, and WA Marcel Brekelmans.  Multi-scale constitutive modelling of heterogeneous materials with a gradient-enhanced computational homogenization scheme.  _International journal for numerical methods in engineering_ , 54(8):1235–1260, 2002.
  * Forest [2002] Samuel Forest.  Homogenization methods and mechanics of generalized continua-part 2.  _Theoretical and applied mechanics_ , (28-29):113–144, 2002.
  * Lew et al. [2004] TL Lew, F Scarpa, and K Worden.  Homogenisation metamodelling of perforated plates.  _Strain_ , 40(3):103–112, 2004.
  * Scarpa et al. [2009] Fabrizio Scarpa, Sondipon Adhikari, and A Srikantha Phani.  Effective elastic mechanical properties of single layer graphene sheets.  _Nanotechnology_ , 20(6):065709, 2009.
  * De Bellis and Addessi [2011] Maria Laura De Bellis and Daniela Addessi.  A cosserat based multi-scale model for masonry structures.  _International Journal for Multiscale Computational Engineering_ , 9(5):543, 2011.
  * Forest and Trinh [2011] Samuel Forest and Duy Khanh Trinh.  Generalized continua and non-homogeneous boundary conditions in homogenisation methods.  _ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift fur Angewandte Mathematik und Mechanik_ , 91(2):90–109, 2011\.
  * Chen et al. [2013] YJ Chen, F Scarpa, YJ Liu, and JS Leng.  Elasticity of anti-tetrachiral anisotropic lattices.  _International Journal of Solids and Structures_ , 50(6):996–1004, 2013.
  * Monchiet et al. [2020] Vincent Monchiet, Nicolas Auffray, and Julien Yvonnet.  Strain-gradient homogenization: a bridge between the asymptotic expansion and quadratic boundary condition methods.  _Mechanics of Materials_ , 143:103309, 2020.
  * Temizer and Wriggers [2011] I Temizer and Peter Wriggers.  Homogenization in finite thermoelasticity.  _Journal of the Mechanics and Physics of Solids_ , 59(2):344–372, 2011.
  * L’Hostis and Devries [1998] Gildas L’Hostis and François Devries.  Characterization of the thermoelastic behavior of syntactic foams.  _Composites Part B: Engineering_ , 29(4):351–361, 1998.
  * Yu and Tang [2007] Wenbin Yu and Tian Tang.  A variational asymptotic micromechanics model for predicting thermoelastic properties of heterogeneous materials.  _International Journal of Solids and Structures_ , 44(22-23):7510–7525, 2007.
  * Pettermann and Suresh [2000] Heinz E Pettermann and Subra Suresh.  A comprehensive unit cell model: a study of coupled effects in piezoelectric 1–3 composites.  _International Journal of Solids and Structures_ , 37(39):5447–5464, 2000.
  * Berger et al. [2005] Harald Berger, Sreedhar Kari, Ulrich Gabbert, Reinaldo Rodríguez-Ramos, Julián Bravo-Castillero, and Raúl Guinovart-Díaz.  A comprehensive numerical homogenisation technique for calculating effective coefficients of uniaxial piezoelectric fibre composites.  _Materials Science and Engineering: A_ , 412(1-2):53–60, 2005.
  * Schröder and Keip [2012] Jörg Schröder and Marc-André Keip.  Two-scale homogenization of electromechanically coupled boundary value problems.  _Computational mechanics_ , 50(2):229–244, 2012\.
  * Zäh and Miehe [2013] Dominic Zäh and Christian Miehe.  Computational homogenization in dissipative electro-mechanics of functional materials.  _Computer Methods in Applied Mechanics and Engineering_ , 267:487–510, 2013.
  * Fantoni et al. [2017] Francesca Fantoni, Andrea Bacigalupo, and Marco Paggi.  Multi-field asymptotic homogenization of thermo-piezoelectric materials with periodic microstructure.  _International Journal of Solids and Structures_ , 120:31–56, 2017.
  * Ohno et al. [2000] N Ohno, X Wu, and T Matsuda.  Homogenized properties of elastic–viscoplastic composites with periodic internal structures.  _International Journal of Mechanical Sciences_ , 42(8):1519–1536, 2000.
  * Abdessamad et al. [2009] Zouhair Abdessamad, Ilya Kostin, Grigory Panasenko, and Valery P Smyshlyaev.  Memory effect in homogenization of a viscoelastic kelvin–voigt model with time-dependent coefficients.  _Mathematical Models and Methods in Applied Sciences_ , 19(09):1603–1630, 2009.
  * Chen et al. [2017] Qiang Chen, Guannan Wang, Xuefeng Chen, and Jia Geng.  Finite-volume homogenization of elastic/viscoelastic periodic materials.  _Composite Structures_ , 182:457–470, 2017.
  * Del Toro et al. [2019] Rosaria Del Toro, Andrea Bacigalupo, and Marco Paggi.  Characterization of wave propagation in periodic viscoelastic materials via asymptotic-variational homogenization.  _International Journal of Solids and Structures_ , 172:110–146, 2019.
  * Salvadori et al. [2014] A Salvadori, E Bosco, and D Grazioli.  A computational homogenization approach for li-ion battery cells: Part 1–formulation.  _Journal of the Mechanics and Physics of Solids_ , 65:114–137, 2014.
  * Bacigalupo et al. [2016] Andrea Bacigalupo, Lorenzo Morini, and Andrea Piccolroaz.  Multiscale asymptotic homogenization analysis of thermo-diffusive composite materials.  _International Journal of Solids and Structures_ , 85:15–33, 2016.
  * Fantoni and Bacigalupo [2020] Francesca Fantoni and Andrea Bacigalupo.  Wave propagation modeling in periodic elasto-thermo-diffusive materials via multifield asymptotic homogenization.  _International Journal of Solids and Structures_ , 196-197:99–128, 2020.
  * Aboudi et al. [2001] J Aboudi, M-J Pindera, and SM Arnold.  Linear thermoelastic higher-order theory for periodic multiphase materials.  _Journal of Applied Mechanics_ , 68(5):697–707, 2001.
  * Kanouté et al. [2009] P Kanouté, DP Boso, JL Chaboche, and BA Schrefler.  Multiscale methods for composites: a review.  _Archives of Computational Methods in Engineering_ , 16(1):31–75, 2009.
  * Zhang et al. [2007] HW Zhang, Sheng Zhang, Jin Ying Bi, and BA Schrefler.  Thermo-mechanical analysis of periodic multiphase materials by a multiscale asymptotic homogenization approach.  _International Journal for Numerical Methods in Engineering_ , 69(1):87–113, 2007.
  * Sixto-Camacho et al. [2013] Lázaro M Sixto-Camacho, Julian Bravo-Castillero, Renald Brenner, Raúl Guinovart-Díaz, Houari Mechkour, Reinaldo Rodríguez-Ramos, and Federico J Sabina.  Asymptotic homogenization of periodic thermo-magneto-electro-elastic heterogeneous media.  _Computers & Mathematics with Applications_, 66(10):2056–2074, 2013.
  * Floquet [1883] G. Floquet.  Sur les équations différentielles linéaires à coefficients périodiques.  _Annales de l’École Normale Supérieure_ , 12:47–88, 1883\.