1 Introduction

Although real elastic media are inhomogeneous, the classical theory of linear elasticity is always established on the assumption of homogeneity [1]. However, the effort on seeking the appropriate equations for inhomogeneous media never stops, especially when we are interested in the detailed information of the wave propagation in composites [2] or materials with micro-defects [3]. However, the corresponding theories develop very slowly.

Forty years ago, Willis established new linear elastodynamic equations for inhomogeneous media based on the variational theory and the perturbation method [4]:

$$ \left\langle {{\varvec{\sigma}}} \right\rangle = {\varvec{C}}^{{{\rm eff}}} * \left\langle {e} \right\rangle + {\varvec{S}}^{{{\text{eff}}}} \circ \left\langle {{\dot{\varvec{u}}}} \right\rangle , $$
(1)
$$\nabla \cdot \left\langle \varvec\sigma \right\rangle + {{f}} = \widehat{{{S}}}^{\text{eff}}\circ \left\langle {\mathop {{e}}\limits^{.} } \right\rangle + \rho ^{{{\text{eff}}}} \odot \left\langle {\mathop {{u}}\limits^{{..}} } \right\rangle,$$
(2)

where σ is the stress; u is the displacement; \({{e}} = \left( {{{u}}\nabla + \nabla {{u}}} \right)/2\) is the linear strain; f is the body force; the overhead dot denotes the differentiation with respect to time \(t; \left\langle {} \right\rangle\) denotes the ensemble average operation; \({{C}}^{{{\text{eff}}}} *\), \({{S}}^{{{\text{eff}}}} \circ\), \({\hat{{S}}}^{{{\text{eff}}}} \circ\), and \({{{\rho}}}^{{{\text{eff}}}} \odot\) are non-local operators that involve temporal-spatial convolutions. The derivation of these equations is briefly outlined in Appendix A according to reference [2]. However, researchers had not paid so much attention to Eqs. 1 and 2, until they tried to design metamaterials [5, 6] by using transformation methods.

Metamaterials are man-made inhomogeneous media that can control wave propagation in the desired manner and achieve many exotic phenomena, such as cloaking [7], negative refraction [8], etc. An important problem in designing metamaterials is how to determine the distribution of material properties independent of the wave incident direction. To solve this problem, Pendry et al. [9] and Leonhardt [10] proposed to use the transformation method [11], by which a wave equation in the virtual homogeneous configuration Bv is transformed into the physical inhomogeneous configuration B0 through local mappings of coordinates and field variables (see Fig. 1). If the wave equation does not change its form after this local spatial transformation, the distribution of the material properties can be obtained by comparing the corresponding terms in the original and transformed equations. Encouraged by the great success of the transformation methods in optics [12] and acoustics [13], researchers tried to implement it to elasticity. However, Milton et al. [14] found out that the classical linear elastodynamic equations are not form-invariant after the local spatial transformation in frequency ___domain. At the same time, it is surprising that the transformed new equations have the similar forms as Eqs. 1 and 2. This discovery manifests the importance of Willis’ original work, and then Eqs. 1 and 2 are customarily called the Willis equations. Since we will discuss other “Willis equations” later, Eqs. 1 and 2 are called the Classical Willis Equations (CWEs) in this paper.

Fig. 1
figure 1

The initial physical configuration B0 is resulted from the transformation from the virtual configuration Bv. Then, B0 changes to the deformed configuration B1 through a small perturbation

Full size image

Since the CWEs play non-fungible roles to design elastic metamaterials, more and more researchers are interested in investigating their properties [15,16,17,18,19,20,21,22,23] and even finding related applications [24]. However, there are still three questions lingering in their minds:

  1. (1)

    The first question is why the CWEs use velocity and strain rate coupling terms to describe lossless material (viscoelastic effects are not considered in the CWEs)? One expedient is regarding \({{S}}^{{{\text{eff}}}}\) as a complex-valued term [15,16,17,18]. Alternatively, one may suggest a coupling term related to accelerations [19,20,21].

  2. (2)

    Milton et al. [14] discuss the form-invariance under local spatial transformations in frequency ___domain. So, the second question is: what is the corresponding transformation in time ___domain that results in the CWEs? Actually, applying pure local spatial transformations to the classical elastodynamic equations in time ___domain, we can only obtain the following equations with displacement and strain coupling terms [25]:

    $${{\varvec{\sigma}}} = {{C}}:{{e}} + {{S}}^{0} \cdot {{u}},$$
    (3)
    $$\nabla \cdot \varvec \sigma + {{f}} = \left( {{{S}}^{0} } \right)^{{\text{T}}} :{{e}} + {{K}}^{0} \cdot {{u}} + \rho \cdot {{\ddot{{u}}}}.$$
    (4)

If transforming the above equations into frequency ___domain, we can get the same results of [14]. It should be noted that all material properties and field variables in Eqs. 3 and 4 are homogenized values of real media. Since the homogenization involves the average over a microstructure with many sub-components [26], which could move in different directions independently, we can generally obtain an anisotropic density ρ represented in tensorial form [27, 28]. This also conforms to the idea of the Jacobi metric for a general kinetic energy [29]. We have proved that S0 and K0 are related to the pre-stress σ0 in the initial configuration B0 [30]. As Appendix B shows, this can also be understood from energy view point. Since Eqs. 3 and 4 have the similar forms as those of Eqs. 1 and 2, they are called the Willis-Form Equations (WFEs) in this paper. It has been proved that the WFEs can meet the requirement of time synchronization after applying only local spatial transformations [31]. The static version of the WFEs is not only verified by the experiment of rotational springs [32], but also implemented to identify the distribution of residual stresses [33] and even helpful to the buckling analysis of thin-walled shells [34, 35]. Although we have showed that the WFEs and the CWEs are mathematically equivalent regarding the differential property of time convolution operator [32], the underlying physical mechanism is not very clear.

  1. (3)

    The third question is why the CWEs involve temporal-spatial convolution operators [2, 4,5,6], while the “Willis equations” derived by other researchers [15,16,17,18,19,20,21,22,23,24,25, 30,31,32,33,34,35] use only algebraic operators?

The above questions imply that the Willis equations may present in different forms. This can be easily understood by reviewing Appendix A, in which Eqs. 48 and 49 just come from a single Eq. 44. Willis has realized this non-uniqueness problem and tried to solve it by introducing eigen-strains [36]. He also demonstrated that the same mean stress and momentum could be obtained from equations in different forms [37]. However, it is necessary to clarify the underlying mechanism through a more general approach.

Actually, the form-invariance under transformations is not only the basis of the transformation method to design metamaterials, but also a general requirement for all physical laws. This is because coordinates are just man-made tools to facilitate calculating or measuring [38]. Therefore, the law of nature should be independent of coordinates, and we have much freedom to choose a coordinate system. It is well known that Newtonian mechanics are form-invariant under Galilean transformation, Maxwell equations are form-invariant under Lorentz transformation that coincides with the special theory of relativity. In continuum mechanics, we have the counterpart of the principle of material frame indifference or objectivity [39, 40], which requires the form-invariance under global Euclidean transformations or mappings of rigid rotations and translations. The underlying mechanism of these invariances is that a correct theory should present in the same form in any reference frame or look the same from any view point, i.e., the equation is form-invariant after coordinate transformations. This is always true not only for the aforementioned global transformations independent of coordinates, but also for local transformations that change at different temporal or spatial positions. A famous example of the form-invariance under local transformations is the general theory of relativity, which has local Lorentz invariance [41,42,43]. In 1918 Weyl even tried to unify electromagnetism and gravitation by investigating the invariance under the scale change of the metric tensor [44], and thus he introduced the concept of gauge transformation [41,42,43,44,45,46]. However, this attempt was not successful until it was adapted to the phase transformation of the wave function in quantum mechanics since 1928 [44]. The invariance of global phase transformation corresponds to the conservation of electric charge [46] and the invariance of local phase transformation coincides with electromagnetic interactions [43,44,45]. The principle of gauge invariance was further extended into the Yang-Mills theory [47] and the related Standard Model of particle physics by introducing more sophisticated transformations to describe the interactions inside atoms. This principle has become a dominant concept in theoretical physics in the second half of the twentieth century [41,42,43,44,45,46].

Almost at the dawning of the Yang–Mills theory, Kondo [48], Bilby [49] and Kröner [50] initiated the pioneering work of connecting local gauge transformations with various deformation incompatibilities due to defects [51]. After years of developing, it has become a new branch of theory called geometric continuum mechanics [3, 52,53,54,55,56]. A related concept is the principle of generalized objectivity, which has been proposed to describe the form-invariance of Lagrangian under local rigid body transformations [57]. Based on this idea, a recent work studied the dynamic crack branching in brittle materials [58]. Since this kind of research is closely related to advanced differential geometry theory, strong mathematical background is required to grasp its essence. This could be a daunting challenge for many interested engineers and material scientists.

Regarding the form-invariance under transformation, we should mention the Nöether theorems [42, 46, 59,60,61,62] published in 1918, which prove the equivalence among the invariance, symmetry and conservation laws. The Nöether’s first theorem talks about the conservation laws under global transformations, which has been extensively discussed in mechanics for homogeneous media [62]. The Nöether’s second theorem focuses on the conservation laws under local transformations, which is closely related to inhomogeneous media and will be emphasized in this paper.

This paper tries to answer the aforementioned three questions by using local gauge theories. For this purpose, Sect. 2 gives a brief review of the Nöether theorems and emphasizes the form-invariance of Lagrangian under local transformations. Then, different “Willis equations” are obtained in Sect. 3 by using the minimal replacement method [52, 57] regarding local temporal-spatial transformations. In this way, it demonstrates that the non-uniqueness of the Willis equations comes from the gauge freedom and they can be presented in unique forms that are easy to calculate by selecting appropriate gauge fixing schemes with clear physical meanings. Finally, the answers to the three questions posed in this section are given in Sect. 4.

2 A brief review of the Nöether theorems

As Fig. 1 shows, in linear elasticity we examine the deformation of a medium from the initial configuration B0 to the perturbed configuration B1. Correspondingly, a material point X initially located at the spatial position x0 in the global Cartesian coordinate system moves to the spatial position x1, resulting in a small displacement u = x1−x0. This deformation process should also satisfy some conservation laws, whose equivalence to the invariance or symmetry of the Hamiltonian action is revealed by the Nöether theorems [42, 46, 59,60,61,62].

Ignoring dissipation, the Hamiltonian action corresponding to the deformation from B0 to B1 is defined as:

$$\begin{array}{*{20}c} {A = \int_{\Omega } {L\left( {u_{i} ,u_{i,\alpha } } \right){\text{d}}\Omega } }\,\, & {\alpha = 1,2,3,t} & {i = 1,2,3} \\ \end{array} ,$$
(5)

where Ω is a temporal-spatial ___domain; Roman indices only represent spatial dimensions ranging from 1 to 3; and Greek indices can also represent temporal dimension denoted by t. In this paper, the spatial and temporal coordinates are not coupled with each other, i.e., ignoring the relativistic effect. The comma presented in the right subscript denotes the partial differentiation with respect to a temporal-spatial coordinate \(x_{\alpha }\).

In Eq. 5, L denotes the incremental Lagrangian from B0 to B1:

$$L\left( {u_{i} ,u_{i,\alpha } } \right) = W + \Phi - T,$$
(6)

where W, Φ and T are the incremental density of strain energy, potential of external body forces and kinetic energy, respectively. In this paper, we discuss only the first-order elasticity, so that L is the explicit function up to the gradient of displacement and the coordinate \(x_{\alpha }\) is implicitly included.

The stationary condition of Eq. 5 results in the Euler–Lagrange equation [61, 62]:

$$\frac{\partial L}{{\partial u_{i} }} - \left( {\frac{\partial L}{{\partial u_{i,\alpha } }}} \right)_{,\alpha } = \frac{\partial L}{{\partial u_{i} }} - \left( {\frac{\partial L}{{\partial u_{i,j} }}} \right)_{,j} - \left( {\frac{\partial L}{{\partial u_{i,t} }}} \right)_{,t} = 0,$$
(7)

where Einstein’s notation of index summation is adopted. Since B1 is unknown before solving Eq. 7, L should be represented in B0 with known quantities in linear elasticity. In this sense, each term in Eq. 7 has its clear physical meaning:

$$\overline{f}_{i} = - \frac{\partial L}{{\partial u_{i} }},\overline{s}_{ij} = \frac{\partial L}{{\partial u_{i,j} }},\overline{p}_{i} = - \frac{\partial L}{{\partial u_{i,t} }},$$
(8)

where \(\overline{f}_{i}\) is the effective body force; \(\overline{s}_{ij}\) is the total effective second Piola–Kirchhoff stress, which is usually not distinguished from the total Cauchy stress \(\overline{\sigma }_{ij}\) in the classical linear elasticity; and \(\overline{p}_{i}\) is the effective momentum. Therefore, Eq. 7 actually represents the equation of motion:

$$\overline{\sigma }_{ij,j} + \overline{f}_{i} = \overline{p}_{i,t} .$$
(9)

According to the Nöether theorems, a conservation law corresponds to a symmetry transformation of elasticity equations, or the invariance of the action A after infinitesimal continuous perturbations of both coordinates and field variables:

$$\begin{array}{*{20}l} {x_{\alpha } } & \to & {x^{\prime}_{\alpha } = x_{\alpha } + {\text{d}}x_{\alpha } } \\ \end{array} \equiv \varphi_{\alpha \beta } x_{\beta } ,$$
(10a)
$$\begin{array}{*{20}l} {u_{i} \left( {x_{\alpha } } \right)} & \to & {u^{\prime}_{i} \left( {x^{\prime}_{\alpha } } \right) = u_{i} + {\text{d}}u_{i} } \\ \end{array} \equiv \psi_{ij} u_{j} ,$$
(10b)

where \(\varphi_{\alpha \beta }\) and \(\psi_{ij}\) are the coordinate transformation operator and the gauge transformation operator, respectively. The coordinate transformation describes the change of viewpoint (or reference frame). The gauge transformation describes the corresponding change of field variables along with the coordinate transformation [43], or examines the effect of “transferring the magnitude of a vector to an infinitesimally close point” [44].

Since physical laws are independent of coordinates, the action A should be invariant under the transformation defined in Eq. (10):

$$\int_{ \Omega } {L\left( {u_{i} ,u_{i,\alpha } } \right){\text{d}}\Omega } = \int_{{\Omega^{\prime}}} {L\left( {u^{\prime}_{i} ,u^{\prime}_{{i,\alpha^{\prime}}} } \right){\text{d}}\Omega^{\prime}} .$$
(11)

In this paper, “\(,\alpha^{\prime}\)” in the right subscript denotes the partial differentiation with respect to a temporal-spatial coordinate \(x^{\prime}_{\alpha }\).

According to Eq. 10, we have:

$${\text{d}}{\Omega}^{\prime} = \det \left( {x^{\prime}_{\alpha ,\beta } } \right){\text{d}}{\Omega} \approx \left[ {1 + \left( {{\text{d}}x_{\alpha } } \right)_{,\alpha } } \right]{\text{d}}\Omega ,$$
(12a)
$$\begin{aligned} u^{\prime}_{{i,\alpha^{\prime}}} &= \left( {u_{i} + {\text{d}}u_{i} } \right)_{,\beta } x_{{\beta ,\alpha^{\prime}}} \\ & = \left( {u_{i} + {\text{d}}u_{i} } \right)_{,\beta } \left[ {\delta_{\beta \alpha } - \left( {{\text{d}}x_{\beta } } \right)_{,\alpha } } \right]\\ & \approx u_{i,\alpha } + \left( {{\text{d}}u_{i} } \right)_{,\alpha } - u_{i,\beta } \left( {{\text{d}}x_{\beta } } \right)_{,\alpha } ,\end{aligned}$$
(12b)

where \(\delta_{\alpha \beta}\) is the Kronecker delta symbol.

Substituting Eqs. 10 and 12 into Eq. 11, yields:

$$\begin{aligned} & \int_{\Omega } {L( {u_{i} ,u_{i,\alpha } } ){\text{d}}\Omega } \approx \\ & \int_{\Omega } {L\left( {u_{i} + {\text{d}}u_{i} ,u_{i,\alpha } + \left( {{\text{d}}u_{i} } \right)_{,\alpha } - u_{i,\beta } \left( {{\text{d}}x_{\beta } } \right)_{,\alpha } } \right)\left[ {1 + \left( {{\text{d}}x_{\alpha } } \right)_{,\alpha } } \right]{\text{d}}\Omega } \hfill \\ & \begin{array}{*{20}c} \Rightarrow & {\int_{\Omega } {\left\{ {\frac{\partial L}{{\partial u_{i} }}{\text{d}}u_{i} + \frac{\partial L}{{\partial u_{i,\alpha } }}\left[ {\left( {{\text{d}}u_{i} } \right)_{,\alpha } - u_{i,\beta } \left( {{\text{d}}x_{\beta } } \right)_{,\alpha } } \right] + L\left( {{\text{d}}x_{\alpha } } \right)_{,\alpha } } \right\}{\text{d}}\Omega } = 0} \\ \end{array} . \hfill \\ \end{aligned}$$
(13)

It should be noted that Eq. (13) is “accurate” in the sense of linear elasticity by ignoring higher order small terms. The following equations also follow this rule without further repeating.

Referring to the Euler–Lagrange equation defined in Eqs. 7, 13 can be rewritten as:

$$\int_{\Omega } {\left[ {\left( {\frac{\partial L}{{\partial u_{i,\alpha } }}{\text{d}}u_{i} } \right)_{,\alpha } - \frac{\partial L}{{\partial u_{i,\alpha } }}u_{i,\beta } \left( {{\text{d}}x_{\beta } } \right)_{,\alpha } + L\left( {{\text{d}}x_{\alpha } } \right)_{,\alpha } } \right]{\text{d}}\Omega } = 0.$$
(14)

Since

$$\begin{aligned} - \frac{{\partial L}}{{\partial u_{{i,\alpha }} }}u_{{i,\beta }} \left( {{\text{d}}x_{\beta } } \right)_{{,\alpha }} & = - \left( {\frac{{\partial L}}{{\partial u_{{i,\alpha }} }}u_{{i,\beta }} {\text{d}}x_{\beta } } \right)_{{,\alpha }} + \left( {\frac{{\partial L}}{{\partial u_{{i,\alpha }} }}} \right)_{{,\alpha }} u_{{i,\beta }} {\text{d}}x_{\beta } + \frac{{\partial L}}{{\partial u_{{i,\alpha }} }}u_{{i,\alpha \beta }} {\text{d}}x_{\beta }\\ & = - \left( {\frac{{\partial L}}{{\partial u_{{i,\alpha }} }}u_{{i,\beta }} {\text{d}}x_{\beta } } \right)_{{,\alpha }} + \left( {\frac{{\partial L}}{{\partial u_{i} }}u_{{i,\beta }} + \frac{{\partial L}}{{\partial u_{{i,\alpha }} }}u_{{i,\alpha \beta }} } \right){\text{d}}x_{\beta }\\ & = - \left( {\frac{{\partial L}}{{\partial u_{{i,\alpha }} }}u_{{i,\beta }} {\text{d}}x_{\beta } } \right)_{{,\alpha }} + L_{{,\alpha }} {\text{d}}x_{\alpha } , \\ \end{aligned}$$
(15a)
$$L\left( {{\text{d}}x_{\alpha } } \right)_{,\alpha } = \left( {L{\text{d}}x_{\alpha } } \right)_{,\alpha } - L_{,\alpha } {\text{d}}x_{\alpha } ,$$
(15b)

Equation 14 finally becomes:

$$\int_{\Omega } {P_{\alpha ,\alpha } {\text{d}}\Omega } = 0,$$
(16)

where \(P_{\alpha }\) is the Nöether current:

$$P_{\alpha } = \left( {L\delta_{\alpha \beta } - \frac{\partial L}{{\partial u_{i,\alpha } }}u_{i,\beta } } \right){\text{d}}x_{\beta } + \frac{\partial L}{{\partial u_{i,\alpha } }}{\text{d}}u_{i} .$$
(17)

Since Ω is arbitrary, the integrand in Eq. 16 should be:

$$P_{\alpha ,\alpha } = 0.$$
(18)

This defines certain conservation laws. Especially, \(L\delta_{\alpha \beta } - \partial L/\partial u_{i,\alpha } u_{i,\beta }\) corresponds to the famous Eshelby energy–momentum tensor (or the material momentum tensor) [61,62,63,64]; and \(\partial L/\partial u_{i,\alpha }\) corresponds to the stress tensor (or the physical momentum tensor) [62] as indicated in Eq. 8.

If \(\varphi_{\alpha \beta }\) and \(\psi_{ij}\) are independent of \(x_{\alpha }\), Eq. 10 defines global transformations, otherwise it defines local transformations. The conservation laws corresponding to global or local transformations are related to the Nöether’s first or second theorem, respectively [42, 46].

For homogeneous media, the corresponding Lagrangian \(L_{0}\) is usually easy to construct, so that researchers have tried to find out a finite number of global symmetry transformations that are related to classical conservation laws. For example, the conservation of energy is related to the global time translation; the conservation of momentum and angular momentum are related to the global physical translation and rotation, respectively; and the conservation of J-integral, L-integral and M-integral in fracture mechanics are related to the global material translation, rotation and scaling, respectively [62].

For inhomogeneous media, there could be an infinite number of local symmetry transformations. An important task is to construct the corresponding Lagrangian \(L_{1}\), which is invariant to local transformations. But the related conservation laws may have obscure physical interpretations. Although the local gauge theory has achieved great success in theoretical physics [41,42,43,44,45,46,47], its application to mechanics is just limited in geometric continuum mechanics [3, 48,49,50,51,52,53,54,55,56,57,58], which regards defects in the initial configuration B0 as transformed results from the virtual space Bv (see Fig. 1) according to their geometric characteristics such as the Burgers vector of dislocations. We attempt to slightly generalize this concept in Sect. 3.

3 Local gauge transformations

3.1 Minimal replacement method

Along with the coordinate transformation, if \(\psi_{ij}\) is a function of \(x_{\alpha }\), the Lagrangian of homogeneous media will change its form after gauge transformation, i.e., \(L_{0} \left( {u^{\prime}_{i} ,u^{\prime}_{{i,\alpha^{\prime}}} } \right) = L_{0} \left( {\psi_{ij} u_{j} ,\psi_{{ij,\alpha^{\prime}}} u_{j} + \psi_{ij} u_{{j,\alpha^{\prime}}} } \right)\), so that Eq. 11 is not valid any more. In order to construct a new Lagrangian \(L_{1}\) that is form-invariant to this local gauge transformation, we can use the minimal replacement method [52, 57] to define:

$$L_{1} = L_{0} \left( {u_{i} ,u_{i;\alpha } } \right).$$
(19)

The semicolon in Eq. 19 denotes the covariance derivative operator:

$$u_{i;\alpha } \equiv u_{i,\alpha } + \Gamma\,_{i\alpha }^{j} u_{j} ,$$
(20)

where the connection \(\Gamma\,_{i\alpha }^{j}\) describes the interaction among local freedoms in inhomogeneous media.

In Eq. 19, \(L_{0}\) is the Lagrangian of homogeneous media, which may contain a constant pre-stress \(\sigma_{ij}^{0}\). Since \(\sigma_{ij}^{0} = \sigma_{ji}^{0}\) and the elasticity tensor satisfies \(C_{ijkl} = C_{klij} = C_{ijlk} = C_{jikl}\), \(L_{0}\) can be explicitly written as:

$$L_{0} = \underbrace {{\sigma_{ij}^{0} u_{i,j} + \frac{1}{2}C_{ijkl} u_{i,j} u_{k,l} }}_{W}\underbrace {{ - \left( {\overline{f}_{i}^{0} + f_{i} } \right)u_{i} }}_{\Phi } - \underbrace {{\frac{1}{2}\rho_{ij} u_{i,t} u_{j,t} }}_{T},$$
(21)

where \(f_{i}\) is the incremental body force from B0 to B1; and \(\overline{f}_{i}^{0}\) is the effective body force in B0 that satisfies the equilibrium equation:

$$\sigma_{ij,j}^{0} + \overline{f}_{i}^{0} = 0.$$
(22)

Substituting Eq. 21 into Eq. 7 and referring to Eq. 22, we obtain a special form of Eq. 9, i.e., the well-known equation of motion for homogeneous media:

$$\sigma_{ij,j} + f_{i} = p_{i,t} ,$$
(23)

where \(\sigma_{ij} = \overline{\sigma }_{ij} - \sigma_{ij}^{0}\) is the incremental stress.

According to Eq. (19), the Lagrangian for inhomogeneous media can be constructed by replacing the differentiation operator in Eq. (21) with the covariance derivative operator defined in Eq. 20, i.e.:

$$\begin{aligned} L_{1} & = \sigma _{{ij}}^{0} u_{{i;j}} + \frac{1}{2}C_{{ijkl}} u_{{i;j}} u_{{k;l}} - \left( {\bar{f}_{i}^{0} + f_{i} } \right)u_{i} - \frac{1}{2}\rho _{{ij}} u_{{i;t}} u_{{j;t}}\\ & = \sigma _{{ij}}^{0} u_{{i,j}} + \sigma _{{ij}}^{0} {\Gamma\,} _{{ij}}^{r} u_{r} + \frac{1}{2}C_{{ijkl}} \left( {u_{{i,j}} u_{{k,l}} + 2u_{{i,j}} {\Gamma\,} _{{kl}}^{s} u_{s} +{\Gamma\,} _{{ij}}^{r} u_{r} {\Gamma\,} _{{kl}}^{s} u_{s} } \right) \\ & \quad- \left( {\bar{f}_{i}^{0} + f_{i} } \right)u_{i} - \frac{1}{2}\rho _{{ij}} \left( {u_{{i,t}} u_{{j,t}} + 2u_{{i,t}} {\Gamma\,} _{{jt}}^{s} u_{s} + {\Gamma\,} _{{it}}^{r} u_{r} {\Gamma\,} _{{jt}}^{s} u_{s} } \right), \\ \end{aligned}$$
(24)

where the conditions of \(C_{ijkl} = C_{klij}\) and \(\rho_{ij} = \rho_{ji}\) are used in the deduction.

It is easy to check that \(L_{1}\) is form-invariant with arbitrary \(\Gamma\,_{ij}^{r}\) and \(\Gamma\,_{it}^{r}\), so that we have much freedom to construct \(L_{1}\), which is certainly not unique. This is because we can freely choose a local coordinate system to describe the interaction among local freedoms. This also coincides with the observation that the Nöether current is always conservative if we add any divergence free vector in Eq. 17. For example, if \(P_{\alpha ,\alpha } = 0\), then \(\left( {P_{\alpha } + \varepsilon_{\alpha \beta \gamma } \,q_{\beta ,\gamma } } \right)_{,\alpha } = 0\) holds for any \(q_{\beta }\), where \(\varepsilon_{\alpha \beta \gamma }\) is the Levi–Civita symbol.

However, it could be dangerous if we have free rein to take any form of \(L_{1}\), which may lack physical meanings and difficult to calculate or measure. Fortunately, as what will be discussed in the following sub-sections, \(L_{1}\) can be constructed with slight freedoms by the understanding of the physical meanings of temporal-spatial transformations.

3.2 Physical meanings of temporal-spatial transformations

Mathematically, we can apply arbitrary coordinate and gauge transformations. However, it could be difficult to verify the resulted theory in experiments if the transformation lacks clear physical meanings [40]. In this section, we will discuss two types of local transformations that are easy to understand and can be verified in experiments.

Before proceeding to the detailed discussions of local temporal-spatial transformations, it should be noted that in order to describe the interaction among local freedoms in inhomogeneous media, we have to study the connection between adjacent material points in a gauge theory [43, 44]. This requirement of nonlocal theories for inhomogeneous media is different from local theories, where we trace the movement of a single material point. However in linear elasticity, all equations should be established in the reference configuration B0. Therefore in the following deductions, the local temporal transformation is discussed at a fixed spatial position related with different material points due to the infinitesimal change of time; and the local spatial transformation is discussed at a fixed time related with different material points at infinitely close spatial positions.

We firstly discuss the temporal transformation at a fixed spatial position \({{x}}^{0}\) in Fig. 1. Initially, two adjacent material points X and X′ located at spatial positions \({{x}}^{0}\) and \({{x^{\prime}}}^{0}\) in B0, respectively. With the advance of time from t to infinitely close tʹ, the material point X′ moves to \({{x}}^{0}\) in B1 (i.e., let \({{x^{\prime}}}^{1} = {{x}}^{0}\) in Fig. 1), so that the displacement at \({{x}}^{0}\) changes from \({{u}}\left( {{{x}}^{0} ,t} \right) = u_{i} {{e}}_{i}\) to \({{u^{\prime}}}\left( {{{x}}^{0} ,t^{\prime}} \right)\). The difference between \({{u^{\prime}}}\left( {{{x}}^{0} ,t^{\prime}} \right)\) and \({{u}}\left( {{{x}}^{0} ,t} \right)\) contains \(u_{i,t} {\text{d}}t{{e}}_{i}\) due to the change of time, and \(\partial u_{i} {/}\partial X_{j} {\text{d}}X_{j} {{e}}_{i}\) due to the change of material points from X to X′, where \({\text{d}}{{X}} = {{X^{\prime}}} - {{X}} = {{x^{\prime}}}^{0} - {{x}}^{0}\). Since the change of the spatial coordinate for X′ is \({\text{d}}{{x}} = {{x}}^{0} - {{x^{\prime}}}^{0} = - {\text{d}}{{X}}\), \(\partial u_{i} {/}\partial X_{j} = - \partial u_{i} {/}\partial x_{j}\). Then, we can write

$$\begin{aligned}{{u^{\prime}}}\left( {{{x}}^{0} ,t^{\prime}} \right) &= \left( {u_{i} + u_{i,t} {\text{d}}t - u_{i,j} {\text{d}}X_{j} } \right){{e}}_{i}\\& = \left( {u_{i} + u_{i,t} {\text{d}}t - u_{i,j} X_{j,t} {\text{d}}t} \right){{e}}_{i} .\end{aligned}$$

Notice that \({\text{d}}X_{j} = x_{j}^{{{\prime }0}} - x_{j}^{0} = u_{j}^{0}\), where \(u_{j}^{0}\) is the relative position change between these two material points from the stress-free state Bv to B0. Therefore, the material velocity \(X_{j,t}\) actually equals to the initial relative velocity \(u_{j,t}^{0}\) between adjacent material points. Thus, the incremental displacement at the spatial position \({{x}}^{0}\) due to the infinitesimal change of time is:

$${{u^{\prime}}}\left( {{{x}}^{0} ,t^{\prime}} \right) - {{u}}\left( {{{x}}^{0} ,t} \right) \approx \left( {u_{i,t} - u_{i,j} u_{j,t}^{0} } \right){\text{d}}t{{e}}_{i} .$$
(25)

In this case, it is easy to find that \(\left. {{\text{d}}u_{i} } \right|_{{{x}}} = \left( {u_{i,t} - u_{i,j} u_{j,t}^{0} } \right){\text{d}}t\) in Eq. 10b.

Then, we discuss the spatial transformation at a fixed time t. As Fig. 1 shows, the displacements of the two adjacent material points X and X′ can be represented as \({{u}}\left( {{{x}}^{0} ,t} \right) = u_{i} {{g}}_{i}\) and \({{u^{\prime}}}\left( {{{x^{\prime}}}^{0} ,t} \right) = u^{\prime}_{i} {{g^{\prime}}}_{i}\), respectively. \({{g}}_{i}\) and \({{g^{\prime}}}_{i}\) are basis vectors in material manifold or Lagrange coordinate system. Without loss of generality, \({{g}}_{i}\) is simply set as \({{e}}_{i}\), so that \({{g^{\prime}}}_{i} = x^{\prime}_{k,i} {{e}}_{k} \approx \left( {\delta_{ki} + u_{k,i}^{0} } \right){{e}}_{k}\). The initial displacement gradient \(u_{k,l}^{0}\) can represent the pre-strain \(e_{kl}^{0}\) that connects to the pre-stress \(\sigma_{ij}^{0}\) through the general Hooke’s law \(\sigma_{ij}^{0} = C_{ijkl} u_{k,l}^{0} = C_{ijkl} e_{kl}^{0}\). Therefore, the incremental displacement at a fixed time t due to the infinitesimal change of spatial ___location is:

$$\begin{aligned}{{u^{\prime}}}\left( {{{x^{\prime}}}^{0} ,t} \right) - {{u}}\left( {{{x}}^{0} ,t} \right) &\approx \left( {u_{i} + u_{i,j} {\text{d}}x_{j} } \right)\left( {\delta_{ki} + u_{k,ij}^{0} {\text{d}}x_{j} } \right){{e}}_{k} - u_{i} {{e}}_{i}\\& \approx \left( {u_{i,j} + u_{i,jk}^{0} u_{k} } \right){\text{d}}x_{j} {{e}}_{i} .\end{aligned}$$
(26)

In this case, it is easy to find that \(\left. {{\text{d}}u_{i} } \right|_{t} = \left( {u_{i,j} + u_{i,jk}^{0} u_{k} } \right){\text{d}}x_{j}\) in Eq. 10b.

For temporally homogeneous media, the initial relative displacement \(u_{j}^{0}\) does not change with time, so that \(u_{j,t}^{0} = 0\). For spatially homogeneous media, the initial displacement gradient \(u_{i,j}^{0}\) is a constant over space, so that \(u_{i,jk}^{0} = 0\). Hence, Eqs. 25 and 26 represent global gauge transformations, and the transformed Lagrangian \(L_{0}\) does not change its form.

For inhomogeneous media, \(u_{j}^{0}\) is time variant and \(u_{i,j}^{0}\) is position dependent, so that Eqs. 25 and 26 represent local gauge transformations. In these cases, \(L_{0}\) changes its form unless we can construct a form-invariant \(L_{1}\) according to Eq. 19. For these gauge transformations, the specific formulations of Eq. 20 can be written as:

$$u_{i;t} = u_{i,t} - u_{i,j} u_{j,t}^{0} ,$$
(27a)
$$u_{i;j} = u_{i,j} + u_{i,jk}^{0} u_{k} .$$
(27b)

In spatial transformations, \(u_{i,jk}^{0}\) is the connection \(\Gamma\,_{ij}^{k}\) (or gauge potential) in the material manifold and \(u_{i,j}^{0}\) serves as a gauge field. It is easy to check that \(u^{\prime}_{{i;j^{\prime}}} = \left( {\psi_{ik} u_{k} } \right)_{{;j^{\prime}}} = \psi_{ik} u_{{k;j^{\prime}}}\) when we adopt the transformation operator \(\psi_{ik} = x_{{k,i^{\prime}}}\) as that in [14]. Notice that \(u_{k,l}^{0}\) is closely related to the pre-stress \(\sigma_{ij}^{0}\) through the general Hooke’s law \(\sigma_{ij}^{0} = C_{ijkl} u_{k,l}^{0} = C_{ijkl} e_{kl}^{0}\). Actually, it is reasonable to take the pre-strain or the pre-stress to measure spatial inhomogeneity, because any spatial inhomogeneity (including defects, phase transformations, plastic deformations, constrained thermal deformations, etc.) is companied with inhomogeneous pre-strains or pre-stresses.

Although the explicit formulation of the connection \(\Gamma\,_{it}^{j}\) (and the corresponding \(\psi_{ij}\)) in temporal transformations is not clear, we can conveniently regard \(- u_{i,j} u_{j,t}^{0}\) as \(\Gamma\,_{it}^{j} u_{j}\), and the initial relative velocity \(u_{j,t}^{0}\) serves as a gauge field.

In a word, \(u_{j,t}^{0}\) and \(u_{i,j}^{0}\) describe the temporal and spatial interactions among local freedoms, so that the covariance derivative operators defined in Eq. 27 have clear physical meanings. In addition, it is possible to measure \(u_{j,t}^{0}\) indirectly [20, 21, 23], and pre-strains or pre-stresses can be routinely measured (e.g. by X-ray method).

The minimal replacement method has also been implemented to geometric continuum mechanics to construct the specific Lagrangian with respect to certain micro-defects, such as dislocations or disclinations [52] and damage [58], etc. In contrast to the connections constructed in this theory, Eq. 27 is just a small step forward not only using \(u_{j,t}^{0}\) to describe temporal transformations, but also using \(u_{i,j}^{0}\) to describe the general spatial transformations regardless of the detailed history from Bv to B0.

3.3 Local gauge transformation for the classical Willis equations

Ignoring spatial transformation, a simple case of Eq. 24 is just setting \({\Gamma\,_{ij}^{r}} = 0\). According to Eq. 27a, and neglecting the higher order small term \(\rho_{ij} \Gamma\,_{it}^{r} u_{r} \Gamma\,_{jt}^{s} u_{s}\), Eq. 24 becomes:

$$L_{1} = \sigma_{ij}^{0} u_{i,j} + \frac{1}{2}C_{ijkl} u_{i,j} u_{k,l} - \left( {\overline{f}_{i}^{0} + f_{i} } \right)u_{i} - \frac{1}{2}\rho_{ij} \left( {u_{i,t} u_{j,t} - 2u_{i,t} u_{s,t}^{0} u_{j,s} } \right).$$
(28)

Substituting Eq. 28 into Eq. 8 and neglecting the difference between the Cauchy stress and the second Piola–Kirchhoff stress, we obtain:

$$\overline{f}_{i} = - \frac{{\partial L_{1} }}{{\partial u_{i} }} = \overline{f}_{i}^{0} + f_{i} ,$$
(29)
$$\overline{p}_{i} = - \frac{{\partial L_{1} }}{{\partial u_{i,t} }} = \rho_{ij} u_{j,t} - \rho_{ij} u_{j,k} u_{k,t}^{0} ,$$
(30)
$$\sigma_{ij} = \overline{\sigma }_{ij} - \sigma_{ij}^{0} = \frac{{\partial L_{1} }}{{\partial u_{i,j} }} - \sigma_{ij}^{0} = C_{ijkl} u_{k,l} + \rho_{ik} u_{j,t}^{0} u_{k,t} .$$
(31)

According to Eqs. 9 and 22 and ignoring \(u_{k,tt}^{0}\), we get the linear equation of motion:

$$\sigma_{ij,j} + f_{i} = - \rho_{ij} u_{k,t}^{0} u_{j,kt} + \rho_{ij} u_{j,tt} .$$
(32)

Similar to the CWEs, the constitutive Eq. 31 contains a velocity coupling term and the equation of motion  32 has a term coupled with \(u_{j,kt}\). However, they are different in two aspects. One difference is that the coupling terms in the CWEs involve temporal-spatial convolution operators \(S_{{\left( {ij} \right)k}}^{{{\text{eff}}}}\) and \(\hat{S}_{{i\left( {jk} \right)}}^{{{\text{eff}}}}\) defined in Eqs.  50c and  50d in Appendix A, while Eqs. 31 and 32 just have algebraic coupling terms. Another difference is that the coupling terms in Eqs. 31 and 32 do not present such symmetries as \(S_{{\left( {ij} \right)k}}^{{{\text{eff}}}}\) and \(\hat{S}_{{i\left( {jk} \right)}}^{{{\text{eff}}}}\).

As Eq.  44 shows, the dynamic Green’s function is used to calculate displacements, so that the CWEs are actually in differential-integral formulations. Since a differential formulation could have its integral counterpart such as the boundary integration equations [65], it is not difficult to understand the first difference. This understanding can also give the answer to the third question posed in Sect. 1.

The second difference shows the non-uniqueness of the elastic equations for inhomogeneous media. This can be easily understood from the fact that Eqs. 48 and 49 in Appendix A just come from a single Eq. 44. This can also be regarded as the result of gauge freedoms mentioned at the end of Sect. 3.1. However, similar to electrodynamics [45], we can take the advantage of gauge freedoms to obtain proper equations that are easy to calculate by adopting an appropriate gauge fixing scheme. This is just like choosing an appropriate local coordinate system to facilitate calculating. For example, we can choose:

$$\begin{aligned}\rho_{ij} {\Gamma}\,_{jt}^{s} u_{s} &= - \frac{1}{4}\left( {\rho_{ij} u_{s,t}^{0} + \rho_{is} u_{j,t}^{0} } \right)\left( {u_{j,s} + u_{s,j} } \right) \\ & = - \frac{1}{2}\left( {\rho_{ij} u_{s,t}^{0} + \rho_{is} u_{j,t}^{0} } \right)e_{js} .\end{aligned}$$
(33)

Substituting Eq. 33 into Eq. 24, we obtain the following equations with the similar deduction:

$$\sigma_{ij} = C_{ijkl} u_{k,l} + \frac{1}{2}\left( {\rho_{jk} u_{i,t}^{0} + \rho_{ik} u_{j,t}^{0} } \right)u_{k,t} ,$$
(34)
$$\sigma_{ij,j} + f_{i} = - \frac{1}{2}\left( {\rho_{ij} u_{k,t}^{0} + \rho_{ik} u_{j,t}^{0} } \right)e_{kj,t} + \rho_{ij} u_{j,tt} .$$
(35)

These equations can be regarded as the differential counterparts of the CWEs, because their coupling terms present the same symmetries as \(S_{{\left( {ij} \right)k}}^{{{\text{eff}}}}\) and \(\hat{S}_{{i\left( {jk} \right)}}^{{{\text{eff}}}}\). This understanding implies that the CWEs could be regarded as the results from local temporal transformations defined in Eq. 25. Thus, it answers the second question posed in Sect. 1. It is also noticed that these equations will not change after adding an arbitrary constant \(\Delta_{i}\) to \(u_{i}^{0}\). This means they are form-invariant under arbitrary rigid translations, just like choosing a different origin of the local coordinate system, which is the essence of gauge freedoms.

Willis has tried to solve the non-uniqueness problem by introducing eigen-strains [36], which can be regarded as a special gauge fixing scheme. He also demonstrated that the equations in different forms can result in the same response [37]. Other “Willis equations” with complex-valued coupling terms [15,16,17,18] and acceleration coupling terms [19,20,21] may also be understood in the sense of gauge freedoms. However, it is fair to say that the CWEs or their differential forms of Eqs. 34 and 35 are elegant due to their symmetric forms with clear physical meanings and easy to calculate.

So far, we have shown that the CWEs can be obtained by either the local temporal gauge transformation illustrated in this section or the homogenization process briefly reviewed in Appendix A. However, in order to further clarify the underlying physical mechanism, we should go deeper into the dynamic behavior of microstructures.

As aforementioned, the tensorial density comes from the homogenization of the disordered movements of sub-components inside the microstructure. In the similar logic, we notice that the total kinetic energy density of all sub-components is \(T^{\prime} = \sum\limits_{M} {\rho^{M} u_{j,t}^{M} u_{j,t}^{M} /2}\) (M denotes the Mth sub-components), while the homogenized kinetic energy density is \(T = \rho_{ij} u_{i,t} u_{j,t} /2\). Usually \(T^{\prime} \ne T\) due to the presence of relative velocities among sub-components. In order to compensate the lost kinetic energy density \(T^{\prime} - T\) in the Lagrangian, we can convert it into the effective strain energy density \(\overline{W} = W - \left( {T^{\prime} - T} \right)\). Generally, we can use the initial relative velocity \(u_{j,t}^{0}\) to characterize the disorder pattern of a microstructure, and \(u_{i,j}\) to measure the magnitude of this disorder. Therefore, \(\overline{W}\) is the function of \(u_{i,t}\), \(u_{j,t}^{0}\) and \(u_{i,j}\), and consequently we have the velocity coupling term in Eqs. 31 and 34. In this sense, it is clear that the velocity coupling term does not represent damping. It is also easy to understand that \(u_{j,t}^{0}\) must be frequency dependent and vanish at zero frequency. Therefore, it is easy to obtain large coupling terms in unsymmetrical microstructures at local resonances [20,21,22,23]. However, \(\rho_{ij} \Gamma\,_{it}^{r} u_{r} \Gamma\,_{jt}^{s} u_{s}\) and \(u_{k,tt}^{0}\) could not be ignored near the resonance and the corresponding theory could be nonlinear [61]. In addition, even the microstructure is geometrically symmetric, we can still have nonzero coupling terms at certain frequencies with unsymmetrical modes. This understanding gives the answer to the first question posed in Sect. 1.

It is very interesting that since \(u_{j,t}^{0}\) serves as the inhomogeneous background velocity, \(\rho_{ik} u_{j,t}^{0} u_{k,t}\) in Eq. 31 and 34 could represent a generalized Coriolis effect. Therefore, the CWEs have intrinsic local directionality. This property has been utilized to design topological sonic crystals [66].

Referring to Eq. 25 we can substitute \(\left. {{\text{d}}u_{i} } \right|_{{{x}}} = \left( {u_{i,t} - u_{i,j} u_{j,t}^{0} } \right){\text{d}}t\) into Eqs. 17 and 18, and obtain:

\(\begin{gathered} \left[ {\left( {L_{1} \delta_{\alpha t} - \frac{{\partial L_{1} }}{{\partial u_{i,\alpha } }}u_{i,t} } \right){\text{d}}t + \frac{{\partial L_{1} }}{{\partial u_{i,\alpha } }}\left( {u_{i,t} - u_{i,j} u_{j,t}^{0} } \right){\text{d}}t} \right]_{,\alpha } = 0 \\ \qquad \qquad \qquad \left( {L_{1} {\text{d}}t + p_{i} u_{i,j} u_{j,t}^{0} {\text{d}}t} \right)_{,t} - \left( {\overline{\sigma }_{ik} u_{i,j} u_{j,t}^{0} {\text{d}}t} \right)_{,k} = 0. \\ \end{gathered}\)

Integrating the above relation over a volume V, and noticing \({\text{d}}t\) is arbitrary, we have the following conservation law:

$$\left[ {\int_{V} {\left( {L_{1} - p_{i} u_{i,j} u_{j,t}^{0} } \right){\text{d}}V} } \right]_{,t} = \int_{\partial V} {\overline{\sigma }_{ik} u_{i,j} u_{j,t}^{0} {\text{d}}S_{k} } ,$$
(36)

where \(\partial V\) denotes the boundary of V, and \(p_{i} u_{i,j}\) is the wave momentum [61, 62].

3.4 Local gauge transformation for the Willis-form equations

Ignoring temporal transformation, we can have another alternative of Eq. 24 by setting \({ \Gamma}\,_{jt}^{s} = 0\). According to Eq. 27b, we can set:

$${ \Gamma}\,_{ij}^{r} = u_{i,jr}^{0} .$$
(37)

Since \(\sigma_{ij}^{0} = C_{ijkl} u_{k,l}^{0}\), \(W^{0}\) defined in Appendix B can be rewritten as \(W^{0} = W\left( {x_{i}^{0} ,0} \right) = \sigma_{ij}^{0} u_{i,j}^{0} /2 = C_{ijkl} u_{i,j}^{0} u_{k,l}^{0} /2\). With this definition, \(\sigma_{ij}^{0} = \partial W^{0} /\partial u_{i,j}^{0}\) and \(C_{ijkl} = \partial^{2} W^{0} /\partial u_{i,j}^{0} \partial u_{k,l}^{0}\). Therefore, we have the following relations:

$$\sigma_{ij}^{0} { \Gamma}\,_{ij}^{r} = \frac{{\partial W^{0} }}{{\partial u_{i,j}^{0} }}\frac{{\partial u_{i,j}^{0} }}{{\partial x_{r} }} = W_{,r}^{0} ,$$
(38a)
$$C_{ijkl} { \Gamma}\,_{kl}^{s} = \frac{{\partial^{2} W^{0} }}{{\partial u_{i,j}^{0} \partial u_{k,l}^{0} }}\frac{{\partial u_{k,l}^{0} }}{{\partial x_{s} }} = \sigma_{ij,s}^{0} ,$$
(38b)
$$C_{ijkl} { \Gamma}\,_{ij}^{r} { \Gamma}\,_{kl}^{s} = \frac{{\partial^{2} W^{0} }}{{\partial u_{i,j}^{0} \partial u_{k,l}^{0} }}\frac{{\partial u_{i,j}^{0} }}{{\partial x_{r} }}\frac{{\partial u_{k,l}^{0} }}{{\partial x_{s} }} = W_{,rs}^{0} .$$
(38c)

Substituting Eqs. 37 and  38 into Eq. 24, we immediately get the same Lagrangian as that in Appendix B, so that the corresponding constitutive equation and the equation of motion should be the same as Eqs.  56 and 57 respectively.

In this case, the gauge field \(u_{i,j}^{0}\) is resulted from the interaction of material points via the gauge potential \({ \Gamma}\,_{ij}^{k} = u_{i,jk}^{0}\). It should be noted that \({ \Gamma}\,_{ij}^{k} \ne { \Gamma}\,_{ji}^{k}\). This implies that \({ \Gamma}\,_{ij}^{k}\) certainly contains torsion and can naturally represent the incompatible deformations induced by micro-defects such as dislocations and disclinations, etc. [3, 51,52,53,54,55]. In addition, \(\sigma_{ij,r}^{0}\) is invariant after \(u_{k,l}^{0}\) changes to \(u_{k,l}^{0} + \Delta_{kl}\) as long as \(\left( {C_{ijkl} \Delta_{kl} } \right)_{,r} = 0\). This implies that there are infinite choices of \(u_{k,l}^{0}\) without affecting the WFEs, presenting the gauge freedom under the condition of \(\left( {C_{ijkl} \Delta_{kl} } \right)_{,r} = 0\). As illustrated in Appendix B, the inhomogeneously distributed pre-stresses result in configurational forces, which are the sources of the coupling terms in the WFEs.

In addition, referring to Eq. 26 we can substitute \(\left. {{\text{d}}u_{i} } \right|_{t} = \left( {u_{i,j} + u_{i,jk}^{0} u_{k} } \right){\text{d}}x_{j}\) into Eqs. 17 and 18, and obtain:

$$\begin{aligned} \left[ {\left( {L_{1} \delta_{\alpha \beta } - \frac{{\partial L_{1} }}{{\partial u_{i,\alpha } }}u_{i,\beta } } \right){\text{d}}x_{\beta } + \frac{{\partial L_{1} }}{{\partial u_{i,\alpha } }}\left( {u_{i,j} + u_{i,jk}^{0} u_{k} } \right){\text{d}}x_{j} } \right]_{,\alpha }& = 0 \\ \left( {L_{1} {\text{d}}x_{k} + \overline{\sigma }_{ik} u_{i,jr}^{0} u_{r} {\text{d}}x_{j} } \right)_{,k} - \left( {p_{i} u_{i,jr}^{0} u_{r} {\text{d}}x_{j} } \right)_{,t} &= 0. \\ \end{aligned}$$

Integrating the above relation over a volume V, and noticing \({\text{d}}x_{j}\) is arbitrary, we have the following conservation law:

$$\left( {\int_{V} {p_{i} u_{i,jr}^{0} u_{r} {\text{d}}V} } \right)_{,t} = \int_{\partial V} {\left( {L_{1} \delta_{kj} + \overline{\sigma }_{ik} u_{i,jr}^{0} u_{r} } \right){\text{d}}S_{k} } .$$
(39)

4 Concluding remarks

In linear elasticity, since the deformed configuration B1 is unknown before solving the problem, all equations should be presented in the initial configuration B0. In addition, it is not necessary to either distinguish B1 from B0, or consider the non-local effect for homogeneous media, so that we have the classical linear elasticity theory in simple forms. However for inhomogeneous media, we have to consider the subtle effects due to the change of configurations and the nonlocal interactions among material points, so that we can find those coupling terms as discussed in Sect. 3 and Appendixes A and B.

With the discussion in Sect. 3, it is clear that the first-order inhomogeneous linear elasticity can be understood through the theory of local temporal-spatial transformations. Thus, the three questions posed in Sect. 1 can be easily answered:

  1. (1)

    The coupling terms related to velocity and strain rate in the CWEs have nothing to do with damping, but are the results of the material change at a fixed spatial position during the local temporal transformation or the compensation of the lost kinetic energy density during the homogenization process.

  2. (2)

    The CWEs are the results of local temporal transformations. The WFEs are the results of local spatial transformations.

  3. (3)

    The CWEs are differential-integral equations. Their pure differential counterparts are Eqs. 34 and 35.

With the help of gauge freedoms, it is also easy to understand why there are many different forms of “Willis equations” reported in literatures. Although all of them are mathematically valid, the CWEs could be the most elegant equations considering the temporal effect, and the WFEs could reveal the contribution of pre-strain gradients or pre-stress gradients. Of course, new equations will be obtained if temporal and spatial transformations are considered simultaneously.

The coupling terms in the CWEs and the WFEs are usually very weak, so that we can safely use the classical linear elasticity equations. However, these coupling terms should not be ignored at local resonances or locations with strong inhomogeneities, such as notches, crack tips or other interfaces.

Since the CWEs and the WFEs have already been partially validated by some experimental evidences [20, 21, 23, 32, 33], the classical principle of material frame indifference or objectivity could be extended to generalized versions that consider local transformations [40, 55, 57].

Weyl always used local gauges to describe the inhomogeneity due to the interaction of local freedoms in a field [46]. Although his original idea of applying gauge transformation to the scale of field variables failed in electromagnetic and gravitational fields [44], it seems that it has been revitalized in elastic field [68].