© The Authors 2026

1 Introduction

Tidal deformation and energy dissipation play an important role in the long-term orbital evolution of planets and moons. The magnitude and frequency of tidal dissipation depend on the internal structure of the deformable body, in particular on the distribution of density, elasticity, and viscosity.

Tidal deformation is modeled with the frequency-dependent degree-2 tidal Love number k₂(s), which can be obtained by solving the governing equations for a radially stratified viscoelastic body. Here, the stratified forward model follows the propagatormatrix formulation for incompressible, spherically symmetric stratified bodies in Sabadini et al. (2016, Chapter 2). Solid layers are modeled with Maxwell rheology and liquid layers through the appropriate small-shear limit. This approach is a consistent way to evaluate the tidal response of the body, but it can become too costly to solve repeatedly within long-term orbital or spin-evolution integrations.

Long-term dynamical evolution codes often rely on simplified homogeneous rheologies or empirical prescriptions, such as constant time lag, constant phase lag, or single-element Maxwell models. These approximations are computationally convenient but generally fail to reproduce the multiple relaxation timescales and dissipation peaks that arise in layered bodies Gevorgyan (2021); Bolmont et al. (2020). In Gevorgyan et al. (2023), the authors showed that, under the assumptions of a simple multilayer rheology, the tidal response of a stratified body can be modeled by a homogeneous body with a sufficiently complex generalized rheology, such as a generalized Voigt model. This motivated us to develop a relatively simple workflow that allows one to replace the stratified problem by an effective homogeneous model that preserves the frequency-dependent tidal response.

In this work, we present an open-source Wolfram Language package¹ that takes a user-specified multilayered model and constructs the corresponding effective homogeneous generalized Voigt rheology. The input file contains the layer radii, densities, shear moduli, viscosities, and forcing frequency, and the main routine performs the full reduction pipeline: computation of k₂(s) for the layered model, identification of secular poles and residues, inversion to the equivalent compliance function, and extraction of the generalized Voigt parameters. An optional reduction step identifies dominant relaxation modes over a user-defined frequency interval, yielding a compact model when an exact full representation is unnecessary.

Our implementation was partly motivated by the need to connect detailed geophysical interior models with modern time-domain orbital integrators such as the RheoVolution (de Oliveira et al. 2025a,b), which can evolve deformable bodies with complex homogeneous rheologies, but require as input the parameters of an effective rheological model rather than a full multilayer interior description. The package converts the tidal response of a stratified viscoelastic body into the parameters of a homogeneous generalized Voigt model that can be used in such simulations.

The paper is organized as follows. In Section 2, we summarize the theoretical framework for effective homogeneous rheology and its relation to the layered tidal response. In Section 3, we describe the structure of the software package, its inputs and outputs, and the computational workflow. In Section 4, we illustrate the method with an application to a five-layer model of the Moon and show how the package recovers the equivalent generalized Voigt rheology and, if desired, an optional reduced model that captures the tidal response over the frequency interval of interest. In Section 5 we conclude the paper with observations on the purpose and scope of the work and a few perspectives for improvement.

2 Theoretical background

This section summarizes the theoretical elements relevant to the software pipeline. The forward computation of the layered tidal response follows the propagator-matrix formulation for incompressible, stratified viscoelastic bodies presented in Sabadini et al. (2016), while the reduction to an equivalent homogeneous rheology is based on Gevorgyan et al. (2023).

2.1 Layered viscoelastic model

If we consider a spherically symmetric body of N incompressible, homogeneous layers, each with outer radius r_k, density ρ_k, rigidity μ_k, and viscosity η_k, solid layers are assumed to have Maxwell rheology, and their complex shear modulus in the Laplace domain is $${\hat{\mu }}_k(s)=\frac{{\mu }_k{\eta }_{k}s}{{\mu }_k+{\eta }_{k}s}.$$(1)

The tidal response is described by the degree-2 Love number k₂(s), obtained from the propagator-matrix solution of the governing equations for a stratified body. All layers are treated within a unified propagator formulation, and liquid layers are treated in the small-shear limit; see Section 3.3.

For a finite multilayer Maxwell body, k₂(s) is a rational function of the Laplace variable and admits the mode decomposition (Gevorgyan et al. 2023, Eq. (26)) $$k_2(s)=k_E+\sum _{i=1}^{n+1}\frac{r_i}{s-s_i},$$(2)

where k_E is the elastic Love number, s_i < 0 are the secular poles, r_i are the corresponding residues and ${\tau }_i=-s_i^{-1}$ are the relaxation times. Here n denotes the number of Kelvin–Voigt elements of the equivalent homogeneous model, so the layered body has n + 1 normal modes, in agreement with Gevorgyan et al. (2023, Eq. (22)). Thus, the layered tidal response is fully characterized by the set of parameters {k_E, s_i, r_i}.

Fig. 1 Schematic representation of the generalized Voigt model.

2.2 The homogeneous model

The equivalent homogeneous representation of stratified bodies is given by a generalized Voigt rheology in Figure 1. The Love number of the body is (Gevorgyan et al. 2023, Eq. (2)) $$k_2(s)=\frac{3{\mathrm{I}}_{\circ }G}{R^5}C(s),$$(3)

where I_° is the mean moment of inertia, C(s) is the complex compliance, and R is the mean radius of the body. The gravitational rigidity γ is defined from the fluid Love number k_f as (Gevorgyan et al. 2023, Eq. (4)) $$\gamma =\frac{3{\mathrm{I}}_{\circ }G}{R^5}\frac{1}{k_f}.$$(4)

For generalized Voigt rheology, complex compliance is decomposed as (Gevorgyan et al. 2023, Eq. (14)) $$C(s)=\frac{1}{\alpha +\gamma }+C_v(s),$$(5)

where α is the elastic rigidity and C_v(s) contains the viscous contribution of the dashpot and a finite number of Kelvin–Voigt elements in Figure 1. The equivalent homogeneous model is specified by the parameters {γ, α, η, α_i, η_i}. In these expressions γ, α, and α_i have dimensions of s⁻², as does the mapping constant Θ = 3I_°G/R⁵, while η and η_i have dimensions of s⁻¹. To convert to SI units of Pa and Pa s, we divide by the rescaling factor $\frac{152\pi }{15}\frac{R}{m_1}$ (Gevorgyan et al. 2020, Section 7.1.1), where m₁ is the mass of the body. The outputs reported by the package (for example, in Table 6) are given in these units.

2.3 Extraction of the homogeneous model parameters

The parameters of the equivalent homogeneous generalized Voigt rheology are obtained from the layered pole–residue data ${\{k_E,s_i,r_i\}}_{i=1}^{n+1}$ by a purely algebraic procedure. Writing the layered compliance $$C_L(s)=\frac{R^5}{3I_{\circ }G}{k}_2(s)=k_E^C+\sum _{i=1}^{n+1}\frac{r_i^C}{s-s_i},$$(6)

where $k_E^C$ and $r_i^C$ are k_E and r_i rescaled by R⁵/(3I_°G), the elastic and gravitational parameters follow from the asymptotic limits of compliance, $$\gamma ={(k_E^C-\sum _{i=1}^{n+1}\frac{r_i^C}{s_i})}^{-1},\alpha =\frac{1}{k_E^C}-\gamma .$$(7)

The remaining parameters are η and the n Kelvin–Voigt pairs ${\left\{\right({\alpha }_i,{\eta }_i\left)\right\}}_{i=1}^n$, totaling 2n + 1 unknowns. Attempting to extract these parameters from the n + 1 pole conditions of the generalized Voigt model would be both underdetermined and nonlinear. Instead, we construct the complex compliance of the viscoelastic element $\hat{J}(s)$ (Gevorgyan et al. 2023, Eq. (13)), $$\hat{J}(s)=\frac{1}{\alpha }+\frac{1}{\eta s}+\sum _{i=1}^n\frac{1}{{\alpha }_i+s{\eta }_i},$$(8)

directly from the layered data through the relation $C^{-1}(s)=\gamma +{\hat{J}}^{-1}(s)$, which yields $$\hat{J}(s)=\frac{C_L(s)}{1-\gamma C_L(s)}.$$(9)

All Kelvin–Voigt parameters are then read off from the partial-fraction decomposition of $\hat{J}(s)$.

Writing the layered compliance (6) as a single fraction C_L(s) = N(s)/D(s), where $D(s)={\prod }_{i=1}^{n+1}(s-s_i)$ and N(s) is a polynomial of degree n + 1, Equation (9) gives $$\hat{J}(s)=\frac{N(s)}{D(s)-\gamma N(s)}=\frac{N(s)}{sQ(s)},$$(10)

where the denominator vanishes at s = 0 as a direct consequence of the definition of γ in Eq. (7), and Q(s) is a polynomial of degree n. The poles of $\hat{J}(s)$ are therefore s = 0 and the n roots of Q(s), which are real, negative, and simple (from the interlacing property in Eq. (22) of Gevorgyan et al. 2023).

The residue of $\hat{J}(s)$ at s = 0 comes entirely from the dashpot term in Eq. (8), the only singular contribution at the origin. Matching the residue computed from Eq. (10) gives $$\eta =\frac{Q(0)}{N(0)}.$$(11)

At each remaining pole q_i, only the i-th Kelvin–Voigt term is singular, and writing $1/({\alpha }_i+s{\eta }_i)={\eta }_i^{-1}/(s-q_i)$ with q_i = −α_i/η_i, residue matching yields $${\eta }_i=\frac{q_i{Q}^{\mathrm{\prime }}(q_i)}{N(q_i)},{\alpha }_i=-q_i{\eta }_i=\frac{-q_i^2{Q}^{\mathrm{\prime }}(q_i)}{N(q_i)}.$$(12)

No nonlinear system is solved. The only numerical step is the polynomial root-finding of Q(s), a standard operation for a polynomial whose roots are real, negative, and simple. For a layered body with n + 1 normal modes, the 2n + 3 data points $\{k_E^C,r_i^C,s_i\}$ uniquely determine the 2n + 3 parameters {γ, α, η, α_i, η_i} of the equivalent homogeneous model.

3 Software implementation and workflow

The package is written in Wolfram Language as a modular pipeline for constructing the homogeneous generalized Voigt rheology of a layered viscoelastic body. The model is not hardcoded to any specific planetary interior: all physical and numerical inputs are supplied by the user through a configuration file, and the same workflow applies to any incompressible, spherically symmetric multilayer model satisfying the assumptions discussed in Section 2.

3.1 Package structure

The code is organized into a small set of modules with distinct roles. The file ModelInput.wl contains the user-specified physical model and numerical settings. The core routines are called by computeHomogeneous, implemented in HomogeneousRheology.wl, which evaluates the layered tidal response, extracts the secular poles and residues, and returns the parameters of the equivalent generalized Voigt model. A driver script, RunHomogeneous.wl, loads the modules, executes the full pipeline, prints summary tables, and generates comparison plots. A separate script, VisualizeLayers.wl, produces a cross-sectional visualization of the layered body directly from the same input file.

3.2 Input configuration

The file ModelInput.wl contains both the physical description of the layered body and a small set of numerical controls for the inversion pipeline. The same workflow can be applied to any spherically symmetric multilayer configuration once the corresponding layer properties and forcing frequency are provided.

The physical model is specified through the arrays of outer radii, densities, rigidities, and viscosities of the layers, together with the tidal forcing frequency. An optional list of layer names may be supplied for visualization purposes. All physical quantities are given in SI units, and all arrays must have the same length, which determines the number of layers in the model. Liquid layers are identified by the condition 0 < μ < 1 Pa; for such layers, the method uses an effective elasticity value to regularize the small-shear limit (Sabadini et al. 2016).

The same input file contains numerical settings that control the precision of the calculations, the search interval used to detect the secular poles, and the frequency range and threshold adopted in the dominant-mode selection step. These parameters are not part of the planetary model but determine how the software resolves the relaxation spectrum and how aggressively the final generalized Voigt model is reduced. Their default values are listed in Table 1.

The single tidal forcing frequency ω and the frequency interval [ω_min, ω_max] play distinct roles. The rheological parameters of the equivalent homogeneous model are frequencyindependent, determined entirely by the secular poles and residues of the layered k₂(s). The single ω is used only to evaluate and report k₂(iω) at the frequency of interest for the body under consideration. For the Moon, this corresponds to the orbital frequency ω = 2.662 × 10⁻⁶ rad s⁻¹, associated with the 27.32-day orbital period, which we use throughout the case study in Section 4. The interval [ω_min, ω_max] enters only in the optional dominant-mode reduction step, where it defines the band over which each Kelvin–Voigt element is assessed for retention in the reduced model.

3.3 Main routine and computational workflow

The main routine

computeHomogeneous[rk,rhok,muk,etak,omega,opts]

performs the parameter extraction through the sequence of operations illustrated in Fig. 2. It follows the same three-step structure as in Sections 1 and 2: forward computation of the layered Love number, mode decomposition, and rheological inversion. The code first evaluates the elastic Love number k_E and the complex tidal response k₂(iω) for the layered body. It then searches for the zeros of the secular determinant on the negative real axis, refines them numerically, and computes the corresponding residues. From these quantities, it evaluates the fluid Love number k_f and the mapping constant Θ = 3I_°G/R⁵, and applies the algebraic procedure of Section 2.3 to extract the parameters {γ, α, η, α_i, η_i} of the equivalent generalized Voigt model.

The forward layered problem is evaluated through a unified 6 × 6 propagator matrix approach for all layers. Liquid layers are treated by small-shear regularization, such that solid and liquid regions are handled within the same numerical framework. The current package does not switch to a separate Clairaut propagator for liquids, but instead uses a single propagation formalism for arbitrary combinations of solid and liquid layers. A true liquid has μ = 0, but the Y-matrix (Sabadini et al. 2016, Eq. (2.42)) contains terms of the form μ/r in rows 3, 4, and 6. Setting μ = 0 exactly would therefore cause these rows to degenerate, causing the linear solver to fail when computing the propagator $$P(r_2,r_1)=Y(r_2)Y^{-1}(r_1).$$(13)

Instead, we assign a very small effective rigidity $\hat{\mu }\approx {10}^{-20}\mathrm{P}\mathrm{a}$ to the liquid layers. This keeps the matrix full-rank and invertible, while making the shear-stress rows effectively zero to machine precision. This value is negligible compared to the physical shear moduli, which are typically of the order of 10⁹−10¹¹ Pa, but remains large enough to preserve numerical stability at all reasonable (actually, even at quite stringent) working precisions. This is a standard regularization in viscoelastic normal-mode calculations that avoids introducing a separate 2 × 2 or 4 × 4 propagator for liquid layers while yielding results indistinguishable from the inviscid limit to within working precision.

After the full generalized Voigt model is obtained, the package performs an optional dominant-mode reduction, following the criterion introduced in Gevorgyan et al. (2023, Appendix A). For each Voigt element, it evaluates the maximum fractional contribution to – Im[C(iω)] over a user-defined frequency interval [ω_min, ω_max]. Elements whose contribution exceeds the prescribed threshold are kept in the reduced model. This produces a compact, effective rheology targeted to the frequency interval most relevant for the intended dynamical application.

Fig. 2 Workflow of the computeHomogeneous pipeline.

3.4 Output

The main routine computeHomogeneous returns the layered response and its homogeneous representation as a Wolfram Association, and also exports the corresponding results to a JSON file. The package also returns the intermediate modes used in the inversion. In this way, the user has access both to the physical outputs of interest – the effective generalized Voigt parameters – and to the pole–residue decomposition that connects the layered model to the homogeneous model. The quantities returned are summarized in Table 2.

The output falls into three groups: the Love numbers of the layered body, the constants used in the layered-to-homogeneous mapping, and the relaxation modes and rheological parameters of the equivalent generalized Voigt model. This organization separates diagnostic quantities from the final effective model. The Love-number and mode outputs are used to check the inversion, and the generalized Voigt parameters are used in downstream applications. The package also generates summary tables, dominant-mode diagnostics, and comparison plots. The user can inspect the mode decomposition, compare the full and reduced models, or export the reduced generalized Voigt parameters to an external integrator such as RheoVolution (de Oliveira et al. 2025a,b).

4 Case study: A five-layer model of the Moon

To test the package, we used the five-layer lunar model of Gevorgyan et al. (2023) to illustrate the equivalence between a multilayer Maxwell body and an effective homogeneous rheology. The model parameters are summarized in Table 3, and its cross-section is illustrated in Fig. 3; the figure was automatically produced by the VisualizeLayers.wl module. The Moon is the default example configuration in the package. The package was also tested for Io, Europa, and Enceladus using models and data from Hussmann & Spohn (2004), Wahr et al. (2006), Bierson & Nimmo (2016), and Thomas et al. (2016), although we are only going to report the results for the Moon. The ModelInput files for these additional test cases can be found in the software package repository.

Fig. 3 Cross-section depiction of the five-layer lunar interior model automatically generated by the VisualizeLayers.wl routine. The numbered list indicates the thickness of each corresponding layer, along with the position (radius) of its outer boundary in parentheses.

4.1 The homogeneous model

The software extracts the secular relaxation poles and residues listed in Table 4. These modes provide the normal-mode decomposition of the lunar tidal response and are used to construct the equivalent homogeneous rheology. The corresponding inversion yields the full equivalent generalized Voigt model, whose parameters are listed in Tables 5 and 6.

In Figure 4, we compare the dissipative part – Im[k₂(iω)] and the elastic part, Re[k₂(iω)], of the layered moon model to the full homogeneous model. The agreement shows that the full equivalent generalized Voigt model reproduces both the real and imaginary parts of the tidal Love number within numerical precision.

4.2 Dominant-mode reduction over a frequency interval

In many applications, it is unnecessary to retain the full effective rheology throughout the frequency range. The package includes a dominant-mode selection step that evaluates the maximum fractional contribution of each Voigt element to – Im[C(iω)] over a prescribed interval of forcing frequencies and retains only those elements whose contribution exceeds a given chosen threshold.

For the Moon, the relevant frequency interval is ω ∈ [10⁻¹⁰, 10⁻⁵] rad s⁻¹ (the default package setting), and we set the dominant-mode threshold to 2%. Under these assumptions, the full effective rheology is reduced to the subset of dominant Voigt elements marked with an asterisk in Table 5. As can be seen in Fig. 5, the reduced model preserves the main features of the layered model in the selected frequency interval.

The reduction step is important from a practical perspective. It allows the user to tailor the effective homogeneous model to the frequency range relevant for a specific dynamical problem rather than to retain all relaxation modes of the exact inversion. The reduced model is therefore specific to the chosen frequency interval and dominant-mode threshold.

5 Conclusions

We have presented a Wolfram Language package that constructs effective homogeneous generalized Voigt rheologies for incompressible, spherically symmetric layered viscoelastic bodies. The package computes the tidal Love number of a multilayer model, extracts its secular poles and residues, and converts this information into the parameters of an equivalent homogeneous rheology. Starting from a user-defined layered model, the package returns the parameters of the corresponding homogeneous generalized Voigt body in a format suitable for validation, comparison, and downstream numerical use and applications.

The lunar case study shows how the code converts a five-layer interior model into an equivalent homogeneous rheology. When all extracted Voigt elements are maintained, the homogeneous model reproduces the layered tidal response within numerical precision. The dominant-mode selection then gives a smaller model adapted to the prescribed frequency interval. The reduction applies when the relevant forcing frequencies lie in a limited range, since it preserves the tidal response of the layered model while reducing the number of rheological elements.

The package allows detailed interior-structure models to be translated into generalized Voigt parameters that can be supplied directly to time-domain tidal evolution codes, enabling their dissipation properties to be incorporated in long-term orbital evolution calculations without repeatedly solving the full layered deformation problem during integration. The same procedure can be applied to other incompressible multilayer models, including models of terrestrial planets and icy satellites, provided that the layer rheologies satisfy the assumptions used in the reduction.

The current implementation assumes Maxwell rheology for each solid layer. Maxwell bodies underestimate tidal dissipation at forcing frequencies close to or above the inverse Maxwell time, where laboratory and seismological evidence favors rheologies with a broader distribution of relaxation times, such as Andrade or Sundberg–Cooper. Our reduction pipeline is exact for a finite multilayer Maxwell body because the resulting k₂(s) is rational with a finite set of secular poles, a property that rheologies with a continuous relaxation spectrum do not share.

However, Andrade and similar rheologies can be approximated within a prescribed frequency range by a finite set of dominant modes, as done in Gevorgyan et al. (2020) for an extended Burgers model. Extending the package in this direction, together with tighter integration with external dynamical solvers, will require extending the current pole–residue reduction beyond the finite Maxwell case.

Fig. 4 Comparison between the exact five-layer lunar model and the equivalent homogeneous generalized Voigt model. The real and imaginary parts of the tidal Love number coincide within numerical precision when all extracted Voigt elements are retained.

Fig. 5 Same as Fig. 4, but for the reduced equivalent homogeneous generalized Voigt model.

Acknowledgements

The author would like to thank the anonymous referee for careful reading and constructive comments that improved the manuscript. The author would also like to thank Clodoaldo Ragazzo (USP) for many discussions and long-standing collaboration on this topic, Diogo Gomes and Alessandro Astolfi (KAUST) for their kind support, and J. Ricardo G. Mendonça (USP) for technical advice and a careful reading of the manuscript that greatly improved its presentation. This work was partially supported by the King Abdullah University of Science and Technology (KAUST), Saudi Arabia.

References

All Tables

All Figures

	Fig. 1 Schematic representation of the generalized Voigt model.
In the text

	Fig. 2 Workflow of the computeHomogeneous pipeline.
In the text

	Fig. 3 Cross-section depiction of the five-layer lunar interior model automatically generated by the VisualizeLayers.wl routine. The numbered list indicates the thickness of each corresponding layer, along with the position (radius) of its outer boundary in parentheses.
In the text

	Fig. 4 Comparison between the exact five-layer lunar model and the equivalent homogeneous generalized Voigt model. The real and imaginary parts of the tidal Love number coincide within numerical precision when all extracted Voigt elements are retained.
In the text

	Fig. 5 Same as Fig. 4, but for the reduced equivalent homogeneous generalized Voigt model.
In the text