© The Authors 2026

1 Introduction

Determination of the temperature averaged over the surface of bodies orbiting the Sun (thereafter, the mean temperature) is a useful information for studies of numerous physical processes in planetary science. We consider an analysis of the non-gravitational accelerations due to absorption and subsequent thermal reemission of sunlight, termed the Yarkovsky effect (review of its multiple applications in dynamics of small bodies in the Solar system may be found in Bottke et al. 2006; Vokrouhlický et al. 2015). The high-accuracy models of the Yarkovsky effect that can account for shape irregularities at all length scales, variable thermal parameters, or complex rotation state, require a fully numerical setup. This approach is needed, but generally requires a significant computational time and often does not provide information about the sensitivity of the results on multiple physical parameters. In contrast, analytical models of the Yarkovsky effect resort to approximations and thus only provide low-accuracy results, but are generally fast and offer the possibility to directly understand the dependence on the physical parameters. The aim of this paper is to generalize previous analytical models of the Yarkovsky effect.

The problem of the heat conduction and reemission is at the core of modelling the Yarkovsky effect. This requires solving the heat diffusion equation in the body with appropriate boundary conditions. In spite of approximating the shape by simple geometries (e.g., ellipsoidal or even spherical), the boundary condition that expresses the energy balance at the surface still presents a difficult problem because the thermally reemitted energy depends on the fourth power of the surface temperature and is therefore strongly non-linear. The analytical methods used to determine the Yarkovsky effect generally resort to a linearization of the thermal reemission term at the surface by assuming that the temperature can be split to a suitably chosen constant value and a small variation¹ (e.g., Afonso et al. 1995; Rubincam 1998; Farinella et al. 1998; Vokrouhlický 1998a,b; Mysen 2008). While providing reasonable results, the assumption of a constant mean temperature is fairly restrictive and adequate to the case of bodies on circular or low-eccentricity orbits.

The aim of this work is to provide a fully analytical formulation of the mean temperature that is suitable for either a linear or a non-linear analysis of the Yarkovsky effect that perturbs the motion of a body on an eccentric heliocentric orbit. While the principal part of the thermal effects modelling is intended for a forthcoming publication, the mean temperature should ideally contain information about part of the seasonal component (including the delay between absorption and reemission due to finite thermal inertia). In order to achieve this goal, we defined the mean temperature T_av as an exact solution of the heat diffusion problem and neglected diurnal effects. This allowed us to assume that T_av only has a radial profile inside the body and that the solar radiation flux is simply represented by its average value over the whole surface at a given heliocentric distance (Sect. 2). This setup resembles earlier thermal models of comets (see, e.g., Kuehrt 1984), with a few exceptions: no activity is assumed in our case, or the physical parameters are considered constant. On the other hand, cometary models are often restricted to modelling subsurface thermal processes assuming an isothermal core (the large-body limit in our terminology), while we treated all regimes of the ratio of the body size and penetration depth of the thermal effects. The cometary models mostly used a numerical approach, while we provide an analytical solution here. The method of its construction is discussed Sect. 3. Finally, we provide in Sect. 4 examples of our results and test the limits of their applicability. Because of the practical necessity of truncating the infinite series representing our solution, we are mostly concerned with the maximum value of the orbital eccentricity for which it can be used.

2 Formulation of the problem

We considered a spherical body with a radius R orbiting the Sun on an elliptical orbit with an eccentricity e. We assumed uniform rotation with the frequency ω about the spin axis s fixed in the inertial space. Adopting a body frame ℬ with the z-axis along the spin vector s, we wished to determine the temperature T(t; r) at an arbitrary time t and location r within the body. To do this, we needed to solve the heat conduction equation, $$\rho C\frac{\mathrm{\partial }T}{\mathrm{\partial }t}=K{\mathrm{\nabla }}^{2}T.$$(1)

For the sake of simplicity, we assumed all physical parameters, such as the bulk density ρ, the specific heat C, and the thermal conductivity K, to be constant. Given the assumed spherical shape of the body, we adopted spherical coordinates (r, θ, φ) to describe r in the body frame ℬ.

The linear dependence of Eq. (1) on T allowed us an analytical solution in terms of its development in a system of properly chosen eigenfunctions (Sect. 3). The free coefficients of this decomposition need to be determined by the boundary conditions. We assumed periodicity of the solution with a period P_rev of the revolution about the Sun for the time domain. Following Vokrouhlický (1999), we assumed that the rotation period P_rot of the body was commensurable with P_rev, such that P_rev/P_rot = m ≫ 1. Moreover, as shown in Farinella & Vokrouhlický (1996), it is possible to extend this result to the case of real m when the timescales of the variation in the spin-axis orientation are longer than the rotation and revolution periods. We assumed a constant spin rate and principal-axis rotation for simplicity. The boundary conditions of the space domain included (i) the regularity of T in the centre of the body r = 0, and (ii) the energy conservation at its surface, $$ϵ\sigma T^4+{K\frac{\mathrm{\partial }T}{\mathrm{\partial }r}|}_{r=R}=\alpha \mathcal{E},$$(2)

which expresses balance between the energy radiated by an arbitrary surface element using a grey-body approximation with emissivity ϵ (σ being the Stefan-Boltzmann constant), the energy conducted to the subsurface layers, and the solar radiation energy absorbed by the surface element. The latter was decomposed into the absorption coefficient α and the incident solar radiation flux projected onto the normal vector of the surface element.

The non-linear dependence of the boundary condition (2) on T represents a major obstacle to the analytic solution of the problem. This issue is traditionally overcome by assuming that T is close enough to a certain average value T_av, such that T(t; r = T_av + ΔT(t; r) and ΔT/T_av ≪ 1 for arbitrary t and r in the solution domain. Taking the approximation $T^4\approx T_{\mathrm{a}\mathrm{v}}^4+4T_{\mathrm{a}\mathrm{v}}^3\mathrm{\Delta }T$ implies linearity of the boundary condition (2) in ΔT, allowing us to handle the problem in analytical terms. The choice of T_av is an important part of the procedure. In this respect, most of the previous analyses adopted the simplest assumption of a constant value T_av. In our context, this is well justified for a body revolving about the Sun on a circular orbit, or on an orbit with a very low eccentricity e. However, when e becomes significant, it is difficult to choose a constant T_av value globally, such that the linearization of the T⁴ term in the boundary conditions (2) still holds.

To solve this problem, we propose to keep the linearization scheme, but we dropped the assumption of a constant T_av value. Instead, we considered T_av depending on time t and radial distance r from the body centre. This brings certain algebraic difficulty in expressing the solution analytically to already the average temperature part, but it leaves the principles how to solve ΔT as before. We postpone the solution of ΔT to a future work and focus here on the novel aspect of the approach, namely the determination of T_av.

T_av as a function of t and r must satisfy the following set of equations: $$\rho C\frac{\mathrm{\partial }{T}_{\mathrm{a}\mathrm{v}}}{\mathrm{\partial }t}=\frac{K}{r^2}\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(r^2\frac{\mathrm{\partial }}{\mathrm{\partial }r}\right)T_{\mathrm{a}\mathrm{v}},$$(3) $$ϵ\sigma T_{\mathrm{a}\mathrm{v}}^4+{K\frac{\mathrm{\partial }{T}_{\mathrm{a}\mathrm{v}}}{\mathrm{\partial }r}|}_{r=R}=\alpha \overline{\mathcal{E}},$$(4)

where $\overline{\mathcal{E}}=\overline{\mathcal{E}}(t)$ is a characteristic mean insolation of the body at time t. The previous analyses of the insolation of a sphere (e.g., Vokrouhlický 1998a; Čapek & Vokrouhlický 2010) have shown that $\overline{\mathcal{E}}=F(t)/4$, where F(t) is the solar radiation flux at time t. The time dependence in F is simply due to the changing heliocentric distance d along the elliptic orbit and F ∝ 1/d².

Before we proceed to the analytical solution of T_av, we simplified Eqs. (3) and (4) by appropriately scaling all variables. This step allowed us to highlight the single fundamental parameter of the problem (aside from the orbital eccentricity e). As in Vokrouhlický (1998a, 1999), we adopted the following replacement rules:

• t → ζ, with a complex variable ζ = e^ιℓ and ℓ the mean anomaly of the elliptic motion about the Sun and $\iota =\sqrt{-1}$;

• r → r′, defined by r′ = r/h_s, where $h_{\mathrm{s}}=\sqrt{K/\rho Cn}$ is the depth of the seasonal thermal wave and n the mean motion of the body about the Sun;

• F → F′, defined by F′ = F/F_a, where F_a is the solar radiation flux at a distance d equal to the semi-major axis a of the heliocentric orbit;

• $T_{\mathrm{a}\mathrm{v}}\to T_{\mathrm{a}\mathrm{v}}^{\mathrm{\prime }}=T_{\mathrm{a}\mathrm{v}}/T_a$, with the reference temperature T_a defined by $ϵ\sigma T_a^4=\alpha F_a$.

Plugging these new variables into Eqs. (3) and (4), we obtained (R′ = R/h_s), $$ı\frac{\mathrm{\partial }{T}_{\mathrm{a}\mathrm{v}}^{\mathrm{\prime }}}{\mathrm{\partial }\zeta }=\frac{1}{r^{\mathrm{\prime }2}}\frac{\mathrm{\partial }}{\mathrm{\partial }{r}^{\mathrm{\prime }}}\left(r^{\mathrm{\prime }2}\frac{\mathrm{\partial }}{\mathrm{\partial }{r}^{\mathrm{\prime }}}\right)T_{\mathrm{a}\mathrm{v}}^{\mathrm{\prime }},$$(5) $$T_{\mathrm{a}\mathrm{v}}^{\mathrm{\prime }4}+{\mathrm{\Theta }\frac{\mathrm{\partial }{T}_{\mathrm{a}\mathrm{v}}^{\mathrm{\prime }}}{\mathrm{\partial }{r}^{\mathrm{\prime }}}|}_{r^{\mathrm{\prime }}=R^{\mathrm{\prime }}}={\overline{\mathcal{E}}}^{\mathrm{\prime }},$$(6)

indicating that the problem depends on two fundamental parameters. The first, $\mathrm{\Theta }=\sqrt{K\rho C}\sqrt{n}/\left(ϵ\sigma T_a^3\right)$ is explicitly displayed in (6) and is usually called the thermal parameter. The second, the eccentricity e of the heliocentric orbit, is implicitly present in the expression of the scaled insolation term on the right side of (6), which reads $${\overline{\mathcal{E}}}^{\mathrm{\prime }}=\frac{1}{4}{\left(\frac{a}{d}\right)}^2.$$(7)

In order to express it using ζ and e, we used an elliptic expansion, $${\overline{\mathcal{E}}}^{\mathrm{\prime }}=\frac{1}{4}\sum _{n\in Z}{a}_n(e){\zeta }^n,$$(8)

where a_n are the Hansen coefficients $a_n=X_n^{-2,0}(e)$, obeying an important property a_n ≃ 𝒪(e^|n|). The zero-order coefficient is $a_0=1/\sqrt{1-e^2}$. Because of the faster convergence, we chose to expand the a_n coefficients in terms of $\beta =\frac{e}{1+\sqrt{1-e^2}}$ instead of e. The leading coefficient thus reads a₀ = (1 + β²)/(1 − β²), while in general, we have $$a_n={\beta }^{|n|}\sum _{\ell =0}^{\mathrm{\infty }}{a}_{n,|n|+2\ell }{\beta }^{2\ell },$$(9)

with a certain set of numerical constants a_i,j (e.g., Tisserand 1894; Breiter et al. 2004). The lowest-degree coefficients a_n (n ≤ 8) read $$\begin{array}{rl}& a_0=1+2{\beta }^2+2{\beta }^4+2{\beta }^6+2{\beta }^8,\\ & a_1=2\beta +{\beta }^3+\frac{23}{6}{\beta }^5+\frac{73}{72}{\beta }^7,\\ & a_2=5{\beta }^2-\frac{22}{3}{\beta }^4+\frac{76}{3}{\beta }^6-\frac{2254}{45}{\beta }^8,\\ & a_3=13{\beta }^3-\frac{181}{4}{\beta }^5+\frac{1267}{8}{\beta }^8,\\ & a_4=\frac{103}{3}{\beta }^4-\frac{2834}{15}{\beta }^6+\frac{36074}{45}{\beta }^8,\\ & a_5=\frac{1097}{12}{\beta }^5-\frac{49531}{72}{\beta }^7,\\ & a_6=\frac{1223}{5}{\beta }^6-\frac{81722}{35}{\beta }^8,\\ & a_7=\frac{47273}{72}{\beta }^7,\\ & a_8=\frac{556403}{310}{\beta }^8,\end{array}$$

where terms to the order eight in β were included.

We point out that the above expansion are convergent only when e < e_L, where e_L ~ 0.66 is the Laplace limit (Da Silva Fernandes 1995). The iterative procedure described in the following section relies on the convergence of the Hansen coefficients expansion, and we therefore consider an object with eccentricity e < e_L in the following.

3 The average temperature solution

In this section, we provide the algorithm for determining $T_{\mathrm{a}\mathrm{v}}^{\mathrm{\prime }}$. The linearity of the heat conduction equation (5) makes the solution still simple, but the eigenfunctions of the differential operator on the right side should be observed. These are the spherical Bessel functions of the first kind and degree zero j₀ (the Bessel functions of the second kind are excluded by regularity of the solution at r′ = 0), while the eigenfunctions of the time operator on the right side of (5) are simply powers ζⁿ (or Fourier series in ℓ if set in real notation). The general solution thus reads $$T_{\mathrm{a}\mathrm{v}}^{\mathrm{\prime }}(\zeta ,r^{\mathrm{\prime }})=\sum _{n\in \mathbb{Z}}{c}_n{\zeta }^n{j}_0\left(z_n\right),$$(10)

where $z_n=\sqrt{-in}{r}^{\mathrm{\prime }}$. We note that $T_{\mathrm{a}\mathrm{v}}^{\mathrm{\prime }}$ must be real, which implies that c_n are complex coefficients and c_−n must be complex-conjugated to c_n. Their value has to be determined from the boundary condition (6), with the algebraic challenge to deal with non-linearity of the $T_{\mathrm{a}\mathrm{v}}^{\mathrm{\prime }4}$ term. In particular, plugging in the solution (10) directly to (6), and separating the n = 0 terms, we have $$\begin{array}{rl}& \stackrel{T_{\mathrm{a}\mathrm{v}}^{\mathrm{\prime }4}}{\overbrace{{[c_0+\sum _{n\ne 0}{c}_n{\zeta }^n{j}_0\left(Z_n\right)]}^4}}+\stackrel{{\mathrm{\Theta }\frac{\mathrm{\partial }}{\mathrm{\partial }{r}^{\mathrm{\prime }}}{T}_{\mathrm{a}\mathrm{v}}^{\mathrm{\prime }}|}_{r=R^{\mathrm{\prime }}}}{\overbrace{\mathrm{\Theta }\frac{\mathrm{\partial }}{\mathrm{\partial }{r}^{\mathrm{\prime }}}{\left[\sum _{n\ne 0}{c}_n{\zeta }^n{j}_0\left(z_n\right)\right]}_{z_n=\sqrt{-in}{R}^{\mathrm{\prime }}}}}\\ & =\frac{a_0}{4}+\frac{1}{4}\sum _{n\ne 0}{a}_n{\zeta }^n.\end{array}$$(11)

Further algebra simplifies with replacing c_n by C_n = c_n j₀ (Z_n), where $z_n=\sqrt{-in}{R}^{\mathrm{\prime }}$, and by introducing $${\psi }_n={\frac{z}{j_0(z)}\frac{d{j}_0}{dz}|}_{z=Z_n}.$$(12)

The boundary condition (11) then reads $${[C_0+\sum _{n\ne 0}{C}_n{\zeta }^n]}^4+\frac{\mathrm{\Theta }}{R^{\mathrm{\prime }}}\sum _{n\ne 0}{C}_n{\psi }_n{\zeta }^n=\frac{a_0}{4}+\frac{1}{4}\sum _{n\ne 0}{a}_n{\zeta }^n.$$(13)

The remaining algebra, at higher orders performed with the aid of computer software, stems from (i) comparing terms of the same power in ζ on either side of Eq. (13), and (ii) at the same time developing C_n in power series of β, namely² $$C_n={\beta }^{|n|}\sum _{\ell =0}^{\mathrm{\infty }}{C}_{n,|n|+2\ell }{\beta }^{2\ell }$$(14)

similar to the series expansion of a_n in (9). Again, comparison of terms of the same power in β yields a general scheme to determine the numerical coefficients C_i,j. We first illustrate the procedure in Sect. 3.1 by pushing the algorithm to the lowest order in both ζ and β, and we then present the general algorithm in Sect. 3.2.

3.1 First-degree approximation

We describe details of the procedure for computing the average temperature coefficients in the case of a first-degree approximation, namely, approximating $$T_{\mathrm{a}\mathrm{v}}^{\mathrm{\prime }}=C_0+C_1\zeta +C_{-1}{\zeta }^{-1}.$$(15)

Expanding the fourth-power in the boundary condition, we obtain $$\begin{array}{rl}{T}_{\mathrm{a}\mathrm{v}}^{\mathrm{\prime }4}=& C_0^4+12C_0^2{C}_1{C}_{-1}+6{\left(C_1{C}_{-1}\right)}^2\\ & +4C_0(C_0^2+3C_1{C}_{-1})(C_1\zeta +C_{-1}{\zeta }^{-1})+\dots ,\end{array}$$(16)

and the boundary condition (13) reads $$\begin{array}{rl}{C}_0^4& +6C_1{C}_{-1}[2C_0^2+C_1{C}_{-1}]+4C_0(C_0^2+3C_1{C}_{-1})(C_1\zeta +C_{-1}{\zeta }^{-1})\\ & +\frac{\mathrm{\Theta }}{R^{\mathrm{\prime }}}(C_1{\psi }_1\zeta +C_{-1}{\psi }_{-1}{\zeta }^{-1})+\dots =\frac{1}{4}\frac{1-{\beta }^2}{1+{\beta }^2}+\frac{1}{4}(a_1\zeta +a_{-1}{\zeta }^{-1})+\dots .\end{array}$$(17)

Recalling that a_n are power series of β such that a_n ~ β^|n|, we compare terms of the same powers in (ζ, β) on either side of Eq. (17). A brief algebra yields to the order β³ $$\begin{array}{rl}& C_{00}^4+4C_{00}^3{C}_{02}{\beta }^2+12C_{00}^2{C}_{11}^2{C}_{-11}^2{\beta }^2+4C_{00}^3{C}_{11}\beta \zeta +\frac{\mathrm{\Theta }}{R^{\mathrm{\prime }}}{C}_{11}{\psi }_1\zeta \\ & +C.C.=\frac{1}{4}+\frac{1}{2}{\beta }^2+\frac{1}{2}\beta \zeta +C.C.,\end{array}$$(18)

where C.C. stands for the complex conjugate terms. Observing equal powers in β and ζ on either side of (18), we obtain $$C_{00}=\frac{1}{\sqrt{2}}$$(19)

at the zero order (completing C₀ to 𝒪(β²)), and $$C_{11}=\frac{1}{2\sqrt{2}}\frac{1}{1+\frac{\mathrm{\Theta }}{\sqrt{2}{R}^{\mathrm{\prime }}}{\psi }_1}$$(20)

at the first order (completing C₁ to 𝒪(β³)). Clearly, the C₋₁₁ coefficient is obtained by taking the complex conjugate of C₁₁.

Remark. Although in principle we could fix the degree of the approximation and increase the order, the formulas behave better when we keep the order of the approximation lower than the maximum degree, as shown in the following.

3.2 Algorithm

In this section, we describe an iterative algorithm for computing the coefficients C_nj, which are functions of the thermal parameter Θ and of the normalized radius R′. The following properties imply that the algorithm is well-defined.

Proposition.

1. Linearity, that is, the coefficient C_nj, can be computed through a first-order linear relation that only involves terms of lower degree and order.

2. C_nj only depends on the terms C_mk for which the following relations hold: $$0\le k+m\le j+n,$$(21) $$0\le k-m\le j-n.$$(22)

Proof.

1. We focus on a generic coefficient C_nj. First of all, we note that by definition, we have that $$j\ge |n|,$$(23)

   so that j + n ≥ 0 and j − n ≥ 0. Moreover, we have that j ≡ n (mod 2). The lowest degree and order at which the coefficient C_nj appears in the boundary condition (18) is degree n and order j. We focus on the left side of the boundary condition equation. Because of the definition of $T_{\mathrm{a}\mathrm{v}}^{\mathrm{\prime }}$ and of C_nj, the above statement is straightforward for the second sum. On the other hand, we consider the form of a generic term in the multinomial expansion of the fourth-power term. This term is proportional to $$C_{i_1}^{k_1}{C}_{i_2}^{k_2}{C}_{i_3}^{k_3}{C}_{i_4}^{k_4}{\zeta }^{k_1{i}_1+k_2{i}_2+k_3{i}_3+k_4{i}_4},$$(24)

   where the positive integers k₁, k₂, k₃, and k₄ verify $$k_1+k_2+k_3+k_4=4.$$(25)

   Therefore, we can focus on terms of the form $$\begin{array}{l}{C}_p^4{\zeta }^{4p}{C}_p^3{C}_q{\zeta }^{3p+q}{C}_p^2{C}_q^2{\zeta }^{2(p+q)}\\ C_p^2{C}_q{C}_r{\zeta }^{2p+q+r}{C}_p{C}_q{C}_r{C}_s{\zeta }^{p+q+r+s}.\end{array}$$(26)

   These relations can be readapted in terms of the expansion in series of β, $$\begin{array}{rl}& C_{pg}^4{\zeta }^{4p}{\beta }^{4g}{C}_{pg}^3{C}_{qh}{\zeta }^{3p+q}{\beta }^{3g+h}{C}_{pg}^2{C}_{q,h}^2{\zeta }^{2(p+q)}{\beta }^{2(g+h)}\\ & C_{pg}^2{C}_{qh}{C}_{ri}{\zeta }^{2p+q+r}{\beta }^{2g+h+i}{C}_{pg}{C}_{qh}{C}_{ri}{C}_{sl}{\zeta }^{p+q+r+s}{\beta }^{g+h+i+l}.\end{array}$$(27)

   Therefore, the lowest degree and order term in which C_nj appears must have the form $$C_{00}^3{C}_{nj}{\zeta }^n{\beta }^j.$$(28)

   Therefore, because of the previous considerations, it is possible to obtain C_nj as the solution of an equation of the form $$F{C}_{nj}=G,$$(29)

   where F and G are functions of the coefficients C_mk with m < n and k < j.

2. We present a maximal set of coefficient that appear in the definition of C_nj. To do this, we focus on the terms of degree n and order j in the boundary condition. Selecting these terms is trivial in the case of the right side of the equation, while it is more challenging in the case of the left side because of the fourth-power term. By focusing on the form of the generic terms appearing in the expansion (which we presented above), we develop formulas that identify the required coefficients. Without loss of generality, we focus on a term of the form $$C_{00}^2{C}_{mk}{C}_{n-m,j-k},$$(30)

   which is a term of degree n and order j. The relation involving the indices of the required coefficients is obtained by remarking that because of condition (23), we have that it must hold that k ≥ |m| and $$j-k\ge |n-m|.$$(31)

   Therefore, we have that $$k\pm m\ge 0$$(32)

   and that $$j-k\pm |n-m|\ge 0.$$(33)

   These two relations imply the thesis. Focusing on terms C_mk with m, k ≠ 0 that appear with an exponent higher than 1 would naturally yield only a subset of the coefficients described by the relation above.

3.2.1 Description of the algorithm

We started by selecting the maximum degree n_max at which we truncated the series for the normalized average temperature, namely $$T_{\mathrm{a}\mathrm{v}}^{\mathrm{\prime }}\simeq \sum _{n=-n_{max}}^{n_{max}}{C}_n{\zeta }^n.$$(34)

Because of the properties of the second part of the above proposition, an efficient method of computing the coefficient is the following.

1. First, compute all the diagonal coefficients C_nn up to n = n_max and with n ≥ 0. Because of the properties of the C_nj, by taking the complex conjugate of the diagonal coefficients, we obtain C_−n,n n = 1, . . ., n_max.

2. Compute the coefficients in the above diagonal as shown by the blue arrows in Fig. 1, i.e. the C_n,n+2 up to n = n_max −2.

3. Iterate the procedure by focusing on the next diagonal until all the coefficients are computed.

In principle, the procedure described above could be extended to arbitrary degree and order, but it is well defined only within the radius of convergence of the series expansion of the Hansen coefficients, that is, e < e_L.

Remark. A key feature of this procedure is that when we wish to increase the fidelity of the model by increasing the degree of the truncation, it is not necessary to recompute all coefficients, but only the new ones, as shown by the shaded areas in Fig. 1.

3.2.2 Example: computation of the coefficient up to degree and order 2

After outlining the general algorithm, we found it useful to illustrate it in the case of pushing the first-order solution from Sect. 3.1 to the second order, determining thus coefficients C₀₂ and C₂₂. To compute the diagonal coefficient C₂₂, we focused on the boundary condition. Given that $C_{00}=1/\sqrt{2}$ and C₁₁ in (20), the terms of degree and order 2 of the boundary condition read $$3C_{11}^2+\sqrt{2}{C}_{22}+\frac{\mathrm{\Theta }}{R^{\mathrm{\prime }}}{C}_{22}=\frac{5}{4},$$(35)

from which we obtained $$C_{22}=\frac{5-12C_{11}^2}{4\sqrt{2}}\frac{1}{1+\frac{\mathrm{\Theta }}{\sqrt{2}{R}^{\mathrm{\prime }}}{\psi }_2}.$$(36)

To compute C₀₂, we proceeded similarly, but focused this time on the terms of degree 0 and order 2 of the boundary conditions, which read $$\sqrt{2}{C}_{02}+6C_{11}{C}_{-11}=\frac{1}{2}.$$(37)

Solving for C₀₂, we obtained $$C_{02}=\frac{1-12C_{11}{C}_{-11}}{2\sqrt{2}}.$$(38)

In the next section, we describe the asymptotic behaviour of the coefficients with respect to the normalized radius R′ in order to simplify formulas similar to Eq. (35).

Fig. 1 Flow chart for determining coefficients Cnk up to nmax = 4. The algorithm proceeds in the order indicated by the blue arrows. The shaded areas show that when nmax is increased to 5, most of the needed coefficients are already known, and it suffices to only compute three new ones: C15, C35, and C55.

3.3 Limits of large and small bodies

The surface temperature, determined by coefficients C_n, is a non-trivial function of the scaled size R′, principally due to their dependence on complex factors ψ_n. However, as a general rule of thumb, there are two regimes with a transition close to R′ ≃ 1. This is when the physical size is comparable to the penetration depth of the seasonal thermal wave³. For lower R′ values, the imaginary part of C_n decreases inversely proportionally to R′, while the real part of C_n converges to particular values independent of Θ. In contrast, for higher R′ values, the real and imaginary parts of C_n both converge to constant values dependent on Θ. The asymptotic values in both limits also depend on the orbital eccentricity. We now describe the two limiting regimes of small (R′ ≪ 1) and large (R′ ≫ 1) bodies with a few more details. To do this, we introduce the following auxiliary functions (n ≥ 0, the case of negative n values is not needed since ψ_−n is complex conjugate to ψ_n): $$G_n^{\pm }\equiv 1-e^{-4x_n}\pm 2e^{-2x_n}\sin 2x_n,$$(39) $$F_n\equiv 1+e^{-4x_n}-2e^{-2x_n}\cos 2x_n,$$(40)

which are related to the real and imaginary parts of ψ_n as follows: $${\psi }_n=\frac{x_n(G_n^{+}+ıG_n^{-})}{F_n}-1,$$(41)

with $x_n=\sqrt{n/2}{R}^{\mathrm{\prime }}$. The imaginary part of ψ_n, namely $x_n{G}_n^{-}/F_n$, plays the key role in the time offset of the thermal response to orbital forcing represented by elliptic expansion in ζⁿ (see Eq. (13)). In physical terms, it expresses a delay between absorption of sunlight and thermal reemission by the body.

Small-body limit. In the limit of small bodies, x_n → 0, we have $$\frac{x_n{G}_n^{+}}{F_n}-1\approx -\frac{4}{15}{x}_n^4+\mathcal{O}\left(x_n^5\right)$$(42) $$\frac{x_n{G}_n^{-}}{F_n}\approx -\frac{2}{3}{x}_n^2+\mathcal{O}(x_n^3),$$(43)

and thus, $${\frac{{\psi }_n}{x_n}|}_{x_n\to 0}\to 0.$$(44)

This implies that the heat conduction term in the boundary condition (13) is proportional to Θ R′. For a bound value of the thermal inertia, thus Θ, such that Θ R′ → 0, the boundary condition takes a simpler form, $${(C_0+\sum _{n\ne 0}{C}_n{\zeta }^n)}^4=\frac{1}{4}\frac{1-{\beta }^2}{1+{\beta }^2}+\frac{1}{4}\sum _{n\ne 0}{a}_n{\zeta }^n,$$(45)

from which the coefficients C_n can again be determined using the above outlined algorithm. Obviously, (45) also implies $$T_{\mathrm{a}\mathrm{v}}^{\mathrm{\prime }4}={\overline{\mathcal{E}}}^{\mathrm{\prime }},$$(46)

which may be used to verify the result. In physical terms, this limit simply means that very small bodies become isothermal (due to efficient heat conduction across their volume) and the absorbed sunlight energy is instantly re-radiated (since these small bodies cannot store the energy and radiate it with a certain delay). Remarkably, the temperature is just the equilibrium temperature corresponding to the heliocentric distance.

Large-body limit. In the limit of large bodies, x_n → ∞, we have G^± ≈ F ≈ 1, and therefore, $$\frac{{\psi }_n}{x_n}\approx 1+ı-\frac{1}{x_n}\underset{x_n\to \mathrm{\infty }}{\overset{}{\to }}1+ı,$$(47)

and the boundary condition (13) takes a simplified form, $${(C_0+\sum _{n\ne 0}{C}_n{\zeta }^n)}^4+\frac{\mathrm{\Theta }}{\sqrt{2}}\sum _{n\ne 0}{C}_n\sqrt{|n|}(1\pm i){\zeta }^n=\frac{1}{4}\frac{1-{\beta }^2}{1+{\beta }^2}+\frac{1}{4}\sum _{n\ne 0}{a}_n{\zeta }^n,$$(48)

with the positive and negative sign in the bracket of the second term forming positive and negative values of n. Explicit formulas for the C_nj coefficients, computed up to degree and order eight in the large-body limit are presented in Appendix B. Because the penetration depth h_s becomes very small compared to the radius R in this limit, all temperature variations become limited to a very narrow skin below the body surface. A description of this situation would benefit from changing the spherical coordinate r, measured from the centre, to a subsurface depth z = R − r (e.g., Vokrouhlický 1998b).

4 Surface mean temperature in the large-body limit

We illustrate the mean-temperature algorithm outlined above for a large body with a fixed elliptic orbit. For the sake of the example, we used the orbit and physical parameters of asteroid (65803) Didymos (see Table 1) because of the growing interest in this object as a target of the ESA Hera mission Michel et al. (2022) and Tortora et al. (2025). Our analytical approach is limited, and we therefore only replaced the true Didymos shape by a sphere of an equivalent volume. Assuming thermal conductivity in a plausible range 0.01–0.1 W m⁻¹ K⁻¹ (e.g., Delbó et al. 2015), we estimated that the penetration depth h_s of the seasonal thermal wave is between 0.1 and 1 m, which justifies the large-body model for this asteroid well. The data in Table 1 also allowed us to estimate the seasonal thermal parameter Θ ≃ 0.07, which is one of the two principal parameters of the thermal model (Sect. 2). This value is low, but not negligible. Moreover, its analogues associated with higher Fourier terms ${\mathrm{\Theta }}_n=\sqrt{n}\mathrm{\Theta }$ become appreciably larger. The orbital eccentricity e ≃ 0.38 is moderately high and suitable for testing our approach.

Figure 2 shows the mean temperature T_av at the surface as a function of the mean anomaly ℓ during a revolution about the Sun for the nominal value of Θ and four selected values of the orbital eccentricity e (including e = 0.38 of Didymos). We tested the effect of series truncation by different values of maximum degree n_max in (34), comparing the results for n_max increasing from 2 to 8. Clearly, restricting to n_max < 4 is insufficient to express the temperature variation for Didymos orbit, while n_max = 8 seems to be a good compromise between precision and complexity. This is highlighted in Fig. 3, where we represent the differences between consecutive truncation degrees and compare them to the accuracy of the Didymos’ temperature measurements that will be performed by the Thermal Infrared Imager (TIRI) on board the Hera mission, according to Okada et al. (2025). The differences become progressively smaller, and in the case of the nominal Didymos eccentricity e = 0.38, for n_max = 8 they are well below the instrumental accuracy threshold. Pushing the procedure to even higher-degree values will clearly improve the accuracy, but at the expense of an increasing number of terms to be handled. For e = 0.5, n_max = 8 is not sufficient to achieve the required accuracy. By means of an exponential fit, we extrapolated a rough estimate for the degree to which we needed to push the procedure to cross the TIRI accuracy threshold, which in this case was n_max = 11.

As mentioned in Section 2, the iterative procedure that we developed relies on the convergence of the series expansion of the Hansen coefficients in terms of the eccentricity e (or β0). Therefore, in the tests performed in this section, we maintained e < e_L. This constraint does not drastically limit the scope of our intended application because e < e_L for most asteroids.

Next, we considered the dependence of the results on the two principal parameters of the thermal model identified in Sect. 2 after transforming all variables to their non-dimensional scaled images, namely the orbital eccentricity e and the thermal parameter Θ. The results are shown in Figs. 4 and 5. The former indicates that the analytical representation of T_av with the series truncation at n_max = 8 is adequate for e ≤ 0.4. This is confirmed by Fig. 6, where we represent the differences between consecutive truncation degrees with varying eccentricity (left). The dashed line represents the predicted accuracy of TIRI. By taking this accuracy as a threshold value, we see that the difference between degrees 8 and 7 might be acceptable up to e = 0.45, but a more conservative value of e = 0.4 is considered. The right panel of Fig. 6 shows the differences between the analytically computed local equilibrium temperature and the truncated average temperature T_av,n obtained for Θ = 0 and varying e. Since we can express T_eq analytically, T_eq − T_av,n is the remainder, obtained as a result of the truncation. By requiring that the remainder is smaller than a user-defined threshold, a suitable value for n_max can be inferred. This choice is motivated by the fact that the remainder in the case Θ = 0 is the only case that can be treated analytically and by continuity provides a good candidate for n_max when Θ is small. For higher eccentricity values, an extension to larger n_max would be needed. This is a matter of more algebraic work, which we postpone to future work. For the dependence on Θ, Fig. 5 illustrates how the higher thermal inertia drives the T_av variation over the revolution cycle away from the local equilibrium value (shown by the dashed black line) in two aspects: (i) first, the temperature difference between perihelion and aphelion decreases for higher Θ values (note especially the drop in the perihelion value), and (ii) the temperature maximum and minimum postdates perihelion and aphelion location. The latter is also demonstrated in Fig. 7, where only n_max = 8 solutions are compared for different Θ values. The minimum of T_av lags longer past aphelion than the maximum past perihelion. This is because of the lower aphelion temperature. In the Θ → 0 limit, T_av is expected to be identical to the locally equilibrium temperature at any instant along the orbit. This is true when n_max → ∞, but truncating the series can cause a discrepancy. If this is too large, as is the case for the lowest values of n_max, T_av would not be suitable for our future analytical modelling of the Yarkovsky effect.

Finally, we address the absolute accuracy of our defined T_av average temperature. To do this, we used the fully numerical model of 1D thermal diffusion developed by Čapek & Vokrouhlický (2004) and Čapek & Vokrouhlický (2005). The non-linear surface boundary condition (2) was solved iteratively for each surface facet of an irregularly shaped body. In our case, we used a regular polyhedron with 2004 facets closely representing a spherical body and the Didymos heliocentric orbit listed in Table 1. For the sake of definiteness, we used an 8 hr rotation period and spin axis perpendicular to the orbital plane (but we verified that the results do not depend on these values). The thermal capacity C = 600 J kg ⁻¹ K ⁻¹ and the surface density 2170 kg m⁻³ were fixed, and the thermal conductivity value K was varied to correspond to the thermal parameter Θ = 0.033, 0.07 and 0.175. At convergence, the code provides the temperature of all surface elements at any epoch along the heliocentric orbit (we used a one-minute integration timestep, but output the results every 20 minutes). For a given mean anomaly value ℓ, we numerically computed the mean value of the nth power of the surface temperature, $$\overline{T^n}(\ell )=\frac{1}{4\pi }\int d\mathrm{\Omega }{T}^n(\theta ,\varphi ;\ell )$$(49)

(d Ω = sin θ dθdϕ). By definition, $\overline{T^n}$ gives a monopole term of the spherical harmonic series of Tⁿ. Of particular interest were the cases n = 1 and n = 4, namely the mean surface temperature $\overline{T}$ and the mean value of $\overline{T^4}$. In addition to Fourier series in the mean anomaly, Tⁿ(θ, ϕ; ℓ) is also generically developed in spherical harmonic series in (θ, ϕ) and $\overline{T^n}(\ell )$ in (49) represents its monopole term. For n = 4, the respective monopole is directly related to the $\alpha \overline{\mathcal{E}}$ radiation flux term on the right side of Eq. (4) used in the definition of T_av. Therefore, we expect the closest link of T_av to some effective temperature to be related to T⁴.

Figure 8 compares our analytically determined degree 8 mean temperature T_av, shown also in Fig. 7, with $\overline{T}$ (left panel) and ${\left(\overline{T^4}\right)}^{1/4}$ (right panel). As argued above, the latter agrees well, but the former is shifted. For the full temperature formulation $\overline{T}$, there is a non-linear leak of the energy term $\alpha \overline{\mathcal{E}}$ into the higher multipoles of the spherical harmonic series, especially for low Θ values. $\overline{T}$ is therefore somewhat lower than T_av. In the limit of very high Θ values (unphysical situations, however) when the conductive part of Eq. (2) takes a larger share, T_av approaches both $\overline{T}$ and ${\left(\overline{T^4}\right)}^{1/4}$.

In conclusion, since we avoided modelling diurnal cycles at this stage, our formulation of T_av necessarily represents a compromise between full-fledged reality of T (θ, ϕ; ℓ) and limits of the analytical description. In particular, it is typically close to the mean ${\left(\overline{T^4}\right)}^{1/4}$, but for low surface thermal inertia tends to deviate from $\overline{T}$. In our intended second paper of this series, we develop a linear model for T (θ, ϕ; ℓ) and provide tests against the full numerical model. A particular goal is the ability to analytically describe the secular perturbation of the orbital semi-major axis and eccentricity.

Fig. 2 Surface mean temperature Tav (ordinate) for a spherical analogue of (65803) Didymos, computed over one revolution about the Sun (mean anomaly ℓ at the abscissa). The orbital and physical parameters are taken from Table 1, and the coefficients of Tav up to degree nmax = 8 are obtained using our algorithm. The dashed black line shows the local equilibrium temperature Teq given by Eq. (46).

Fig. 3 Differences between consecutive degrees of approximation for four different values of e. The dashed line represents the accuracy of TIRI, and the solid lines are obtained through exponential fits that can be used to extrapolate an estimate of the required truncation degree for which the differences fall below the accuracy threshold.

Fig. 4 Similar to Fig. 2, but now the orbital eccentricity e is varied for the sake of illustration (all other parameters, including the semi-major axis, have the nominal values): (i) e = 0.1 (top and left), (ii) e = 0.3 (top and right), (iii) e = 0.4 (bottom and left), and (iv) e = 0.5 (bottom and right).

Fig. 5 Similar to Fig. 2, but now the thermal parameter Θ is varied for the sake of illustration (the nominal value e = 0.38 for the eccentricity is used here): (i) Θ = 0 (top and left), (ii) Θ = 0.033 (top and right), (iii) Θ = 0.07 (bottom and left), and (iv) Θ = 0.175 (bottom and right). The Θ = 0 solution for Tav is very close to the local equilibrium situation from (46), and increasing Θ values moves Tav away from equilibrium due to finite thermal relaxation.

Fig. 6 Differences between consecutive degrees of approximation for varying e (left) and differences between the analytically computed local equilibrium temperature and the truncated average temperature obtained for Θ = 0 and varying e (right). The dashed lines represent the accuracy of TIRI, and the solid lines are obtained through exponential fits.

Fig. 7 Degree 8 mean temperature profiles for e = 0.38 and varied values of Θ (solid lines). The open circles indicate maxima (red) and minima (green) of Tav. The nested plot shows a zoom of the post-perihelion temperature profile. The higher Θ values result in a smaller difference between maximum and minimum Tav and a longer post-perihelion and -aphelion lag.

Fig. 8 Comparison between our formulation of the mean temperature Tav (degree 8 truncated series; dashed lines) and the mean value of the first (left panels) and fourth power (right panels) of the numerically determined temperature T(θ, ϕ; ℓ) defined by (49) (solid lines). The spherical body on a Didymos heliocentric orbit assumed and the four values of the thermal parameter used are Θ = 0 (black), Θ = 0.033 (blue), Θ = 0.07 (green), and Θ = 0.175 (red).

5 Conclusions

We developed an analytical framework for determining the mean temperature T_av of a spherical body moving in an elliptic orbit about the Sun. The method uses power-series expansion in two parameters ζ = exp(ιℓ), which is equivalent to Fourier series in the mean anomaly ℓ, and $\beta =e/(1+\sqrt{1-e^2})$, where e is the orbital eccentricity. The resulting algorithm is well-behaved and was tested here to degree and order eight in the series expansion. This truncation allowed us to evaluate T_av for an eccentricity up to ≃0.4. In order to extend this range to higher eccentricity values, the approach can be efficiently extended to higher degrees and orders with no need to recompute the necessary coefficients determined so far. However, since the expansions used in our method use powers in eccentricity (or β), we expect a loss of convergence at the Laplace limit of e ~ 0.66.

A possible future application of our results includes an analytical theory implementing seasonal and diurnal Yarkovsky effects, and their mutual interaction, for small Solar System bodies. With the mean temperature properly defined for eccentric orbits, this will generalize the results given by Vokrouhlický (1999), where the source of the coupling between the seasonal and diurnal Yarkovsky effects originated from a finite ratio of the revolution period about the Sun and the rotation period about the spin axis. Having formulated the mean temperature along an eccentric orbit, we will be able to analyse the interaction between the two variants of the Yarkovsky effect originating from the orbital eccentricity using an analytical approach.

Acknowledgements

RP and MM are grateful to the Italian Space Agency (ASI) for financial support through Agreement No. 2022-8-HH.0 in the context of ESA’s Hera mission. RP also acknowledges ASI’s support through the project “MONitoring ASTERoids” (grant 2022-33-HH.0). The work of DV was partially supported by the Czech Science Foudation (grant 25-16507S).

References

Appendix A Symbols and notation

Constants and parameters:
nmax	maximum degree of truncation
σ	Stefan–Boltzmann constant
ρ	density
C	thermal capacity
K	thermal conductivity
Γ	thermal inertia
α	absorption coefficient
ϵ	emissivity
R	effective radius of the body
Θ	thermal parameter
Θn	auxiliary thermal parameter defined by ${\mathrm{\Theta }}_n=\sqrt{n}\mathrm{\Theta }$, with n = 1, ..., nmax
hs	penetration depth of the seasonal thermal wave
n	mean motion
a	semimajor axis
e	eccentricity
eL	Laplace limit
β	$\beta =e/(1+\sqrt{1-e^2})$
Fa	solar radiation flux at distance equal to the semimajor axis
Ta	reference temperature defined by $ϵ\sigma T_a^4=\alpha F_a$
Variables:
r	radial position with respect to the body centre
r′	radial position scaled by hs
d	distance from the Sun
t	time
ℓ	mean anomaly
ζ	time-related variable, ζ = eiℓ
T	temperature
Tav	mean temperature
$T_{av}^{\mathrm{\prime }}$	mean temperature scaled by Ta
Teq	local equilibrium temperature obtained from $T_{\text{av\hspace{0.17em}}}^{\mathrm{\prime }4}={\overline{\mathcal{E}}}^{\mathrm{\prime }}$
ΔT	temperature variation with respect to the average value
ℰ	insolation
$\overline{T}$	mean surface temperature
${\overline{\mathcal{E}}}^{\mathrm{\prime }}$	average insolation scaled by the reference flux Fa
xn	$x_n=\sqrt{n/2}{R}^{\mathrm{\prime }}$
zn	$z_n=\sqrt{-\iota n}{r}^{\mathrm{\prime }}$
Coefficients:
an	Hansen coefficients
Cn	coefficients of the expansion in terms of ζ
Cnk	coefficients of the expansion of Cn in terms of β
Functions:
j0	zeroth spherical Bessel function of the first kind
ψn	function defined by (z/j0(z)) dz(j0(z)|z=Zn
$G_n^{+},G_n^{-},F_n$	auxiliary functions

Appendix B Coefficients C_nk up to degree and order eight

Here we provide explicit formulas for coefficients C_nk up to degree and order eight for the large-body limit (R′ → ∞). We introduce the following quantities $$\begin{array}{rl}{\mathrm{\Theta }}_n=& \sqrt{n}\mathrm{\Theta }\text{\hspace{0.17em}and\hspace{0.17em}}{Z}_n=\frac{1+{\mathrm{\Theta }}_n/2-i{\mathrm{\Theta }}_n/2}{\sqrt{2}(1+{\mathrm{\Theta }}_n+{\mathrm{\Theta }}_n^2/2),}\text{\hspace{0.17em}with which we obtain:\hspace{0.17em}}\\ C_{11}=& Z_1/2C_{02}=(12C_{-11}{C}_{11}-1)/(2\sqrt{2})C_{22}=Z_2(125-3C_{11}^2)C_{33}=Z_3(13-8\sqrt{2}{C}_{11}^3-24C_{11}{C}_{22})/4\\ C_{13}=& Z_1[1-24C_{-11}(\sqrt{2}{C}_{11}^2+C_{22})-24C_{11}{C}_{02}]/4C_{44}=Z_4(103-12C_{11}^4-72\sqrt{2}{C}_{11}^2{C}_{22}-36C_{22}^2-72C_{11}{C}_{33})/12\\ C_{24}=& -Z_2[11+36C_{11}{C}_{13}+12C_{-11}(2C_{11}^3+6\sqrt{2}{C}_{11}{C}_{22}+3C_{33})+36\sqrt{2}{C}_{11}^2{C}_{02}+36C_{22}{C}_{02}]/6\\ C_{04}=& (1-12C_{-13}{C}_{11}-12\sqrt{2}{C}_{-22}{C}_{11}^2-12C_{-11}^2{C}_{11}^2-12C_{-11}{C}_{13}-12C_{-22}{C}_{22}-12\sqrt{2}{C}_{-11}^2{C}_{22}-24\sqrt{2}{C}_{-11}{C}_{11}{C}_{02}-6C_{02}^2)/(2\sqrt{2})\\ C_{55}=& Z_5(1097-192C_{11}^3{C}_{22}-288\sqrt{2}{C}_{11}^2{C}_{33}-288C_{22}{C}_{33}-288C_{11}(\sqrt{2}{C}_{22}^2+C_{44}))Z_5)/48\\ C_{35}=& -Z_3\{181+96C_{13}{C}_{22}+96\sqrt{2}{C}_{-11}{C}_{22}^2+96C_{11}^2(\sqrt{2}{C}_{13}+2C_{-11}{C}_{22})+96C_{-11}{C}_{44}+64C_{11}^3{C}_{02}+96C_{33}{C}_{02}\\ & +96C_{11}[C_{24}+2\sqrt{2}(C_{-11}{C}_{33}+C_{22}{C}_{02})]\}/16\\ C_{15}=& Z_1[23-288\sqrt{2}{C}_{-11}{C}_{11}{C}_{13}-288C_{-11}^2{C}_{11}{C}_{22}-144C_{-13}(\sqrt{2}{C}_{11}^2+C_{22})-144C_{-11}{C}_{24}-144\sqrt{2}{C}_{-11}^2{C}_{33}\\ & -48C_{-22}(2C_{11}^3+6\sqrt{2}{C}_{11}{C}_{22}+3C_{33})-288C_{-11}{C}_{11}^2{C}_{02}-144C_{13}{C}_{02}-288\sqrt{2}{C}_{-11}{C}_{22}{C}_{02}-144\sqrt{2}{C}_{11}{C}_{02}^2-144C_{11}{C}_{04}]/24\\ C_{66}=& Z_6[1223-40\sqrt{2}{C}_{22}^3-80C_{11}^3{C}_{33}-60C_{33}^2-120C_{22}{C}_{44}-120C_{11}^2(C_{22}^2+\sqrt{2}{C}_{44})-120C_{11}(2\sqrt{2}{C}_{22}{C}_{33}+C_{55})]/20\\ C_{46}=& -Z_4\{1417+120C_{11}^3{C}_{13}+180C_{13}{C}_{33}+180C_{22}(C_{24}+2\sqrt{2}{C}_{-11}{C}_{33})+180C_{-11}{C}_{55}+180\sqrt{2}{C}_{22}^2{C}_{02}+180C_{44}{C}_{02}\\ & +180C_{11}^2(\sqrt{2}{C}_{24}+2C_{-11}{C}_{33}+2C_{22}{C}_{02})+180C_{11}[2\sqrt{2}{C}_{13}{C}_{22}+C_{35}+2C_{-11}(C_{22}^2+\sqrt{2}{C}_{44})+2\sqrt{2}{C}_{33}{C}_{02}]\}/30\\ C_{26}=& Z_2\{19-9C_{13}^2-18C_{11}{C}_{15}-36C_{-22}{C}_{11}^2{C}_{22}-18\sqrt{2}{C}_{-22}{C}_{22}^2-36\sqrt{2}{C}_{-22}{C}_{11}{C}_{33}-6C_{-13}(2C_{11}^3+6\sqrt{2}{C}_{11}{C}_{22}+3C_{33})\\ & -18C_{-22}{C}_{44}-18C_{-11}^2(C_{22}^2+2C_{11}{C}_{33}+\sqrt{2}{C}_{44})-36\sqrt{2}{C}_{11}{C}_{13}{C}_{02}-18C_{24}{C}_{02}-18C_{11}^2{C}_{02}^2-18\sqrt{2}{C}_{22}{C}_{02}^2\\ & -18C_{-11}[2C_{11}^2{C}_{13}+2\sqrt{2}{C}_{13}{C}_{22}+C_{35}+2\sqrt{2}{C}_{33}{C}_{02}+2C_{11}(\sqrt{2}{C}_{24}+2C_{22}{C}_{02})]-18\sqrt{2}{C}_{11}^2{C}_{04}-18C_{22}{C}_{04}\}/3\\ C_{06}=& (1-12C_{-15}{C}_{11}-12\sqrt{2}{C}_{-24}{C}_{11}^2-24C_{-11}{C}_{-13}{C}_{11}^2-8C_{-33}{C}_{11}^3-12C_{-13}{C}_{13}-24\sqrt{2}{C}_{-22}{C}_{11}{C}_{13}-24C_{-11}^2{C}_{11}{C}_{13}\\ & -12C_{-11}{C}_{15}-12C_{-24}{C}_{22}-24\sqrt{2}{C}_{-11}{C}_{-13}{C}_{22}-24\sqrt{2}{C}_{-33}{C}_{11}{C}_{22}-48C_{-22}{C}_{-11}{C}_{11}{C}_{22}-12C_{-22}{C}_{24}-12\sqrt{2}{C}_{-11}^2{C}_{24}\\ & -12C_{-33}{C}_{33}+24\sqrt{2}{C}_{-22}{C}_{-11}{C}_{33}-8C_{-11}^3{C}_{33}-24\sqrt{2}{C}_{-13}{C}_{11}{C}_{02}-24C_{-22}{C}_{11}^2{C}_{02}-24\sqrt{2}{C}_{-11}{C}_{13}{C}_{02}-24\sqrt{2}{C}_{-22}{C}_{22}{C}_{02}\\ & -24C_{-11}^2{C}_{22}{C}_{02}-24C_{-11}{C}_{11}{C}_{02}^2-4\sqrt{2}{C}_{02}^3-24\sqrt{2}{C}_{-11}{C}_{11}{C}_{04}-12C_{02}{C}_{04})/(2\sqrt{2})\\ C_{77}=& Z_7\{47273/288-6\sqrt{2}{C}_{22}^2{C}_{33}-4C_{11}^3{C}_{44}-6C_{33}{C}_{44}-6C_{22}{C}_{55}\\ & -6C_{11}^2(2C_{22}{C}_{33}+\sqrt{2}{C}_{55})-2C_{11}[2C_{22}^3+6\sqrt{2}{C}_{22}{C}_{44}+3(\sqrt{2}{C}_{33}^2+C_{66})]\}\\ C_{57}=& -Z_5\{49531+1152C_{-11}{C}_{22}^3+1152C_{11}^3{C}_{24}+1728C_{24}{C}_{33}+1728\sqrt{2}{C}_{-11}{C}_{33}^2+1728C_{22}{C}_{35}+3456\sqrt{2}{C}_{-11}{C}_{22}{C}_{44}\\ & +1728C_{13}(\sqrt{2}{C}_{22}^2+C_{44})+1728C_{-11}{C}_{66}+3456\sqrt{2}{C}_{22}{C}_{33}{C}_{02}+1728C_{55}{C}_{02}+1728C_{11}^2(2C_{13}{C}_{22}+\sqrt{2}{C}_{35}+2C_{-11}{C}_{44}\\ & +2C_{33}{C}_{02})+1728C_{11}[2\sqrt{2}{C}_{13}{C}_{33}+2C_{22}(\sqrt{2}{C}_{24}+2C_{-11}{C}_{33})+C_{46}+2\sqrt{2}{C}_{-11}{C}_{55}+2C_{22}^2{C}_{02}+2\sqrt{2}{C}_{44}{C}_{02}]\}/288\\ C_{37}=& Z_3\{1267-192C_{15}{C}_{22}-192\sqrt{2}{C}_{-13}{C}_{22}^2-192C_{13}{C}_{24}-384\sqrt{2}{C}_{-11}{C}_{22}{C}_{24}-384\sqrt{2}{C}_{-11}{C}_{13}{C}_{33}-384\sqrt{2}{C}_{-22}{C}_{22}{C}_{33}\\ & -384C_{-11}^2{C}_{22}{C}_{33}-192C_{-13}{C}_{44}-192C_{-11}{C}_{46}-192C_{-22}{C}_{55}-192\sqrt{2}{C}_{-11}^2{C}_{55}-384\sqrt{2}{C}_{13}{C}_{22}{C}_{02}-384C_{-11}{C}_{22}^2{C}_{02}\\ & -192C_{35}{C}_{02}-384\sqrt{2}{C}_{-11}{C}_{44}{C}_{02}-192\sqrt{2}{C}_{33}{C}_{02}^2-192C_{11}^2[\sqrt{2}{C}_{15}+2(C_{-13}{C}_{22}+C_{-11}{C}_{24}+C_{-22}{C}_{33}+C_{13}{C}_{02})]\\ & -128C_{11}^3{C}_{04}-192C_{33}{C}_{04}-192C_{11}[\sqrt{2}{C}_{13}^2+4C_{-11}{C}_{13}{C}_{22}+C_{26}+2\sqrt{2}{C}_{-13}{C}_{33}+2\sqrt{2}{C}_{-11}{C}_{35}+2C_{-11}^2{C}_{44}\\ & +2C_{-22}(C_{22}^2+\sqrt{2}{C}_{44})+2\sqrt{2}{C}_{24}{C}_{02}+4C_{-11}{C}_{33}{C}_{02}+2C_{22}{C}_{02}^2+2\sqrt{2}{C}_{22}{C}_{04}\left]\right\}/32\\ C_{17}=& -Z_1(73+3456\sqrt{2}{C}_{-13}{C}_{11}{C}_{13}+3456C_{-22}{C}_{11}^2{C}_{13}+1728\sqrt{2}{C}_{-11}{C}_{13}^2+3456\sqrt{2}{C}_{-11}{C}_{11}{C}_{15}\\ & +6912C_{-11}{C}_{-13}{C}_{11}{C}_{22}+3456C_{-33}{C}_{11}^2{C}_{22}+3456\sqrt{2}{C}_{-22}{C}_{13}{C}_{22}+3456C_{-11}^2{C}_{13}{C}_{22}+1728\sqrt{2}{C}_{-33}{C}_{22}^2\\ & +3456C_{-22}{C}_{-11}{C}_{22}^2+1728C_{-15}(\sqrt{2}{C}_{11}^2+C_{22})+1728C_{-13}{C}_{24}+3456\sqrt{2}{C}_{-22}{C}_{11}{C}_{24}+3456C_{-11}^2{C}_{11}{C}_{24}\\ & +1728C_{-11}{C}_{26}+3456\sqrt{2}{C}_{-11}{C}_{-13}{C}_{33}+3456\sqrt{2}{C}_{-33}{C}_{11}{C}_{33}+6912C_{-22}{C}_{-11}{C}_{11}{C}_{33}+576C_{-24}(2C_{11}^3+6\sqrt{2}{C}_{11}{C}_{22}\\ & +3C_{33})+1728C_{-22}{C}_{35}+1728\sqrt{2}{C}_{-11}^2{C}_{35}+1728C_{-33}{C}_{44}+3456\sqrt{2}{C}_{-22}{C}_{-11}{C}_{44}+1152C_{-11}^3{C}_{44}+3456C_{-13}{C}_{11}^2{C}_{02}\\ & +6912C_{-11}{C}_{11}{C}_{13}{C}_{02}+1728C_{15}{C}_{02}+3456\sqrt{2}{C}_{-13}{C}_{22}{C}_{02}+6912C_{-22}{C}_{11}{C}_{22}{C}_{02}+3456\sqrt{2}{C}_{-11}{C}_{24}{C}_{02}\\ & +3456\sqrt{2}{C}_{-22}{C}_{33}{C}_{02}+3456C_{-11}^2{C}_{33}{C}_{02}+1728\sqrt{2}{C}_{13}{C}_{02}^2+3456C_{-11}{C}_{22}{C}_{02}^2+1152C_{11}{C}_{02}^3+3456C_{-11}{C}_{11}^2{C}_{04}\\ & +1728C_{13}{C}_{04}+3456\sqrt{2}{C}_{-11}{C}_{22}{C}_{04}+3456\sqrt{2}{C}_{11}{C}_{02}{C}_{04}+1728C_{11}{C}_{06})/288\\ C_{88}=& Z_8[556403/1260-C_{22}^4-3C_{44}^2-6C_{22}^2(2C_{11}{C}_{33}+\sqrt{2}{C}_{44})-4C_{11}^3{C}_{55}-6C_{33}{C}_{55}-6C_{22}(\sqrt{2}{C}_{33}^2+2C_{11}^2{C}_{44}\\ & +2\sqrt{2}{C}_{11}{C}_{55}+C_{66})-6C_{11}^2(C_{33}^2+\sqrt{2}{C}_{66})-6C_{11}(2\sqrt{2}{C}_{33}{C}_{44}+C_{77})]\\ C_{68}=& -Z_6\{40861+420C_{22}^2(\sqrt{2}{C}_{24}+2C_{-11}{C}_{33})+280C_{11}^3{C}_{35}+420C_{33}{C}_{35}+420C_{24}{C}_{44}+840\sqrt{2}{C}_{-11}{C}_{33}{C}_{44}\\ & +420C_{13}{C}_{55}+420C_{-11}{C}_{77}+280C_{22}^3{C}_{02}+420\sqrt{2}{C}_{33}^2{C}_{02}+420C_{66}{C}_{02}+420C_{11}^2(2C_{22}{C}_{24}+2C_{13}{C}_{33}+\sqrt{2}{C}_{46}\\ & +2C_{-11}{C}_{55}+2C_{44}{C}_{02})+420C_{11}[2\sqrt{2}{C}_{24}{C}_{33}+2C_{-11}{C}_{33}^2+2\sqrt{2}{C}_{22}{C}_{35}+4C_{-11}{C}_{22}{C}_{44}+2C_{13}(C_{22}^2+\sqrt{2}{C}_{44})+C_{57}\\ & +2\sqrt{2}{C}_{-11}{C}_{66}+4C_{22}{C}_{33}{C}_{02}+2\sqrt{2}{C}_{55}{C}_{02}]+420C_{22}[2\sqrt{2}{C}_{13}{C}_{33}+C_{46}+2\sqrt{2}(C_{-11}{C}_{55}+C_{44}{C}_{02})]\}/70\\ C_{48}=& Z_4\{18037-360C_{11}^3{C}_{15}-540\sqrt{2}{C}_{13}^2{C}_{22}-360C_{-22}{C}_{22}^3-270C_{24}^2-540C_{22}{C}_{26}-540C_{15}{C}_{33}-1080\sqrt{2}{C}_{-13}{C}_{22}{C}_{33}\\ & -1080\sqrt{2}{C}_{-11}{C}_{24}{C}_{33}-540\sqrt{2}{C}_{-22}{C}_{33}^2-540C_{-11}^2{C}_{33}^2-1080\sqrt{2}{C}_{-11}{C}_{22}{C}_{35}-1080\sqrt{2}{C}_{-22}{C}_{22}{C}_{44}-1080C_{-11}^2{C}_{22}{C}_{44}\\ & -540C_{-13}{C}_{55}-540C_{-11}{C}_{57}-540C_{-22}{C}_{66}-540\sqrt{2}{C}_{-11}^2{C}_{66}-1080\sqrt{2}{C}_{22}{C}_{24}{C}_{02}-2160C_{-11}{C}_{22}{C}_{33}{C}_{02}-540C_{46}{C}_{02}\\ & -1080\sqrt{2}{C}_{-11}{C}_{55}{C}_{02}-540C_{22}^2{C}_{02}^2-540\sqrt{2}{C}_{44}{C}_{02}^2-540C_{13}[C_{35}+2C_{-11}(C_{22}^2+\sqrt{2}{C}_{44})+2\sqrt{2}{C}_{33}{C}_{02}]-540\sqrt{2}{C}_{22}^2{C}_{04}\\ & -540C_{44}{C}_{04}-540C_{11}[2\sqrt{2}{C}_{15}{C}_{22}+2\sqrt{2}{C}_{13}{C}_{24}+4C_{-11}{C}_{22}{C}_{24}+4C_{-11}{C}_{13}{C}_{33}+4C_{-22}{C}_{22}{C}_{33}+C_{37}+2C_{-13}(C_{22}^2\\ & +\sqrt{2}{C}_{44})+2\sqrt{2}{C}_{-11}{C}_{46}+2\sqrt{2}{C}_{-22}{C}_{55}+2C_{-11}^2{C}_{55}+4C_{13}{C}_{22}{C}_{02}+2\sqrt{2}{C}_{35}{C}_{02}+4C_{-11}{C}_{44}{C}_{02}+2C_{33}{C}_{02}^2+2\sqrt{2}{C}_{33}{C}_{04}]\\ & -540C_{11}^2[C_{13}^2+\sqrt{2}{C}_{26}+2(C_{-13}{C}_{33}+C_{-11}{C}_{35}+C_{-22}{C}_{44}+C_{24}{C}_{02}+C_{22}{C}_{04})]\}/90\\ C_{28}=& -Z_2\{1127+1080C_{-11}{C}_{11}{C}_{13}^2+1080C_{-11}{C}_{11}^2{C}_{15}+40C_{13}{C}_{15}+540C_{11}{C}_{17}+1080C_{-24}{C}_{11}^2{C}_{22}\\ & +160C_{-22}{C}_{11}{C}_{13}{C}_{22}+1080\sqrt{2}{C}_{-11}{C}_{15}{C}_{22}+540\sqrt{2}{C}_{-24}{C}_{22}^2+1080C_{-33}{C}_{11}{C}_{22}^2+1080C_{-22}{C}_{11}^2{C}_{24}\\ & +1080\sqrt{2}{C}_{-11}{C}_{13}{C}_{24}+1080\sqrt{2}{C}_{-22}{C}_{22}{C}_{24}+1080C_{-11}^2{C}_{22}{C}_{24}+1080\sqrt{2}{C}_{-11}{C}_{11}{C}_{26}+1080\sqrt{2}{C}_{-24}{C}_{11}{C}_{33}\\ & +1080C_{-33}{C}_{11}^2{C}_{33}+1080\sqrt{2}{C}_{-22}{C}_{13}{C}_{33}+1080C_{-11}^2{C}_{13}{C}_{33}+1080\sqrt{2}{C}_{-33}{C}_{22}{C}_{33}+2160C_{-22}{C}_{-11}{C}_{22}{C}_{33}\\ & +180C_{-15}(2C_{11}^3+6\sqrt{2}{C}_{11}{C}_{22}+3C_{33})+1080\sqrt{2}{C}_{-22}{C}_{11}{C}_{35}+1080C_{-11}^2{C}_{11}{C}_{35}+540C_{-11}{C}_{37}+540C_{-24}{C}_{44}\\ & +1080\sqrt{2}{C}_{-33}{C}_{11}{C}_{44}+2160C_{-22}{C}_{-11}{C}_{11}{C}_{44}+540C_{-22}{C}_{46}+540\sqrt{2}{C}_{-11}^2{C}_{46}+540C_{-33}{C}_{55}+1080\sqrt{2}{C}_{-22}{C}_{-11}{C}_{55}\\ & +360C_{-11}^3{C}_{55}+540\sqrt{2}{C}_{13}^2{C}_{02}+1080\sqrt{2}{C}_{11}{C}_{15}{C}_{02}+2160C_{-11}{C}_{13}{C}_{22}{C}_{02}+1080C_{-22}{C}_{22}^2{C}_{02}+2160C_{-11}{C}_{11}{C}_{24}{C}_{02}\\ & +540C_{26}{C}_{02}+2160C_{-22}{C}_{11}{C}_{33}{C}_{02}+1080\sqrt{2}{C}_{-11}{C}_{35}{C}_{02}+1080\sqrt{2}{C}_{-22}{C}_{44}{C}_{02}+1080C_{-11}^2{C}_{44}{C}_{02}+1080C_{11}{C}_{13}{C}_{02}^2\\ & +540\sqrt{2}{C}_{24}{C}_{02}^2+1080C_{-11}{C}_{33}{C}_{02}^2+360C_{22}{C}_{02}^3+540C_{-13}[2C_{11}^2{C}_{13}+2\sqrt{2}{C}_{13}{C}_{22}+2C_{-11}{C}_{22}^2+C_{35}+2\sqrt{2}{C}_{-11}{C}_{44}\\ & +2\sqrt{2}{C}_{33}{C}_{02}+2C_{11}(\sqrt{2}{C}_{24}+2C_{-11}{C}_{33}+2C_{22}{C}_{02})]+1080\sqrt{2}{C}_{11}{C}_{13}{C}_{04}+2160C_{-11}{C}_{11}{C}_{22}{C}_{04}+540C_{24}{C}_{04}\\ & +1080\sqrt{2}{C}_{-11}{C}_{33}{C}_{04}+1080C_{11}^2{C}_{02}{C}_{04}+1080\sqrt{2}{C}_{22}{C}_{02}{C}_{04}+540\sqrt{2}{C}_{11}^2{C}_{06}+540C_{22}{C}_{06}\}/90\\ C_{08}=& (1-12C_{-17}{C}_{11}-12\sqrt{2}{C}_{-26}{C}_{11}^2-12C_{-13}^2{C}_{11}^2-24C_{-11}{C}_{-15}{C}_{11}^2-8C_{-35}{C}_{11}^3-12C_{-15}{C}_{13}-24\sqrt{2}{C}_{-24}{C}_{11}{C}_{13}\\ & -48C_{-11}{C}_{-13}{C}_{11}{C}_{13}-24C_{-33}{C}_{11}^2{C}_{13}-12\sqrt{2}{C}_{-22}{C}_{13}^2-12C_{-11}^2{C}_{13}^2-12C_{-13}{C}_{15}-24\sqrt{2}{C}_{-22}{C}_{11}{C}_{15}-24C_{-11}^2{C}_{11}{C}_{15}\\ & -12C_{-11}{C}_{17}-12C_{-26}{C}_{22}-12\sqrt{2}{C}_{-13}^2{C}_{22}-24\sqrt{2}{C}_{-11}{C}_{-15}{C}_{22}-24\sqrt{2}{C}_{-35}{C}_{11}{C}_{22}-48C_{-24}{C}_{-11}{C}_{11}{C}_{22}\\ & -48C_{-22}{C}_{-13}{C}_{11}{C}_{22}-24C_{-44}{C}_{11}^2{C}_{22}-4\sqrt{2}{C}_{-33}{C}_{13}{C}_{22}-48C_{-22}{C}_{-11}{C}_{13}{C}_{22}-12\sqrt{2}{C}_{-44}{C}_{22}^2-12C_{-22}^2{C}_{22}^2-24C_{-33}{C}_{-11}{C}_{22}^2\\ & -12C_{-24}{C}_{24}-24\sqrt{2}{C}_{-11}{C}_{-13}{C}_{24}-24\sqrt{2}{C}_{-33}{C}_{11}{C}_{24}-48C_{-22}{C}_{-11}{C}_{11}{C}_{24}-12C_{-22}{C}_{26}-12\sqrt{2}{C}_{-11}^2{C}_{26}-12C_{-35}{C}_{33}\\ & -24\sqrt{2}{C}_{-24}{C}_{-11}{C}_{33}-24\sqrt{2}{C}_{-22}{C}_{-13}{C}_{33}-24C_{-11}^2{C}_{-13}{C}_{33}-24\sqrt{2}{C}_{-44}{C}_{11}{C}_{33}-24C_{-22}^2{C}_{11}{C}_{33}-48C_{-33}{C}_{-11}{C}_{11}{C}_{33}\\ & -12C_{-33}{C}_{35}-24\sqrt{2}{C}_{-22}{C}_{-11}{C}_{35}-8C_{-11}^3{C}_{35}-12C_{-44}{C}_{44}-12\sqrt{2}{C}_{-22}^2{C}_{44}-24\sqrt{2}{C}_{-33}{C}_{-11}{C}_{44}-24C_{-22}{C}_{-11}^2{C}_{44}\\ & -24\sqrt{2}{C}_{-15}{C}_{11}{C}_{02}-24C_{-24}{C}_{11}^2{C}_{02}-24\sqrt{2}{C}_{-13}{C}_{13}{C}_{02}-48C_{-22}{C}_{11}{C}_{13}{C}_{02}-24\sqrt{2}{C}_{-11}{C}_{15}{C}_{02}-24\sqrt{2}{C}_{-24}{C}_{22}{C}_{02}\\ & -48C_{-11}{C}_{-13}{C}_{22}{C}_{02}-48C_{-33}{C}_{11}{C}_{22}{C}_{02}-24\sqrt{2}{C}_{-22}{C}_{24}{C}_{02}-24C_{-11}^2{C}_{24}{C}_{02}-24\sqrt{2}{C}_{-33}{C}_{33}{C}_{02}-48C_{-22}{C}_{-11}{C}_{33}{C}_{02}\\ & -24C_{-13}{C}_{11}{C}_{02}^2-24C_{-11}{C}_{13}{C}_{02}^2-24C_{-22}{C}_{22}{C}_{02}^2-2C_{02}^4-24\sqrt{2}{C}_{-13}{C}_{11}{C}_{04}-24C_{-22}{C}_{11}^2{C}_{04}-24\sqrt{2}{C}_{-11}{C}_{13}{C}_{04}\\ & -24\sqrt{2}{C}_{-22}{C}_{22}{C}_{04}-24C_{-11}^2{C}_{22}{C}_{04}-48C_{-11}{C}_{11}{C}_{02}{C}_{04}-12\sqrt{2}{C}_{02}^2{C}_{04}-6C_{04}^2-24\sqrt{2}{C}_{-11}{C}_{11}{C}_{06}-2C_{02}{C}_{06})/(2\sqrt{2})\end{array}$$