Online Material

Appendix A: Solution of the Duhamel-Neumann equation

We outline the main steps needed to solve the Duhamel-Neumann equation in Eq. (18) for our work. The temperature field, which produces the thermal stresses, is assumed to have a linearized form $T= T_{\rm av}+\Delta T$ where $\Delta T$ is given by Eq. (9). The uniqueness of the solution arises from (i) the regularity in the whole volume; and (ii) matching the free boundary conditions (given by Eq. (19)) at the surface r = R.

There are different ways in which we can decompose the displacement vector $\vec{u}$ into spherical-harmonics-type expansion (see, e.g., Thorne 1980, for an insightful review). Here we use a decomposition into spheroidal and toroidal components traditionally used in geophysical analyses (e.g., Kaula 1968; Bullen 1975). With this approach, related to what Thorne (1980) calls pure-spin vector harmonics, we have

$\begin{displaymath}\vec{u} = \sum\limits_{n=0}^\infty\sum\limits_{k=-n}^n\left(\... ...m S}_{nk}+\vec{u}^{\rm T}_{nk} \right)~\exp(\imath k\omega t), \end{displaymath}$

(A.1)

with

$\displaystyle \vec{u}^{\rm S}_{nk}$	=	$\displaystyle \left(\begin{array}{c} U_{nk}(r) \\ V_{nk}(r)~\frac{\partial }{\... ...heta}\frac{\partial }{\partial \phi} \end{array}\right) Y_{nk}(\theta,\phi)\; ,$	(A.2)
$\displaystyle \vec{u}^{\rm T}_{nk}$	=	$\displaystyle W_{nk}(r)\left(\begin{array}{c} 0 \\ -\frac{1}{\sin\theta}\frac{... ...\ \frac{\partial }{\partial \theta} \end{array}\right) Y_{nk}(\theta,\phi)\; .$	(A.3)

The first component in Eq. (A.1), $\vec{u}^{\rm S}_{nk}$ , is actually two separate spheroidal terms, that we call S1 and S2 in Sect. 3.1.2, characterized with radial-profile amplitudes U_nk(r) and V_nk(r). The second component in Eq. (A.1), $\vec{u}^{\rm T}_{nk}$ , is the toroidal component.

The spheroidal character of our source (temperature) term the Duhamel-Neumann equation implies two simplifications. First, the toroidal part of the displacement vector becomes negligible and we have W_nk=0. Second, we can restrict the summation over degrees n in Eq. (A.1) to the dipole and higher-order terms only, ignoring the monopole n=0. This is because the monopole part would correspond to purely radial temperature field, such as has been considered, for instance, in the previous works on our topic (e.g., Tambovtseva & Shestakova 1999; Shestakova & Tambovtseva 1997; Kuehrt 1984). Our temperature representation $T= T_{\rm av}+\Delta T$ does not contain a non-trivial, purely radial profile and the only viable free monopole term must have U₀₀=W₀₀=0to match the boundary conditions. Finally, we note that we also anticipated the Fourier-development structure in Eq. (A.1) as it follows from the source ( $\Delta T$ development).

Substituting the spheroidal-vector representation of $\vec{u}$ into the Duhamel-Neumann Eq. (18) we obtain the following system of equations for the radial profile of the amplitude functions U_nk(r) and V_nk(r) ( $n \geq 1$ ):

		$\displaystyle (\lambda+2~\mu)~r^2 \frac{{\rm d}^2}{{\rm d}r^2}U_{nk}(r) + 2(\lambda+2~\mu)~r \frac{{\rm d}}{{\rm d}r}U_{nk}(r)$
		$\displaystyle \qquad - n (n+1)(\lambda+\mu)~r \frac{\rm d}{{\rm d}r}V_{nk}(r)+ n(n+1)(\lambda+3~\mu)~V_{nk}(r)$
		$\displaystyle \qquad - \left[ ~2(\lambda+2\mu)+\mu~n(n+1)-k^2\omega^2r^2\rho \right]~ U_{nk}(r)$
		$\displaystyle \quad= \alpha~(3\lambda+2\mu)~r^2~\frac{{\rm d}}{{\rm d}r}T_{nk}(r),$	(A.4)

and

		$\displaystyle \mu~r^2\frac{{\rm d}^2}{{\rm d}r^2}V_{nk}(r) +2\mu~r\frac{\rm d}{{\rm d}r}V_{nk}(r) +(\lambda+\mu)~r\frac {\rm d}{{\rm d}r}U_{nk}(r)$
		$\displaystyle \qquad - \left[ n(n+1)(\lambda+2\mu)-k^2\omega^2r^2\rho \right]~V_{nk}(r)$
		$\displaystyle \qquad+ 2(\lambda+2\mu)~U_{nk}(r)=~\alpha~(3\lambda+2\mu)~r~T_{nk}(r).$	(A.5)

Here we found it useful to separate the r- and t-dependences of the t_nk(r,t) amplitudes of the $\Delta T$ development in Eq. (9) and introduce pure radial parts T_nk(r) such that $t_{nk}(r,t)=T_{nk}(r) \exp(\imath k\omega t)$ . These represent source terms in Eqs. (A.4) and (A.5).

While the solution of U_nk(r) and V_nk(r) is coupled by means of Eqs. (A.4) and (A.5), the fundamental implication of the Duhamel-Neumann equation linearity is that amplitude terms of different degrees and orders in the spherical harmonics development as well as the different Fourier modes are not mixed and can be solved separately.

Once we obtain U_nk(r) and V_nk(r), we can readily compute components of the corresponding stress tensor $\vec{\tau}^{\rm U}$ arising from the displacement vector field by using the Hook's law in Eq. (16). Given its linearity, we thus again have

$\begin{displaymath}\vec{\tau}^{\rm U} = \vec{\tau}(\vec{u}) = \sum\limits_{n=1}... ...ty\sum\limits_{k=-n}^n \vec{\tau}_{nk}~\exp(\imath k\omega t), \end{displaymath}$

(A.6)

where $\vec{\tau}_{nk}=\vec{\tau}(\vec{u}_{nk})$ . Projecting components of the stress tensor onto the orthonormal basis $(\vec{e}_{\rm r},\vec{e}_{\rm\theta}, \vec{e}_{\rm\phi})$ from Eq. (33), as outlined in Sect. 3.1, we obtain (for simplicity we dropped here the degree- and order-indexes n and k)

$\displaystyle \tau_{\rm rr}$	=	$$\displaystyle \frac{1}{r}\left[ \left( \lambda+2\mu \right)~r~$
			(A.7)
$\displaystyle \tau_{\rm r\theta}$	=	$\displaystyle \frac{\mu}{r}\left[ U_{nk}+r~\frac{{\rm d} V_{nk}}{{\rm d} r} -V_{nk} \right]\frac{\partial}{\partial\theta}Y_{nk},$	(A.8)
$\displaystyle \tau_{\rm r\phi}$	=	$\displaystyle \frac{\mu}{r}\left[ U_{nk}+r~\frac{{\rm d} V_{nk}}{{\rm d} r} -V_{nk} \right]\frac{\imath k}{\sin\theta} Y_{nk}\; ,$	(A.9)
$\displaystyle \tau_{\rm\theta\theta}$	=	$\displaystyle \frac{1}{r}\Bigl[2(\lambda+\mu)U_{nk}- \lambda n(n+1)V_{nk}$
	$\textstyle \quad +$	$\displaystyle \lambda~r~\frac{{\rm d} U_{nk}}{{\rm d} r}\Bigr]Y_{nk}+2\mu\frac{V_{nk}}{r} \frac{\partial}{\partial\theta}Y_{nk},$	(A.10)
$\displaystyle \tau_{\rm\phi \phi}$	=	$\displaystyle \frac{1}{r}\Bigl[2(\lambda+\mu)U_{nk} -n(n+1)(\lambda+2\mu)V_{nk}$
	$\textstyle \quad +$	$\displaystyle \lambda~r~\frac{{\rm d} U_{nk}}{{\rm d} r}\Bigr]Y_{nk} -2\mu~\frac{V_{nk}}{r}\frac{\partial^2 }{\partial \theta^2}Y_{nk},$	(A.11)
$\displaystyle \tau_{\rm\theta\phi}$	=	$\displaystyle 2\mu\frac{V_{nk}}{r}\frac{\imath k}{\sin\theta} \biggl(\frac{\partial}{\partial\theta}-\frac{\cos\theta}{\sin\theta} \biggr)Y_{nk}.$	(A.12)

The partial derivatives of the spherical functions Y_nk are computed using

$\displaystyle \frac{\partial }{\partial \theta}Y_{nk}$	=	$\displaystyle \frac{1}{2}\Bigl[\sqrt{(n-k)(n+k+1)}~ e^{-\imath\phi}~Y_{nk+1}$
	$\textstyle \quad -$	$\displaystyle \sqrt{(n+k)(n-k+1)}~e^{\imath\phi}~Y_{nk-1}\Bigr],$	(A.13)
$\displaystyle \frac{\partial^2 }{\partial \theta^2}Y_{nk}$	=	$\displaystyle \left[ \frac{k^2}{\sin^2\theta}-n(n+1) \right] Y_{nk}-\frac{\cos\theta}{\sin\theta}\frac{\partial }{\partial \theta}Y_{nk}.$	(A.14)

Equations (A.7) to (A.12) yield components of the stress tensor that explicitly depend on the displacement vector $\vec{u}$ . The part $\vec{\tau}^{\rm T}= -\alpha(3\lambda+2\mu)~ \Delta T~ {\bf 1}$ , which explicitly depends on the temperature (see the generalized Hook's law in Eq. (16), should be added separately to the total stress-tensor field. Because of the explicit analytical solution for $\Delta T$ , this is achieved at no computational expense.

Because of the linearity of the Duhamel-Neumann equation, a general solution is expressed in terms of a linear superposition of (i) a solution of the homogeneous system; and (ii) a particular solution of the inhomogeneous system. The next two sections discuss the two cases separately.

A.1 Solution of the homogeneous Duhamel-Neumann equation

Equations (A.4) and (A.5) with zero right-hand sides represent the homogeneous Duhamel-Neumann equation broken into parts corresponding to the individual spheroidal modes. Its solution is quite complicated, but may be significantly simplified in our case. This is because for the range of material parameters, sizes and rotation frequencies that apply for meteoroids we always have $\omega^2r^2\rho\ll \mu \sim \lambda+2\mu$ . With these we may neglect the troublesome term $k^2\omega^2r^2\rho$ in Eqs. (A.4) and (A.5). A major implication of this is then that the system of solutions of the homogeneous Duhamel-Neumann equation become degenerate in the k (order) index of the spherical-harmonics representation.

Adopting the aforementioned approximation, the homogeneous system of Eqs. (A.4) and (A.5) now has a form of Euler equations. As such, it has a fundamental system of power-law solutions Uⁱ_nk=Q_i r^m_i and Vⁱ_nk = r^m_i with $i=1, \ldots, 4$ , real-valued exponents m_i and amplitudes Q_i. After a straightforward algebra, we obtain

U¹_nk(r)	=	$\displaystyle (n+1)~\frac{n~\lambda+(n-2)~\mu}{(n+3)~\lambda+(n+5) ~\mu}~r^{n+1},$	(A.15)
V¹_nk(r)	=	rⁿ⁺¹,	(A.16)
U²_nk(r)	=	$\displaystyle n~r^{n-1},\quad V^2_{nk}(r) = r^{n-1},$	(A.17)
U³_nk(r)	=	$\displaystyle -(n+1)~r^{-n+2},\quad V^3_{nk}(r) = r^{-n+2},$	(A.18)
U⁴_nk(r)	=	$\displaystyle n~\frac{(n+1)~\lambda+(n+3)~\mu}{(n-2)~\lambda+ (n-4)~\mu}~r^{-n},$	(A.19)
V⁴_nk(r)	=	r^-n.	(A.20)

The last two modes, 3 and 4, diverge at the center r=0 and therefore must be excluded. We are thus left with the first two modes, 1 and 2, that produce the spheroidal modes $\vec{u}^{\rm S1}=\vec{u}^{\rm S} (U^1,V^1)$ and $\vec{u}^{\rm S2}=\vec{u}^{\rm S}(U^2,V^2)$ in Sect. 3.1.2 and whose associated stress field was given in Eqs. (44)-(49) and Eqs. (50)-(55). Obviously, they have also been used to obtain the stationary part of the stress field discussed in Sect. 3.1.1. We note that the second spheroidal mode represents a pure shear with no volumic changes (compression or expansion) because $\nabla\cdot \vec{u}^{\rm S2}=0$ .

A.2 Particular solution of the Duhamel-Neumann equation

We next find a particular solution of the inhomogeneous Duhamel-Neumann equation with the thermal source $T_{nk}(r)\neq 0$ . We divide this task into a discussion of the stationary case (k=0) and time-dependent case ( $k\neq 0$ ). In both cases, we again use the approximation of neglecting the $\propto (k\omega r \rho)^2$ terms in Eqs. (A.4) and (A.5).

A.2.1 Time-independent part

The stationary temperature field is given by $\Delta T=\sum_{n=1}^\infty~ C_{n0}(\theta_0)~r^n~Y_{n0}(\theta,\phi)$ (Eqs. (9) and (10)) and thus $T_{n0}(r)=C_{n0}(\theta_0)~r^n$ . We again search the fundamental system of solutions in a power-law form $U^{\rm P}_{n0}=Q_{\rm U}~r^{m_{\rm U}}$ and $V^{\rm P}_{n0}=Q_{\rm V} r^{m_{\rm V}}$ with some real-valued exponents $(m_{\rm U},m_{\rm V})$ and amplitudes $(Q_{\rm U},Q_{\rm V})$ . After a brief algebraic derivation, we obtain

$\displaystyle U_{n0}^{\rm P}$	=	$\displaystyle \frac{n+2}{2(2n+3)}\frac{\alpha(3\lambda+2\mu)}{ \lambda+2\mu}C_{n0}(\theta_0)~r^{n+1},$	(A.21)
$\displaystyle V_{n0}^{\rm P}$	=	$\displaystyle \frac{1}{2(2n+3)}\frac{\alpha(3\lambda+2\mu)}{ \lambda+2\mu}C_{n0}(\theta_0)~r^{n+1}.$	(A.22)

We note that this mode has the same radial profile as the $\vec{u}^{\rm S1}$ spheroidal model found above.

A.2.2 Time-dependent part

The time-dependent temperature field is given by $\Delta T= \sum_{n=1}^\infty\sum_{k=-n}^n~C_{nk}(\theta_0)~j_n(z_k)~Y_{nk} (\theta,\phi)~\exp(\imath k\omega t)$ with $k\ne 0$ (Eqs. (9) and (11)) and thus $T_{nk}=C_{nk}(\theta_0)~j_n(z_k)$ . We assume that the particular solution has a form $\vec{u}^{\rm P}=\nabla \Phi$ . Substituting this ansätz to the Duhamel-Neumann equation Eq. (18), we obtain

$\begin{displaymath}(\lambda+2\mu)(\nabla\cdot\nabla)\Phi-\rho\frac{\partial^2 }{\partial t^2}\Phi = \alpha(3\lambda+2\mu) \Delta T, \end{displaymath}$

(A.23)

where we have suitably assumed that the arbitrary constant on the right-hand side canceled the monopole (constant) temperature part. This is an inhomogeneous wave equation on a sphere that, however, takes a simple form because of the spherical-harmonic and Fourier structure of the source term $\Delta T$ on the right-hand side. Assuming thus a separation

$\begin{displaymath}\Phi = \sum_{n=1}^\infty\sum_{k=-n}^n \Gamma_{nk} \Sigma_{nk}(r) Y_{nk}(\theta,\phi)\exp(\imath k\omega t), \end{displaymath}$

(A.24)

we obtain

$\begin{displaymath}\Sigma_{nk}(r) = j_n(z_k), \end{displaymath}$

(A.25)

and

$\begin{displaymath}\Gamma_{nk} = \frac{\alpha}{\rho}\frac{3\lambda+2\mu}{k^2\omega^2 +\imath k(v_{\rm P}/\ell_{\rm d})^2}~ C_{nk}(\theta_0)\; , \end{displaymath}$

(A.26)

where $v_{\rm P}=\sqrt{(\lambda+2\mu)/\rho}$ is the elastic P-wave velocity as above. Adopting again the approximation $v_{\rm P}\gg \omega R \geq \omega\ell_{\rm d}$ , we may neglect the first term in the denominator of Eq. (A.26). Translating this solution into the amplitude-functions of the spheroidal-field representation (A.2), we finally obtain

$\displaystyle U^{\rm P}_{nk}$	=	$\displaystyle -\alpha\frac{3\lambda+2\mu}{\lambda+2\mu} C_{nk}(\theta_0)~\frac{\ell_{\rm d}}{\sqrt{-\imath k}} \frac{{\rm d}~ j_n(z_k)}{{\rm d}z},$	(A.27)
$\displaystyle V^{\rm P}_{nk}$	=	$\displaystyle -\alpha\frac{3\lambda+2\mu}{\lambda+2\mu} C_{nk}(\theta_0)~\frac{\ell_{\rm d}}{\sqrt{-\imath k}} \frac{j_n(z_k)}{z_k}.$	(A.28)

The corresponding stress tensor is expressed by Eqs. (58)-(63).

A.3 Complete expression of the thermal stress tensor

The complete solution of the Duhamel-Neumann equation is a linear combination of the free-spheroidal modes $\vec{u}^{\rm S1}$ and $\vec{u}^{\rm S2}$ from Sect. A.1 and the particular mode $\vec{u}^{\rm P}=\vec{u}(U^{\rm P},V^{\rm P})$ from Sect. A.2. In the individual spherical harmonics modes, we have $\vec{u}_{nk}=Q^1_{nk}\vec{u}^{\rm S1}_{nk}+Q^2_{nk} \vec{u}^{\rm S2}_{nk}+\vec{u}^{\rm P}_{nk}$ , where Q¹_nk and Q²_nk are some coefficients. We have to choose them to satisfy the surface boundary condition (19), namely $\vec{\tau}\cdot \vec{e}_{\rm r}=0$ at r=R. Here the total stress tensor is given by

$\begin{displaymath}\vec{\tau} = \vec{\tau}(\vec{u}) + \vec{\tau}^{\rm T}, \end{displaymath}$

(A.29)

or again in the spherical harmonics modes

$\begin{displaymath}\vec{\tau}_{nk} = \vec{\tau}(\vec{u}_{nk}) + \vec{\tau}^{\rm T}_{nk}. \end{displaymath}$

(A.30)

The truly active and independent conditions are $\tau_{\rm rr}=0$ and $\tau_{\rm r\theta}=0$ , from which the two constants Q¹_nk and Q²_nk follow. One easily checks that the third condition, $\tau_{\rm r\phi}=0$ , is always linearly dependent on $\tau_{\rm r\theta}=0$ (Eqs. (A.8) and (A.9)) and thus we do not need to consider it.

We were able to carry out all necessary algedraic manipulations and obtain a close form of the resulting formulae for the case of the stationary (zonal, k=0) part of the stress field. These are given in Eqs. (36)-(39) (Sect. 3.1.1). In the case of the time-dependent part of the stress field, the algebra is more involved and we could not reach as simple and compact results as for the time-independent part. We thus confine ourselves to provide formulae for the stress-tensor components of the individual components and those for the integration constants Q¹_nk and Q²_nk(Sect. 3.1.2).