A non-perturbative approach to computing seismic normal modes in rotating planets

TL;DR

This paper introduces a non-perturbative, high-performance computational method for accurately determining seismic normal modes in rotating planets, incorporating complex physics and validated against known models.

Contribution

It presents a novel finite-element based approach with advanced eigensolvers for seismic mode computation in rotating planets, including fluid and solid regions, validated against analytical and numerical benchmarks.

Findings

Validated methods against analytical solutions for non-rotating cases.

Successfully computed seismic modes in Earth and Mars models.

Analyzed effects of heterogeneity on seismic eigenfrequencies.

Abstract

A Continuous Galerkin method-based approach is presented to compute the seismic normal modes of rotating planets. Special care is taken to separate out the essential spectrum in the presence of a fluid outer core using a polynomial filtering eigensolver. The relevant elastic-gravitational system of equations, including the Coriolis force, is subjected to a mixed finite-element method, while self-gravitation is accounted for with the fast multipole method. Our discretization utilizes fully unstructured tetrahedral meshes for both solid and fluid regions. The relevant eigenvalue problem is solved by a combination of several highly parallel and computationally efficient methods. We validate our three-dimensional results in the non-rotating case using analytical results for constant elastic balls, as well as numerical results for an isotropic Earth model from standard ``radial" algorithms. We also validate the computations in the rotating case, but only in the slowly-rotating regime where perturbation theory applies, because no other independent algorithms are available in the general case. The algorithm and code are used to compute the point spectra of eigenfrequencies in several Earth and Mars models studying the effects of heterogeneity on a large range of scales.

Tables

Table 1. Table 1: Boundary conditions for a hydrostatic planet (cf. ( dahlen1998theoretical, , Table 3.4) ).

boundary types	linearized boundary conditions
free surface, $\partial \tilde{X}$	$T^0\cdot \nu =0;\nu \cdot T^{\mathrm{L1}}=0$
\hdashlinesolid-solid interfaces ${\Sigma }^{\text{SS}}$	${\left[\nu \cdot T^{\mathrm{L1}}\right]}_{-}^{+}=0;{\left[T^0\cdot \nu \right]}_{-}^{+}=0;{\left[u\right]}_{-}^{+}=0$
\hdashlinefluid-solid interfaces ${\Sigma }^{\text{FS}}$	${\left[T^0\cdot \nu \right]}_{-}^{+}=0;{\left[u\cdot \nu \right]}_{-}^{+}=0$
& fluid-fluid interfaces ${\Sigma }^{\text{FF}}$	${\left[\nu \cdot T^{\mathrm{L1}}\right]}_{-}^{+}=\nu {\left[\nu \cdot T^{\mathrm{L1}}\cdot \nu \right]}_{-}^{+}=0$
\hdashlineall interfaces $\Sigma$ & $\partial \tilde{X}$	${\left[S\left(u\right)\right]}_{-}^{+}=0;{\left[\nabla S\left(u\right)\cdot \nu +4\pi G{\rho }^{0}u\cdot \nu \right]}_{-}^{+}=0$

Table 2. Table 2: Bulk parameters of Earth and Mars; ϵ ˙ ​ ( r e ) ˙ italic-ϵ subscript 𝑟 𝑒 \dot{\epsilon}(r_{e}) denotes the derivative of ϵ italic-ϵ \epsilon at a = r e 𝑎 subscript 𝑟 𝑒 a=r_{e} , and ϵ ( r e ) hyd superscript subscript italic-ϵ subscript 𝑟 𝑒 hyd \epsilon_{(r_{e})}^{\text{hyd}} and ϵ ( r e ) obs superscript subscript italic-ϵ subscript 𝑟 𝑒 obs \epsilon_{(r_{e})}^{\text{obs}} denote the computed hydrostatic ellipticity and observed ellipticity, respectively.

parameters	$\Omega$ (s-1)	$r_e$ (km)	$g_{\left(r_e\right)}$ (m/s2)	$\dot{ϵ}\left(r_e\right)$	${ϵ}_{\left(r_e\right)}^{\text{hyd}}$	${ϵ}_{\left(r_e\right)}^{\text{obs}}$
Earth	7.2921$e-5$	6371.0	9.80	3.05$e-5>0$	3.34$e-3$	3.35$e-3$
Mars	7.0882$e-5$	3389.5	3.71	-8.98	N/A	5.89$e-3$

Table 3. Table 3: Implicit definition of the matrices. In the above, ∫ Ω S = ∑ k = 1 N K S ∫ K k S subscript superscript Ω S superscript subscript 𝑘 1 subscript superscript 𝑁 S 𝐾 subscript superscript subscript 𝐾 𝑘 S \int_{\Omega^{\text{S}}}=\sum_{k=1}^{N^{\text{S}}_{K}}\int_{K_{k}^{\text{S}}} , ∫ Ω F = ∑ k = 1 N K F ∫ K k F subscript superscript Ω F superscript subscript 𝑘 1 subscript superscript 𝑁 F 𝐾 subscript superscript subscript 𝐾 𝑘 F \int_{\Omega^{\text{F}}}=\sum_{k=1}^{N^{\text{F}}_{K}}\int_{K_{k}^{\text{F}}} and ∫ Σ FS = ∑ l = 1 N E FS ∫ E l FS subscript superscript Σ FS superscript subscript 𝑙 1 subscript superscript 𝑁 FS 𝐸 subscript superscript subscript 𝐸 𝑙 FS \int_{\Sigma^{\text{FS}}}=\sum_{l=1}^{N^{\text{FS}}_{E}}\int_{E_{l}^{\text{FS}}} .

operations	physical meanings	corresponding formulae
		${\int }_{{\Omega }^{\text{S}}}\nabla {\overline{v}}_h^s:(c:\nabla u_h^s)dx$
		$+{\int }_{{\Sigma }^{\text{FS}}}𝔖\left\{\left({\overline{v}}_h^s\cdot g^{\prime }\right)\left({\nu }^{s\to f}\cdot u_h^s\right){\left[{\rho }^0\right]}^f\right\}d\Sigma$
${\left({\tilde{v}}^s\right)}^{𝖧}{A}_{sg}{\tilde{u}}^s$	solid stiffness matrix with gravity	$+{\int }_{{\Omega }^{\text{S}}}𝔖\{(\nabla \cdot {\overline{v}}_h^s)(g^{\prime }\cdot u_h^s){\rho }^0$
		$-u_h^s\cdot (\nabla g^{\prime })\cdot {\overline{v}}_h^s{\rho }^0-u_h^s\cdot (\nabla {\overline{v}}_h^s)\cdot g^{\prime }{\rho }^0\}dx$
\hdashline		${\int }_{{\Omega }^{\text{F}}}{\rho }^0{N}^2\frac{\left(g^{\prime }\cdot {\overline{v}}_h^f\right)\left(g^{\prime }\cdot u_h^f\right)}{{\Vert g^{\prime }\Vert }^2}dx$
${\left({\tilde{v}}^f\right)}^{𝖧}{A}_f{\tilde{u}}^f$	buoyancy term	$+{\int }_{{\Sigma }^{\text{FF}}}\left(g^{\prime }\cdot \nu \right)\left({\overline{v}}_h^f\cdot \nu \right)\left(u_h^f\cdot \nu \right){\left[{\rho }^0\right]}_{-}^{+}d\Sigma$
\hdashline${\left({\tilde{v}}^p\right)}^{𝖧}{A}_p\tilde{p}$	fluid potential	$-{\int }_{{\Omega }^{\text{F}}}{\overline{v}}_h^p{p}_h{\kappa }^{-1}dx-{\int }_{{\Sigma }^{\text{FS}}}{\overline{v}}^p\left({\nu }^{f\to s}\cdot u^s\right)d\Sigma$
${\left({\tilde{v}}^f\right)}^{𝖧}{A}_{\text{dg}}\tilde{p}$	fluid stiffness matrix with gravity	${\int }_{{\Omega }^{\text{F}}}\left[{\overline{v}}_h^f\cdot \left(\nabla p_h\right)-\left({\overline{v}}_h^f\cdot g^{\prime }\right)p_h{\rho }^0{\kappa }^{-1}\right]dx$
${\left({\tilde{v}}^p\right)}^{𝖧}{A}_{\text{dg}}^{𝖳}{\tilde{u}}^f$	constraint with gravity	${\int }_{{\Omega }^{\text{F}}}\left[\left(\nabla {\overline{v}}_h^p\right)\cdot u_h^f-{\overline{v}}_h^p\left(g^{\prime }\cdot u_h^f\right){\rho }^0{\kappa }^{-1}\right]dx$
${\left({\tilde{v}}^s\right)}^{𝖧}{E}_{\text{FS}}\tilde{p}$	fluid-solid boundary condition	${\int }_{{\Sigma }^{\text{FS}}}\left({\overline{v}}_h^s\cdot {\nu }^{s\to f}\right)p_{h}d\Sigma$
${\left({\tilde{v}}^p\right)}^{𝖧}{E}_{\text{FS}}^{𝖳}{\tilde{u}}^f$	fluid-solid boundary condition	${\int }_{{\Sigma }^{\text{FS}}}-{\overline{v}}_h^p\left({\nu }^{f\to s}\cdot u_h^s\right)d\Sigma$
\hdashline${\left({\tilde{v}}^s\right)}^{𝖧}{R}_s{\tilde{u}}^s$	Coriolis force in ${\Omega }^{\text{S}}$	${\int }_{{\Omega }^{\text{S}}}{\overline{v}}_h^s\cdot \left(\Omega \times u_h^s\right){\rho }^{0}dx$
${\left({\tilde{v}}^f\right)}^{𝖧}{R}_f{\tilde{u}}^f$	Coriolis force in ${\Omega }^{\text{F}}$	${\int }_{{\Omega }^{\text{F}}}{\overline{v}}_h^f\cdot \left(\Omega \times u_h^f\right){\rho }^{0}dx$
\hdashline${\left({\tilde{v}}^s\right)}^{𝖧}{M}_s{\tilde{u}}^s$	solid mass matrix	${\int }_{{\Omega }^{\text{S}}}\left({\overline{v}}_h^s\cdot u_h^s\right){\rho }^{0}dx$
${\left({\tilde{v}}^f\right)}^{𝖧}{M}_f{\tilde{u}}^f$	fluid mass matrix	${\int }_{{\Omega }^{\text{F}}}\left({\overline{v}}_h^f\cdot u_h^f\right){\rho }^{0}dx$

Table 4. Table 4: Implicit definition of the submatrices for perturbation to the gravitational potential.

operations	physical meanings	corresponding formulae
		${\int }_{{\Omega }^{\text{S}}}\nabla \cdot \left({\rho }^0{u}_h^s\right)dx,$
		${\int }_{{\Sigma }^{\text{FS}}}\left({\nu }^{f\to s}\cdot u_h^s\right){\left[{\rho }^0\right]}^{s}dx,$
$C_s{\tilde{u}}^s$	N bodies in $\overline{{\Omega }^{\text{S}}}$	${\int }_{{\Sigma }^{\text{SS}\cup \partial {\tilde{X}}^{\text{S}}}}\left(\nu \cdot u_h^s\right){\left[{\rho }^0\right]}_{-}^{+}dx$
\hdashline		${\int }_{{\Omega }^{\text{F}}}\nabla \cdot \left({\rho }^0{u}_h^f\right)dx,$
		${\int }_{{\Sigma }^{\text{FS}}}\left({\nu }^{s\to f}\cdot u_h^f\right){\left[{\rho }^0\right]}^{f}dx,$
$C_f{\tilde{u}}^f$	N bodies in $\overline{{\Omega }^{\text{F}}}$	${\int }_{{\Sigma }^{\text{FF}\cup \partial {\tilde{X}}^{\text{F}}}}\left(\nu \cdot u_h^f\right){\left[{\rho }^0\right]}_{-}^{+}dx$
\hdashline		$G{\int }_{\tilde{X}}\frac{{\nabla }^{\prime }\cdot \left({\rho }^0\left(x^{\prime }\right)u_h\left(x^{\prime }\right)\right)}{\Vert x-x^{\prime }\Vert }d{x}^{\prime }$
$S\left(C\tilde{u}\right)$	solution for Poisson’s equation	$+G{\int }_{\Sigma \cup \partial \tilde{X}}\frac{\nu \left(x^{\prime }\right)\cdot u_h\left(x^{\prime }\right){\left[{\rho }^0\left(x^{\prime }\right)\right]}_{-}^{+}}{\Vert x-x^{\prime }\Vert }d{x}^{\prime }$
\hdashline		${\int }_{{\Omega }^{\text{S}}}\nabla \cdot \left({\rho }^0{\overline{v}}_h^s\right)S\left(u_h\right)dx$
	incremental gravitational field	$+{\int }_{{\Sigma }^{\text{FS}}}\left({\overline{v}}_h^s\cdot {\nu }^{f\to s}\right)S\left(u_h\right){\left[{\rho }^0\right]}^{s}dx$
${\left({\tilde{v}}^s\right)}^{𝖧}{C}_s^{𝖳}\left(SC\tilde{u}\right)$	in $\overline{{\Omega }^{\text{S}}}$	$+{\int }_{{\Sigma }^{\text{SS}}\cup \partial {\tilde{X}}^{\text{S}}}\left({\overline{v}}_h^s\cdot \nu \right)S\left(u_h\right){\left[{\rho }^0\right]}_{-}^{+}dx$
\hdashline		${\int }_{{\Omega }^{\text{F}}}\nabla \cdot \left({\rho }^0{\overline{v}}_h^f\right)S\left(u_h\right)dx$
	incremental gravitational field	$+{\int }_{{\Sigma }^{\text{FS}}}\left({\overline{v}}_h^f\cdot {\nu }^{s\to f}\right)S\left(u_h\right){\left[{\rho }^0\right]}^{f}dx$
${\left({\tilde{v}}^s\right)}^{𝖧}{C}_f^{𝖳}\left(SC\tilde{u}\right)$	in $\overline{{\Omega }^{\text{F}}}$	$+{\int }_{{\Sigma }^{\text{FF}}\cup \partial {\tilde{X}}^{\text{F}}}\left({\overline{v}}_h^f\cdot \nu \right)S\left(u_h\right){\left[{\rho }^0\right]}_{-}^{+}dx$

Table 5. Table 5: Errors in the gravitational calculation of a constant density ball.

# of elements	116,085	1,136,447	2,019,017	3,081,551	4,035,022
MSE of ${\Phi }^0$	2.133e-6	7.452e-8	1.784e-8	1.545e-8	1.430e-8
MSE of $g$	1.102e-3	1.848e-4	1.156e-4	8.781e-5	7.363e-5

Table 6. Table 6: Errors of three-layer approximations in the gravitational calculation.

# of elements	5,800	57,490	503,882	1,136,447	2,093,055	5,549,390	7,825,918
MSE of ${\Phi }^0$	3.604e-3	2.635e-4	4.071e-5	2.092e-5	1.354e-5	4.059e-6	2.396e-9
MSE of $g$	5.805e-2	5.479e-3	7.320e-4	3.218e-4	2.068e-4	9.524e-5	5.609e-5

Table 7. Table 7: Errors of seven-layer approximations in the gravitational calculation.

# of elements	2,031,729	5,018,249	8,043,617	12,479,828	16,560,615
MSE of ${\Phi }^0$	2.333e-7	4.485e-8	1.286e-8	9.785e-9	5.548e-9
MSE of $g$	1.926e-4	8.606e-5	5.186e-5	4.036e-6	3.394e-5

Table 8. Table 8: Test cases with self-gravitation for different solid models using P1 elements for the frequency range [ 0.1 , 1.0 ] 0.1 1.0 [0.1,1.0] mHz.

Exp.	#elm.	size of $A$	#surf.	size of $S_s$	$[f_1,f_2]$ (mHz)	(deg,#it)	#eigs
C1p1	5,123	2,727	392	5,515	[0.1,1.0]	(14,192)	70
C2p1	21,093	10,644	956	22,049	[0.1,1.0]	(25,232)	92
C3p1	39,273	19,131	956	40,229	[0.1,1.0]	(34,252)	92
C4p1	105,115	51,933	3,608	108,723	[0.1,1.0]	(50,252)	92
C5p1	495,099	242,721	14,888	509,987	[0.1,1.0]	(108,272)	92

Table 9. Table 9: Convergence tests with self-gravitation for different solid models in Table  8 with self-gravitation for P1 elements.

Exp.	${}_0{}^{}T_2^{}$	${}_0{}^{}S_2^{}$	${}_1{}^{}S_1^{}$	${}_0{}^{}S_0^{}$	${}_0{}^{}T_3^{}$	${}_0{}^{}S_3^{}$	${}_1{}^{}S_2^{}$	${}_0{}^{}T_4^{}$	${}_0{}^{}S_4^{}$
C1p1	0.3724	0.4178	0.4600	0.5105	0.5881	0.6322	0.6900	0.7973	0.8359
C2p1	0.3653	0.4112	0.4511	0.5053	0.5692	0.6052	0.6708	0.7587	0.7791
C3p1	0.3643	0.4103	0.4502	0.5053	0.5665	0.6017	0.6680	0.7527	0.7721
C4p1	0.3622	0.4089	0.4472	0.5035	0.5612	0.5932	0.6622	0.7424	0.7526
C5p1	0.3612	0.4086	0.4460	0.5035	0.5587	0.5899	0.6596	0.7374	0.7445
1D	0.3607	0.4087	0.4456	0.5040	0.5574	0.5885	0.6582	0.7348	0.7406

Table 10. Table 10: Test cases with self-gravitation for different solid models using P2 elements for the frequency range [ 0.1 , 1.0 ] 0.1 1.0 [0.1,1.0] mHz.

Exp.	# of elm.	size of $A$	#surf.	size of $S_s$	$[f_1,f_2]$ (mHz)	(deg,#it)	#eigs
C1p2	19,073	75,888	956	20,029	[0.1,1.0]	(44,512)	92
C2p2	40,378	170,025	3,608	43,986	[0.1,1.0]	(58,492)	92
C3p2	80,554	335,103	5,924	86,478	[0.1,1.0]	(81,492)	92
C4p2	152,426	645,687	14,888	167,314	[0.1,1.0]	(129,492)	92
C5p2	334,193	1,360,140	14,888	349,081	[0.1,1.0]	(200,492)	92

Table 11. Table 11: Convergence tests with self-gravitation for the solid models in Table  10 using P2 elements.

Exp.	${}_0{}^{}T_2^{}$	${}_0{}^{}S_2^{}$	${}_1{}^{}S_1^{}$	${}_0{}^{}S_0^{}$	${}_0{}^{}T_3^{}$	${}_0{}^{}S_3^{}$	${}_1{}^{}S_2^{}$	${}_0{}^{}T_4^{}$	${}_0{}^{}S_4^{}$
C1p2	0.3619	0.4100	0.4473	0.5094	0.5594	0.5908	0.6605	0.7376	0.7439
C2p2	0.3610	0.4090	0.4459	0.5042	0.5579	0.5889	0.6587	0.7355	0.7413
C3p2	0.3609	0.4089	0.4463	0.5042	0.5577	0.5888	0.6585	0.7352	0.7410
C4p2	0.3608	0.4088	0.4456	0.5041	0.5575	0.5886	0.6583	0.7349	0.7408
C5p2	0.3608	0.4087	0.4456	0.5041	0.5575	0.5885	0.6583	0.7349	0.7407
1D	0.3607	0.4087	0.4456	0.5040	0.5574	0.5885	0.6582	0.7348	0.7406

Table 12. Table 12: Test cases with self-gravitation for different Earth models E1p1 - E8p1 using P1 elements for the frequency range [ 0.1 , 1.0 ] 0.1 1.0 [0.1,1.0] mHz and Earth model E8p2 using P2 elements for the frequency range [ 0.3 , 1.5 ] 0.3 1.5 [0.3,1.5] mHz.

Exp.	# of elm.	size of $A_G$	size of $A_p$	#surf.	size of $S$	$[f_1,f_2]$ (mHz)	(deg,#it)	#eigs
E1p1	9,721	7,590	887	2,304	12,025	[0.1,1.0]	(187,392)	64
E2p1	20,466	14,736	974	4,956	25,422	[0.1,1.0]	(182,372)	72
E3p1	42,828	30,384	3,171	8,172	51,000	[0.1,1.0]	(342,452)	83
E4p1	83,354	63,225	5,298	22,104	105,458	[0.1,1.0]	(745,452)	88
E5p1	157,057	96,852	6,771	22,104	179,161	[0.1,1.0]	(747,492)	88
E6p1	303,218	164,673	10,077	22,104	325,322	[0.1,1.0]	(685,492)	88
E7p1	639,791	361,587	21,824	60,288	700,079	[0.1,1.0]	(685,492)	88
E8p1	1,972,263	1,086,702	70,429	150,288	2,122,551	[0.1,1.0]	(1565,492)	88
E8p2	1,972,263	8,400,630	522,705	150,288	2,122,551	[0.3,1.5]	(1185,1051)	268

Table 13. Table 13: Convergence tests with self-gravitation for different Earth models in Table  12 .

Exp.	${}_0{}^{}S_2^{}$	${}_0{}^{}T_2^{}$	${}_2{}^{}S_1^{}$	${}_0{}^{}S_3^{}$	${}_0{}^{}T_3^{}$
E1p1	0.3284	0.3953	0.4179	0.5242	0.6241
E2p1	0.3229	0.3921	0.4149	0.5077	0.6146
E3p1	0.3177	0.3884	0.4113	0.4932	0.6062
E4p1	0.3166	0.3842	0.4090	0.4903	0.5980
E5p1	0.3137	0.3845	0.4085	0.4863	0.5962
E6p1	0.3126	0.3840	0.4080	0.4768	0.5945
E7p1	0.3116	0.3834	0.4073	0.4742	0.5933
E8p1	0.3112	0.3829	0.4067	0.4721	0.5920
E8p2	0.3106	0.3826	0.4063	0.4708	0.5912
1D	0.3110	0.3826	0.4063	0.4713	0.5912

Table 14. Table 14: Test case with self-gravitation for an Earth model with a non-spherically symmetric CMB using P2 elements.

Exp.	# of elm.	size of $A_G$	size of $A_p$	#surf.	size of $S$	$[f_1,f_2]$ (mHz)	(deg,#it)	#eigs
CMB8	2,007,479	8,711,940	633,358	177,352	2,184,831	[1.5,2.0]	(3591,1251)	350

Table 15. Table 15: Test cases for four different Earth models using the Cowling approximation.

Exp.	# of elm.	size of $A_G$	size of $A_p$	$[f_1,f_2]$ (mHz)	(deg,#it)	#eigs
E9p2	4,094,031	17,469,666	1,181,103	$[2.0,2.5]$	(4054,1892)	528
MIT_2016May	4,048,932	16,578,945	879,067	$[2.0,2.5]$	(2674,1912)	520
MIT+crust 1.0	4,044,225	16,550,922	878,808	$[2.0,2.5]$	(6984,1912)	550

Table 16. Table 16: Numerical parameter values pertaining to the testing of computational accuracy and estimating the cost in different models.

Exp.	# of elm.	size of $A$	size of $A_p$	size of $S$	$[f_1,f_2]$ (mHz)
Constant (C3kp1)	3,129	1,821	0	3,521	[0.35,0.85]
Earth (E3kp1)	3,330	2,760	392	4,242	[0.3,0.86]
Mars (M2kp1)	1,887	1,677	145	2,539	[0.4,1.14]
Mars (M8kp1)	8,020	7,557	152	12,436	[0.4, 1.14]
Earth (E40kp1)	42,828	30,384	3,171	51,000	[0.1,1.5]

Table 17. Table 17: Numerical parameters of the Earth models used in the benchmark experiments.

Exp.	# of elm.	size of $A$	size of $A_p$	size of $S$	$[f_1,f_2]$ (mHz)
Earth (E1Mp1)	1,011,973	537,198	31,849	1,074,577	[0.04,1.5]
Earth (E2Mp2)	2,015,072	8,569,197	530,721	2,165,360	[0.2,1.5]
Earth (EE1Mp1)	1,003,065	533,064	31,688	1,065,629	[0.04,1.5]
Earth (EE2Mp2)	2,002,581	8,520,432	528,124	2,153,109	[0.2,1.5]

Table 18. Table 18: Numerical parameters for the Mars models.

Exp.	# of elm.	size of $A$	size of $A_p$	size of $S$	$[f_1,f_2]$ (mHz)
Mars (M2Mp2)	1,996,773	8,967,684	579,338	2,257,801	[0.2,2.0]
Mars (EM2Mp2)	2,001,619	8,984,532	579,667	2,262,143	[0.2,2.0]
Mars (TM2Mp2)	2,008,654	8,289,927	323,810	2,158,366	[0.2,2.0]

Table 19. Table 19: Implicit definition of the matrices in ( 87 ) (no summations over k 𝑘 k and m 𝑚 m ). Since the construction of A s ​ g ( r ) superscript subscript 𝐴 𝑠 𝑔 𝑟 A_{sg}^{(r)} is standard, we refer to ( dahlen1998theoretical, , (8.43) & (8.44)) and ( jingchen2018revisiting, , (3.1)) . In the above, ∫ Ω ( r ) S = ∑ q = 1 N L S ∫ L q S subscript superscript subscript Ω 𝑟 S superscript subscript 𝑞 1 superscript subscript 𝑁 𝐿 S subscript superscript subscript 𝐿 𝑞 S \int_{\Omega_{(r)}^{\text{S}}}=\sum_{q=1}^{N_{L}^{\text{S}}}\int_{L_{q}^{\text{S}}} and ∫ Ω ( r ) F = ∑ q = 1 N L F ∫ L q F subscript superscript subscript Ω 𝑟 F superscript subscript 𝑞 1 superscript subscript 𝑁 𝐿 F subscript superscript subscript 𝐿 𝑞 F \int_{\Omega_{(r)}^{\text{F}}}=\sum_{q=1}^{N_{L}^{\text{F}}}\int_{L_{q}^{\text{F}}} .

operations	physical meanings	corresponding formulae
${\left({\tilde{U}}_{km}^{s^{\prime }}\right)}^{𝖳}{A}_{sg}^{\left(r\right)}{\tilde{U}}_{km}^s$	solid stiffness matrix	(jingchen2018revisiting, , (3.1))
\hdashline		${\int }_{{\Omega }_{\left(r\right)}^{\text{F}}}{U}_{km;h}^{f^{\prime }}{U}_{km;h}^f{N}_{\left(r\right)}^2{\rho }_{\left(r\right)}^0{r}^{2}dr$
${\left({\tilde{U}}_{km}^{f^{\prime }}\right)}^{𝖳}{A}_f^{\left(r\right)}{\tilde{U}}_{km}^f$	buoyancy term	$+{{\left[{\rho }_{\left(r\right)}^0\right]}_{-}^{+}{g}_{\left(r\right)}{U}_{km;h}^f{U}_{km;h}^{f^{\prime }}{r}^2|}_{{\Sigma }_{\left(r\right)}^{\text{FF}}}$
\hdashline${\left({\tilde{P}}_{km}^{\prime }\right)}^{𝖳}{A}_p^{\left(r\right)}{\tilde{P}}_{km}$	fluid potential	${\int }_{{\Omega }_{\left(r\right)}^{\text{F}}}{P}_{km;h}^{\prime }{P}_{km;h}{\kappa }_{\left(r\right)}^{-1}{r}^{2}dr$
\hdashline		${\int }_{{\Omega }_{\left(r\right)}^{\text{F}}}{U}_{km;h}^{f^{\prime }}\left({\partial }_r{P}_{km;h}+{\rho }_{\left(r\right)}^0{g}_{\left(r\right)}{\kappa }_{\left(r\right)}^{-1}{P}_{km;h}\right)r^{2}dr$
${\left({\tilde{u}}_{km}^{f^{\prime }}\right)}^{𝖳}{A}_{\text{dg}}^{\left(r\right)}{\tilde{P}}_{km}$	fluid stiffness matrix	$+{\int }_{{\Omega }_{\left(r\right)}^{\text{F}}}\sqrt{l\left(l+1\right)}{P}_{km;h}{V}_{km;h}^{f^{\prime }}rdr$
\hdashline		${\int }_{{\Omega }_{\left(r\right)}^{\text{F}}}\left({\partial }_r{P}_{km;h}^{\prime }+{\rho }_{\left(r\right)}^0{g}_{\left(r\right)}{\kappa }_{\left(r\right)}^{-1}{P}_{km;h}^{\prime }\right)U_{km;h}^f{r}^{2}dr$
${\left({\tilde{P}}_{km}^{\prime }\right)}^{𝖳}A_{\text{dg}}^{\left(r\right)}{}_{}{}^{𝖳}{\tilde{u}}_{km}^f$	constraint	$+{\int }_{{\Omega }_{\left(r\right)}^{\text{F}}}\sqrt{l\left(l+1\right)}{P}_{km;h}^{\prime }{V}_{km;h}^{f}rdr$
\hdashline${\left({\tilde{U}}_{km}^{s^{\prime }}\right)}^{𝖳}{E}_{\text{FS}}^{\left(r\right)}{\tilde{P}}_{km}$	fluid-solid boundary condition	$-{P_{km;h}{U}_{km;h}^{s^{\prime }}{r}^2|}_{{\Sigma }_{\left(r\right)}^{\text{FS}}}$
${\left({\tilde{P}}_{km}^{\prime }\right)}^{𝖳}E_{\text{FS}}^{\left(r\right)}{}_{}{}^{𝖳}{\tilde{U}}_{km}^s$	fluid-solid boundary condition	$-{P_{km;h}^{\prime }{U}_{km;h}^s{r}^2|}_{{\Sigma }_{\left(r\right)}^{\text{FS}}}$
${\left({\tilde{U}}_{km}^{s^{\prime }}\right)}^{𝖳}{M}_s^{\left(r\right)}{\tilde{U}}_{km}^s$	solid mass matrix	${\int }_{{\Omega }_{\left(r\right)}^{\text{S}}}\left(U_{km;h}^{s^{\prime }}{U}_{km;h}^s+V_{km;h}^{s^{\prime }}{V}_{km;h}^s\right){\rho }_{\left(r\right)}^0{r}^{2}dr$
${\left({\tilde{U}}_{km}^{f^{\prime }}\right)}^{𝖳}{M}_f^{\left(r\right)}{\tilde{U}}_{km}^f$	fluid mass matrix	${\int }_{{\Omega }_{\left(r\right)}^{\text{F}}}\left(U_{km;h}^{f^{\prime }}{U}_{km;h}^f+V_{km;h}^{f^{\prime }}{V}_{km;h}^f\right){\rho }_{\left(r\right)}^0{r}^{2}dr$

Table 20. Table 20: Implicit definition of the matrices in ( 87 ) (no summation over k 𝑘 k and m 𝑚 m ). In the above, ∫ Ω ( r ) S = ∑ q = 1 N L S ∫ L q S subscript superscript subscript Ω 𝑟 S superscript subscript 𝑞 1 superscript subscript 𝑁 𝐿 S subscript superscript subscript 𝐿 𝑞 S \int_{\Omega_{(r)}^{\text{S}}}=\sum_{q=1}^{N_{L}^{\text{S}}}\int_{L_{q}^{\text{S}}} and ∫ Ω ( r ) F = ∑ q = 1 N L F ∫ L q F subscript superscript subscript Ω 𝑟 F superscript subscript 𝑞 1 superscript subscript 𝑁 𝐿 F subscript superscript subscript 𝐿 𝑞 F \int_{\Omega_{(r)}^{\text{F}}}=\sum_{q=1}^{N_{L}^{\text{F}}}\int_{L_{q}^{\text{F}}} . In the Poisson’s equation, the computation of the integral ∫ 0 ∞ superscript subscript 0 \int_{0}^{\infty} requires special treatment, see ( jingchen2018revisiting, , Chapter 3.2.2) .

operations	physical meanings	corresponding formulae
		${\int }_{{\Omega }_{\left(r\right)}^{\text{S}}}\left({\partial }_r{S}_{km;h}^{\prime }\right)U_{km;h}^s{\rho }_{\left(r\right)}^0{r}^{2}dr$
${\left({\tilde{S}}_{km}^{\prime }\right)}^{𝖳}{C}_s^{\left(r\right)}{\tilde{u}}_{km}^s$	density changes in $\overline{{\Omega }_{\left(r\right)}^{\text{S}}}$	$+{\int }_{{\Omega }_{\left(r\right)}^{\text{S}}}\sqrt{l\left(l+1\right)}{S}_{km;h}^{\prime }{V}_{km;h}^s{\rho }_{\left(r\right)}^{0}rdr$
\hdashline		${\int }_{{\Omega }_{\left(r\right)}^{\text{F}}}\left({\partial }_r{S}_{km;h}^{\prime }\right)U_{km;h}^f{\rho }_{\left(r\right)}^0{r}^{2}dr$
${\left({\tilde{S}}_{km}^{\prime }\right)}^{𝖳}{C}_f^{\left(r\right)}{\tilde{u}}_{km}^f$	density changes in $\overline{{\Omega }_{\left(r\right)}^{\text{F}}}$	$+{\int }_{{\Omega }_{\left(r\right)}^{\text{F}}}\sqrt{l\left(l+1\right)}{S}_{km;h}^{\prime }{V}_{km;h}^f{\rho }_{\left(r\right)}^{0}rdr$
\hdashline		${\left(4\pi G\right)}^{-1}{\int }_0^{\infty }({\partial }_r{S}_{km;h}^{\prime }{\partial }_r{S}_{km;h}{r}^2$
${\left({\tilde{S}}_{km}^{\prime }\right)}^{𝖳}{S}^{\left(r\right)}{\tilde{S}}_{km}$	Poisson’s equation	$+l(l+1)S_{km;h}^{\prime }{S}_{km;h})dr$
\hdashline	incremental gravitational field	${\int }_{{\Omega }_{\left(r\right)}^{\text{S}}}{U}_{km;h}^{s^{\prime }}\left({\partial }_r{S}_{km;h}\right){\rho }_{\left(r\right)}^0{r}^{2}dr$
${\left({\tilde{u}}_{km}^{s^{\prime }}\right)}^{𝖳}C_s^{\left(r\right)}{}_{}{}^{𝖳}{\tilde{S}}_{km}$	in $\overline{{\Omega }_{\left(r\right)}^{\text{S}}}$	$+{\int }_{{\Omega }_{\left(r\right)}^{\text{S}}}\sqrt{l\left(l+1\right)}{V}_{km;h}^{s^{\prime }}{S}_{km;h}{\rho }_{\left(r\right)}^{0}rdr$
\hdashline	incremental gravitational field	${\int }_{{\Omega }_{\left(r\right)}^{\text{F}}}{U}_{km;h}^{f^{\prime }}\left({\partial }_r{S}_{km;h}\right){\rho }_{\left(r\right)}^0{r}^{2}dr$
${\left({\tilde{u}}_{km}^{f^{\prime }}\right)}^{𝖳}C_f^{\left(r\right)}{}_{}{}^{𝖳}{\tilde{S}}_{km}$	in $\overline{{\Omega }_{\left(r\right)}^{\text{F}}}$	$+{\int }_{{\Omega }_{\left(r\right)}^{\text{F}}}\sqrt{l\left(l+1\right)}{V}_{km;h}^{f^{\prime }}{S}_{km;h}{\rho }_{\left(r\right)}^{0}rdr$

Table 21. Table 21: Implicit definition of the matrices in ( 88 ) for the Cowling approximation.

operations	physical meanings	corresponding formulae
		$\sum _{km}\sum _{k^{\prime }{m}^{\prime }}{y}_{k^{\prime }{m}^{\prime }}^{\prime }\{{\int }_{{\Omega }^{\text{S}}}\nabla u_{k^{\prime }{m}^{\prime };h}:(c:\nabla u_{km;h})dx$
		$+{\int }_{{\Sigma }^{\text{FS}}}𝔖\left\{\left(g\cdot u_{k^{\prime }{m}^{\prime };h}\right)\left({\nu }^{s\to f}\cdot u_{km;h}\right){\left[{\rho }^0\right]}^f\right\}d\Sigma$
		$+{\int }_{{\Omega }^{\text{S}}}𝔖\{(\nabla \cdot u_{k^{\prime }{m}^{\prime };h})({\rho }^0{u}_{km;h}\cdot g)-{\rho }^0{u}_{k^{\prime }{m}^{\prime };h}\cdot (\nabla g)\cdot u_{km;h}$
		$-{\rho }^0{u}_{km;h}\cdot (\nabla u_{k^{\prime }{m}^{\prime };h})\cdot g\}dx$
		$+{\int }_{{\Omega }^{\text{F}}}{\rho }^0{N}^2\frac{\left(g\cdot u_{k^{\prime }{m}^{\prime };h}\right)\left(g\cdot u_{km;h}\right)}{{\Vert g\Vert }^2}dx$
		$+{\int }_{{\Sigma }^{\text{FF}}}\left(g\cdot \nu \right)\left(u_{km;h}\cdot \nu \right)\left(u_{k^{\prime }{m}^{\prime };h}\cdot \nu \right){\left[{\rho }^0\right]}_{-}^{+}d\Sigma$
		$+{\int }_{{\Omega }^{\text{F}}}\kappa \left(\nabla u_{k^{\prime }{m}^{\prime };h}+{\rho }^0{\kappa }^{-1}{u}_{k^{\prime }{m}^{\prime };h}\cdot g\right)$
${\left(y^{\prime }\right)}^{𝖳}{A}_G^{\left(c\right)}y$	stiffness matrix	$(\nabla u_{km;h}+{\rho }^0{\kappa }^{-1}{u}_{km;h}\cdot g)dx\}y{}_{km}$
\hdashline		$\sum _{km}\sum _{k^{\prime }{m}^{\prime }}{y}_{k^{\prime }{m}^{\prime }}^{\prime }\{{\int }_{{\Omega }^{\text{S}}}{u}_{k^{\prime }{m}^{\prime };h}\cdot u_{km;h}{\rho }^{0}dx$
${\left(y^{\prime }\right)}^{𝖳}{M}^{\left(c\right)}y$	mass matrix	$+{\int }_{{\Omega }^{\text{F}}}{u}_{k^{\prime }{m}^{\prime };h}\cdot u_{km;h}{\rho }^{0}dx\}y{}_{km}$

Table 22. Table 22: Implicit definition of the matrices in ( 88 ).

operations	physical meanings	corresponding formulae
		$\sum _{km}\sum _{k^{\prime }{m}^{\prime }}{z}_{k^{\prime }{m}^{\prime }}^{\prime }\{{\int }_{{\Omega }^{\text{S}}}{s}_{k^{\prime }{m}^{\prime };h}\nabla \cdot \left({\rho }^0{u}_{km;h}\right)dx$
		$+{\int }_{{\Sigma }^{\text{SS}\cup \partial {\tilde{X}}^{\text{S}}}}{s}_{k^{\prime }{m}^{\prime };h}\nu \cdot u_{km;h}{\left[{\rho }^0\right]}_{-}^{+}d\Sigma$
		$+{\int }_{{\Sigma }^{\text{FS}}}{s}_{k^{\prime }{m}^{\prime };h}{\nu }^{f\to s}\cdot u_{km;h}{\left[{\rho }^0\right]}^{s}d\Sigma$
		$+{\int }_{{\Sigma }^{\text{FF}\cup \partial {\tilde{X}}^{\text{F}}}}{s}_{k^{\prime }{m}^{\prime };h}\nu \cdot u_{km;h}{\left[{\rho }^0\right]}_{-}^{+}d\Sigma$
		$+{\int }_{{\Omega }^{\text{F}}}{s}_{k^{\prime }{m}^{\prime };h}\nabla \cdot \left({\rho }^0{u}_{km;h}\right)dx$
${\left(z^{\prime }\right)}^{𝖳}{C}^{\left(c\right)}y$	density changes in $\tilde{X}$	$+{\int }_{{\Sigma }^{\text{FS}}}{s}_{k^{\prime }{m}^{\prime };h}{\nu }^{s\to f}\cdot u_{km;h}{\left[{\rho }^0\right]}^{f}d\Sigma \}y{}_{km}$
\hdashline${\left(z^{\prime }\right)}^{𝖳}{S}^{\left(c\right)}z$	Poisson’s equation	$\sum _{km}\sum _{k^{\prime }{m}^{\prime }}{z}_{k^{\prime }{m}^{\prime }}^{\prime }\left\{{\int }_{{ℝ}^3}\left(\nabla s_{k^{\prime }{m}^{\prime };h}\right)\cdot \left(\nabla s_{km;h}\right)dx\right\}{z}_{km}$
\hdashline		$\sum _{km}\sum _{k^{\prime }{m}^{\prime }}{y}_{k^{\prime }{m}^{\prime }}^{\prime }\{{\int }_{{\Omega }^{\text{S}}}\nabla \cdot \left({\rho }^0{u}_{k^{\prime }{m}^{\prime };h}\right)s_{km;h}dx$
		$+{\int }_{{\Sigma }^{\text{SS}}\cup \partial {\tilde{X}}^{\text{S}}}{\left[{\rho }^0\right]}_{-}^{+}\nu \cdot u_{k^{\prime }{m}^{\prime };h}{s}_{km;h}d\Sigma$
		$+{\int }_{{\Sigma }^{\text{FS}}}{\left[{\rho }^0\right]}^s{\nu }^{f\to s}\cdot u_{k^{\prime }{m}^{\prime };h}{s}_{km;h}d\Sigma$
		$+{\int }_{{\Sigma }^{\text{FF}}\cup \partial {\tilde{X}}^{\text{F}}}{\left[{\rho }^0\right]}_{-}^{+}\nu \cdot u_{k^{\prime }{m}^{\prime };h}{s}_{km;h}d\Sigma$
	incremental gravitational field	$+{\int }_{{\Omega }^{\text{F}}}\nabla \cdot \left({\rho }^0{u}_{k^{\prime }{m}^{\prime };h}\right)s_{km;h}dx$
${\left(y^{\prime }\right)}^{𝖳}C_{}^{\left(c\right)}{}_{}{}^{𝖳}z$	in $\tilde{X}$	$+{\int }_{{\Sigma }^{\text{FS}}}{\left[{\rho }^0\right]}^f{\nu }^{s\to f}\cdot u_{k^{\prime }{m}^{\prime };h}{s}_{km;h}d\Sigma \}z{}_{km}$

Equations

$$\tilde{X}={\Omega }^{\text{S}}\cup {\Omega }^{\text{F}}\cup \Sigma \cup \partial \tilde{X},\partial \tilde{X}=\partial {\tilde{X}}^{\text{S}}\cup \partial {\tilde{X}}^{\text{F}},{\Omega }^{\text{S}}={\Omega }^{\text{CM}}\cup {\Omega }^{\text{IC}},{\Omega }^{\text{F}}={\Omega }^{\text{OC}}\cup {\Omega }^{\text{O}}.$$

$$\Sigma ={\Sigma }^{\text{SS}}\cup {\Sigma }^{\text{FF}}\cup {\Sigma }^{\text{FS}},{\Sigma }^{\text{FS}}={\Sigma }_{\text{int}}^{\text{FS}}\cup {\Sigma }_{\text{O}}^{\text{FS}},{\Sigma }_{\text{int}}^{\text{FS}}={\Sigma }^{\text{CMB}}\cup {\Sigma }^{\text{ICB}}.$$

$$\Delta {\Phi }^0=4\pi G{\rho }^0,$$

$$\Delta S(u)=-4\pi G\nabla \cdot ({\rho }^{0}u).$$

$$\psi (x)=-\frac{1}{2}\left[{\Omega }^2{x}^2-{\left(\Omega \cdot x\right)}^2\right],$$

$$g^{\prime }=g-\nabla \psi =-\nabla ({\Phi }^0+\psi ),$$

$$g=-\nabla {\Phi }^0.$$

$$\nabla \cdot T^0=-{\rho }^0{g}^{\prime }.$$

$$-{\omega }^2{\rho }^{0}u+2i\omega {\rho }^0{R}_{\Omega }u=\nabla \cdot T^{L1}-\nabla \cdot (u\cdot \nabla T^0)-{\rho }^{E1}\nabla {\Phi }^0-{\rho }^0\nabla S(u),$$

$${{\rm Y}}_{ijkl}^{T^0}=c_{ijkl}+\frac{1}{2}(-T_{ij}^0{\delta }_{kl}+T_{kl}^0{\delta }_{ij}+T_{ik}^0{\delta }_{jl}-T_{jl}^0{\delta }_{ik}+T_{jk}^0{\delta }_{il}-T_{il}^0{\delta }_{jk}),$$

$${\omega }^2{\rho }^{0}u-2i\omega {\rho }^0{R}_{\Omega }u=-\nabla \cdot (c:\nabla u)-\nabla ({\rho }^{0}u\cdot g^{\prime })+\nabla \cdot ({\rho }^{0}u)g^{\prime }+{\rho }^0\nabla S(u).$$

$$b_H(u,v)={\int }_{{\Omega }^{\text{S}}}({\overline{v}}^s\cdot u^s){\rho }^0\mathrm{d}x+{\int }_{{\Omega }^{\text{F}}}({\overline{v}}^f\cdot u^f){\rho }^0\mathrm{d}x,$$

$$c_r(u,v)={\int }_{{\Omega }^{\text{S}}}{\overline{v}}^s\cdot (\Omega \times u^s){\rho }^0\mathrm{d}x+{\int }_{{\Omega }^{\text{F}}}{\overline{v}}^f\cdot (\Omega \times u^f){\rho }^0\mathrm{d}x,$$

$$a_{2}(u,v)=\int_{\Omega^{\text{S}}}(\nabla\overline{v}^{s})\colon(c\colon\nabla u^{s})\,\mathrm{d}x+\int_{\Sigma^{\text{FS}}}\mathfrak{S}\{(\overline{v}^{s}\cdot g^{\prime})(\nu^{s\rightarrow f}\cdot u^{s})[\rho^{0}]^{f}\}\,\mathrm{d}\Sigma\\ +\int_{\Omega^{\text{S}}}\mathfrak{S}\big{\{}(\nabla\cdot\overline{v}^{s})(g^{\prime}\cdot u^{s})\rho^{0}-u^{s}\cdot(\nabla g^{\prime})\cdot\overline{v}^{s}\rho^{0}-u^{s}\cdot(\nabla\overline{v}^{s})\cdot g^{\prime}\rho^{0}\big{\}}\,\mathrm{d}x\\ +\int_{\Omega^{\text{F}}}\rho^{0}N^{2}\frac{(g^{\prime}\cdot\overline{v}^{f})(g^{\prime}\cdot u^{f})}{\|g^{\prime}\|^{2}}\,\mathrm{d}x+\int_{\Sigma^{\text{FF}}}(g^{\prime}\cdot\nu)(\overline{v}^{f}\cdot\nu)(u^{f}\cdot\nu)[\rho^{0}]^{+}_{-}\,\mathrm{d}\Sigma\\ +\int_{\Omega^{\text{F}}}\kappa(\nabla\cdot\overline{v}^{f}+\rho^{0}\kappa^{-1}g^{\prime}\cdot\overline{v}^{f})(\nabla\cdot u^{f}+\rho^{0}\kappa^{-1}g^{\prime}\cdot u^{f})\,\mathrm{d}x\\ -\frac{1}{4\pi G}\int_{\mathbb{R}^{3}}\nabla S(\overline{v})\cdot\nabla S(u)\,\mathrm{d}x,$$

$$a_2(u,v)={\omega }^2{b}_H(u,v)-i\omega c_r(u,v).$$

$$-p{\kappa }^{-1}=\nabla \cdot u^f+{\rho }^0{\kappa }^{-1}{g}^{\prime }\cdot u^f\text{in}{\Omega }^{\text{F}}.$$

$$0$$

$$-{\int }_{{\Sigma }^{\text{FS}}}{\overline{v}}^p({\nu }^{f\to s}\cdot u^s)\mathrm{d}\Sigma -{\int }_{\partial {\tilde{X}}^{\text{F}}}{\overline{v}}^p(\nu \cdot u^f)\mathrm{d}\Sigma ,$$

$$[\nu \cdot (\kappa \nabla \cdot u^f)]{|}_{\partial {\tilde{X}}^{\text{F}}}=0,$$

$$(\nu \cdot u^f){|}_{\partial {\tilde{X}}^{\text{F}}}=-\Vert g^{\prime }{\Vert }^{-1}(g^{\prime }\cdot u^f){|}_{\partial {\tilde{X}}^{\text{F}}}=({\rho }^0\Vert g^{\prime }\Vert {)}^{-1}p{|}_{\partial {\tilde{X}}^{\text{F}}}.$$

$$0=-\int_{\Omega^{\text{F}}}\overline{v}^{p}p\kappa^{-1}\,\mathrm{d}x+\int_{\Omega^{\text{F}}}\Big{[}(\nabla\overline{v}^{p})\cdot u^{f}-\overline{v}^{p}(g^{\prime}\cdot u^{f})\rho^{0}\kappa^{-1}\Big{]}\,\mathrm{d}x\\ -\int_{\Sigma^{\text{FS}}}\overline{v}^{p}(\nu^{f\rightarrow s}\cdot u^{s})\,\mathrm{d}\Sigma-\int_{\partial\tilde{X}^{\text{F}}}(\rho^{0}\|g^{\prime}\|)^{-1}\overline{v}^{p}p\,\mathrm{d}\Sigma=:c_{g}([u,p],v^{p}).$$

$$\begin{gathered} {\int }_{{\Omega }^{\text{F}}}\kappa (\nabla \cdot {\overline{v}}^f+{\rho }^0{\kappa }^{-1}{g}^{\prime }\cdot {\overline{v}}^f)(\nabla \cdot u^f+{\rho }^0{\kappa }^{-1}{g}^{\prime }\cdot u^f)\mathrm{d}x \\ ={\int }_{{\Omega }^{\text{F}}}[({\overline{v}}^f\cdot \nabla p)-({\overline{v}}^f\cdot g^{\prime })p{\rho }^0{\kappa }^{-1}]\mathrm{d}x-{\int }_{{\Sigma }^{\text{FS}}}({\overline{v}}^f\cdot {\nu }^{f\to s})p\mathrm{d}\Sigma . \end{gathered}$$

$$-{\int }_{{\Sigma }^{\text{FS}}}({\overline{v}}^f\cdot {\nu }^{f\to s})p\mathrm{d}\Sigma ={\int }_{{\Sigma }^{\text{FS}}}({\overline{v}}^s\cdot {\nu }^{s\to f})p\mathrm{d}\Sigma ,$$

$$\begin{gathered} a_2^f([u,p],v)={\int }_{{\Omega }^{\text{F}}}{\rho }^0{N}^2\frac{(g^{\prime }\cdot {\overline{v}}^f)(g^{\prime }\cdot u^f)}{\Vert g^{\prime }{\Vert }^2}\mathrm{d}x+{\int }_{{\Omega }^{\text{F}}}{\overline{v}}^f\cdot (\nabla p-g^{\prime }p{\rho }^0{\kappa }^{-1})\mathrm{d}x \\ +{\int }_{{\Sigma }^{\text{FF}}}(g^{\prime }\cdot \nu )({\overline{v}}^f\cdot \nu )(u^f\cdot \nu )[{\rho }^0{]}_{-}^{+}\mathrm{d}\Sigma . \end{gathered}$$

$$\begin{gathered} a_2^s(u,v)={\int }_{{\Omega }^{\text{S}}}(\nabla {\overline{v}}^s):(c:\nabla u^s)\mathrm{d}x \\ +{\int }_{{\Omega }^{\text{S}}}\mathfrak{S}\{(\nabla \cdot {\overline{v}}^s)(g^{\prime }\cdot u^s){\rho }^0-u^s\cdot (\nabla g^{\prime })\cdot {\overline{v}}^s{\rho }^0-u^s\cdot (\nabla {\overline{v}}^s)\cdot g^{\prime }{\rho }^0\}\mathrm{d}x \\ +{\int }_{{\Sigma }^{\text{FS}}}\mathfrak{S}\{({\overline{v}}^s\cdot g^{\prime })({\nu }^{s\to f}\cdot u^s)[{\rho }^0{]}^f\}\mathrm{d}\Sigma +{\int }_{{\Sigma }^{\text{FS}}}({\overline{v}}^s\cdot {\nu }^{s\to f})p\mathrm{d}\Sigma . \end{gathered}$$

$$\begin{gathered} a_G(u,v)=-\frac{1}{4\pi G}{\int }_{{\mathbb{R}}^3}\nabla S(\overline{v})\cdot \nabla S(u)\mathrm{d}x= \\ -{\int }_{{\Omega }^{\text{S}}}\nabla \cdot ({\rho }^0{\overline{v}}^s)S(u)\mathrm{d}x-{\int }_{{\Sigma }^{\text{SS}}\cup \partial {\tilde{X}}^{\text{S}}}(\nu \cdot {\overline{v}}^s)S(u)[{\rho }^0{]}_{-}^{+}\mathrm{d}\Sigma \\ -{\int }_{{\Omega }^{\text{F}}}\nabla \cdot ({\rho }^0{\overline{v}}^f)S(u)\mathrm{d}x-{\int }_{{\Sigma }^{\text{FF}}\cup \partial {\tilde{X}}^{\text{F}}}(\nu \cdot {\overline{v}}^f)S(u)[{\rho }^0{]}_{-}^{+}\mathrm{d}\Sigma \\ -{\int }_{{\Sigma }^{\text{FS}}}\left\{({\nu }^{f\to s}\cdot {\overline{v}}^s)S(u)[{\rho }^0{]}^s+({\nu }^{s\to f}\cdot {\overline{v}}^f)S(u)[{\rho }^0{]}^f\right\}\mathrm{d}\Sigma , \end{gathered}$$

$$S(u)=G{\int }_{\tilde{X}}\frac{{\nabla }^{\prime }\cdot ({\rho }^0(x^{\prime })u(x^{\prime }))}{\Vert x-x^{\prime }\Vert }\mathrm{d}{x}^{\prime }+G{\int }_{\Sigma \cup \partial \tilde{X}}\frac{\nu (x^{\prime })\cdot u(x^{\prime })[{\rho }^0(x^{\prime }){]}_{-}^{+}}{\Vert x-x^{\prime }\Vert }\mathrm{d}{\Sigma }^{\prime }.$$

$$S(u)=G\bigg{\{}\int_{\Omega^{\text{S}}}\frac{\nabla^{\prime}\cdot(\rho^{0}(x^{\prime})u^{s}(x^{\prime}))}{\|x-x^{\prime}\|}\,\,\mathrm{d}x^{\prime}+\int_{\Omega^{\text{F}}}\frac{\nabla^{\prime}\cdot(\rho^{0}(x^{\prime})u^{f}(x^{\prime}))}{\|x-x^{\prime}\|}\,\,\mathrm{d}x^{\prime}\\ +\int_{\Sigma^{\text{SS}}\cup\partial\tilde{X}^{\text{S}}}\frac{\nu(x^{\prime})\cdot u^{s}(x^{\prime})[\rho^{0}(x^{\prime})]^{+}_{-}}{\|x-x^{\prime}\|}\,\,\mathrm{d}\Sigma^{\prime}+\int_{\Sigma^{\text{FF}}\cup\partial\tilde{X}^{\text{F}}}\frac{\nu(x^{\prime})\cdot u^{f}(x^{\prime})[\rho^{0}(x^{\prime})]^{+}_{-}}{\|x-x^{\prime}\|}\,\,\mathrm{d}\Sigma^{\prime}\\ +\int_{\Sigma^{\text{FS}}}\frac{[\rho^{0}(x^{\prime})]^{s}\nu^{f\rightarrow s}(x^{\prime})\cdot u^{s}(x^{\prime})+[\rho^{0}(x^{\prime})]^{f}\nu^{s\rightarrow f}(x^{\prime})\cdot u^{f}(x^{\prime})}{\|x-x^{\prime}\|}\,\,\mathrm{d}\Sigma^{\prime}\bigg{\}}.$$

$$a_2([u,p],v)=a_2^s(u,v)+a_2^f([u,p],v)+a_G(u,v).$$

$$\{\begin{array}{ll}{a}_2([u,p,S(u)],v)& ={\omega }^2{b}_H(u,v)-2i\omega c_r(u,v),\\ c_g([u,p],v^p)& =0.\end{array}$$

$$a_2([u,p(u^f,u_{{\Sigma }^{\text{FS}}}^s),S(u)],v)={\omega }^2{b}_H(u,v)-2i\omega c_r(u,v).$$

Full text

Declarations

This research was supported by the Simons Foundation under the MATH+X program, the National Science Foundation grant DMS-1815143, the members of the Geo-Mathematical Imaging Group at Rice University, and XSEDE research allocation TG-EAR170019. The work by R.L. was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344 (LLNL-JRNL-780818). Y.X. and Y.S. were supported by NSF-1812695.

The codes are made available via https://github.com/js1019/NormalModes and https://github.com/eigs/pEVSL. The data can be reproduced using the codes in https://github.com/js1019/PlanetaryModels. In addition, the Mars models can be found in agwx-jd58-20 , where the performance and reproducibility were studied in shi2021planetary .

1 Introduction

Planetary normal modes are instrumental for studying the dynamic response to sources including earthquakes along faults and meteorite impacts, as well as tidal forces dahlen1998theoretical ; lognonne2005planetary . The low-angular-order eigenfrequencies contain critical information about the planet’s large-scale structure and provide constraints on heterogeneity in composition, temperature, and anisotropy, while rotation constrains the shapes as well as possible density distributions of planets. The effect of rotation on the seismic point spectrum of the Earth is well understood and has been observed for decades (park2005earth, , Fig.1). The observation of spectral energy of low-frequency toroidal modes in vertical seismic recordings of the 1998 Balleny Islands earthquake zurn2000observation , is a manifestation of the three-dimensional heterogeneity and anisotropy of the mantle structures and rotation.

For a review of Earth’s free oscillations, we refer to woodhouse2007 . Current standard approaches to computing the seismic point spectrum and associated normal modes have several limitations. Assuming spherical symmetry for non-rotating planets, the problem becomes one-dimensional and the computation of normal modes in such models using MINEOS woodhouse1988calculation ; masters2011mineos is still common practice; these are then typically used in perturbation-theory and mode-coupling approaches to include lateral heterogeneities. Full-mode coupling methodology utilizing normal modes in a spherically symmetric model as a basis has been adopted to studying Earth’s interior for decades dahlen1968normal ; dahlen1969normal ; woodhouse1978effect ; woodhouse1980coupling ; park1986synthetic ; park1990subspace ; romanowicz1987multiplet ; lognonne1990modelling ; hara1991inversion ; hara1993inversion ; um1991normal ; lognonne1991normal ; deuss2001theoretical ; deuss2004iteration ; al2012calculation ; yang2015synthetic . This methodology is of Rayleigh-Ritz type, and is justified under the assumption that the space in which the normal modes lie contains the mentioned basis, which requires spherically symmetric fluid-solid and surface boundaries. Here, we remove this limitation. Moreover, a separation of the essential spectrum needs to be carefully carried out, which has been commonly ignored in the “radial” algorithms. We discuss the mode-coupling approach and the conditions under which it applies in Appendix B.

To simulate seismic waves in strongly heterogeneous media, the spectral-element method (SPECFEM) komatitsch1998spectral ; komatitsch1999introduction has been widely used for more than two decades. We mention the software package SPECFEM3D_globe komatitsch2002spectral1 ; komatitsch2002spectral2 , which is capable of modeling relatively high-frequency waveforms in an entire planet while suppressing the perturbation to the gravitational potential. Other implementations of SPECFEM chaljub2003solving ; chaljub2004spectral ; chaljub2007spectral have been developed with alternative numerical approaches pertaining to the fluid outer core. In principle, seismic eigenfrequencies show up by taking a discrete Fourier transform of numerical solutions; however, it is a major computational challenge to control the accuracy at very long time scales. We note that in SPECFEM3D_globe, the fluid displacement is replaced by a scalar potential, which results in a non-symmetric system of discretized equations. Moreover, the (square of the) Brunt-Väisälä frequency is assumed to be zero. Rotation in the fluid regions is unnaturally introduced by means of an additional vector (cf. (komatitsch2002spectral2, , (16) and (17)) and (chaljub2007spectral, , (30))). In addition, current SPECFEM3D_globe does not include the incremental gravitational field, which limits its usage for relatively higher frequency wave propagation.

One may view the computational approach developed in this paper as forming a bridge between SPECFEM3D_globe, and the mode-coupling approaches derived from modes in a spherically symmetric model, involving finer scale heterogeneity and higher seismic eigenfrequencies. Our approach facilitates the studies of the highly heterogeneous crust models and complex three-dimensional models through the planetary spectrum, as well as the naturally efficient computation of seismograms from many different sources. Naturally, we also include the Coriolis force and centrifugal potential and formulate it as a nonlinear eigenvalue problem. We can accommodate arbitrarily shaped fluid-solid boundaries which becomes increasingly important at higher rotation rates. In our formulation, the rotation rate might spatially vary, which is relevant to the future computation of normal modes in gas giants in our solar system.

In this paper, we revisit the work of buland1984computation . Buland and collaborators encountered several complications that we overcome by characterizing and separating the essential spectrum using a polynomial filtering eigensolver and introducing a new formulation that properly models the elastic-gravitational system without simplifications. In our proposed formulation, the displacement, the proper orthonormal condition and the symmetry of the system for non-rotating planets are preserved. We apply fully unstructured meshes to model fully heterogeneous planets, and the mixed finite-element method (FEM) to discretize the elastic-gravitational system. Our method can handle fully heterogeneous planetary models easily, and guarantee that accurate solutions lie in the space to which normal modes associated with the seismic point spectrum belong. In a previous paper DBLP:conf/sc/ShiLXSH18 , we introduced a highly parallel algorithm for solving the generalized eigenvalue problem resulting from our analysis for Cowling approximation using P1 mixed FEM. We achieved high parallel computational and memory scalabilities with demonstrated performance on modern supercomputers. In the following paper shi2021planetary , we extended our algorithm using P2 mixed FEM for better accuracy and discussed the reproducibility of our codes reported from several universities during the student cluster competion at the supercomputing conference.

Self-gravitation manifests itself in the incremental gravitational potential as the density changes with displacement. We utilize the Green’s solution of Poisson’s equation and treat the self-gravitation as an $N$-body problem. We then apply the fast multipole method (FMM) greengard1997new ; gimbutas2011fmmlib3d ; yokota2013fmm , which reduces the algorithmic complexity significantly, to compute both the reference gravitational and the incremental gravitational potentials. Alternatively, one can apply a finite-infinite element method zienkiewicz1983novel ; burnett1994three for modeling unbounded domain problems to approximate the far-field of Poisson’s equation. More recently, the spectral-infinite-element method gharti2018spectral has been developed to incorporate gravity. While our eigensolver DBLP:conf/sc/ShiLXSH18 only takes matrix-vector products, any suitable schemes, including FMM or infinite-element methods, can be used in our computational framework.

To include rotation in the elastic-gravitational system through the Coriolis force and the centrifugal potential, in this work, we utilize extended Lanczos vectors computed in a non-rotating planet – with the shapes of boundaries of a rotating planet and accounting for the centrifugal potential – as a truncated basis to properly facilitate reduction to one of the equivalent linear forms of the quadratic eigenvalue problem (QEP). Here, the separation of the essential spectrum comes into play again and the normal modes computed are guaranteed to lie in the appropriate space of functions. The reduced system can be solved with a standard eigensolver.

We present and validate our three-dimensional computations using constant elastic balls and an isotropic preliminary reference Earth model as non-rotating planets with standard radial codes. The computational accuracy for rotating planets is illustrated and tested but only in the regime where perturbation theory applies as no other independent algorithms are available in the general case. We use our algorithm and code to compute the point spectra of eigenfrequencies in several Earth and Mars models, acknowledging relatively low rotation rates, studying the effects of heterogeneity on a large range of scales. The Mars models are relevant to the InSight (Interior exploration using Seismic Investigations, Geodesy and Heat Transport) banerdt2013insight ; lognonne2019seis mission. It is expected that a set of eigenfrequencies is observable panning2017planned ; bissig2018detectability . Here, we select one Mars model khan2016single from the set of blind tests clinton2017preparing ; van2019preparing and combine it with the topography zuber1992mars ; smith1999global and a three-dimensional crust belleguic2005constraints ; goossens2017evidence to create a realistic Mars model. We compute the low-angular-order eigenfrequencies and study the general effects of rotation and heterogeneity combined.

The outline of this paper is as follows. In Section 2, we revisit the form and physics of the elastic-gravitational system of a rotating planet and establish the weak formulation of the system with a separation of the essential spectrum using a polynomial filtering eigensolver. In Section 3, we discuss the hydrostatic equilibrium of a rotating fluid outer core in the presence of the gravitational and the centrifugal forces. In Section 4, we introduce the Continuous Galerkin mixed FEM and obtain the corresponding matrix equations. In Section 5, we study the computation of the reference gravitational field and the perturbation of the gravitational field using the FMM. In Section 6, we validate the computational accuracy of our work for non-rotating Earth models and quantify the effect on the point spectrum from three-dimensional heterogeneity. In Section 7, we illustrate the computational accuracy of our proposed method and show several computational experiments for different planetary models, including standard Earth and Mars models as well as related effects due to rotation and a three-dimensional crust. In Section 8, we discuss the significance of our results and directions of future research.

2 The elastic-gravitational system with rotation

In this section, we present a modified elastic-gravitational system of equations of a rotating planet to deal with the separation of the essential spectrum in the weak form de2015system (see dahlen1998theoretical for the strong formulation).

2.1 Natural subdomains and computational meshes

Following the notation in de2015system , a bounded set $\tilde{X}\subset\mathbb{R}^{3}$ is used to represent the interior of the Earth, with Lipschitz continuous exterior boundary $\partial\tilde{X}$. The exterior boundary $\partial\tilde{X}$ contains fluid (ocean) surfaces $\partial\tilde{X}^{\text{F}}$ and solid surfaces $\partial\tilde{X}^{\text{S}}$. We subdivide the set $\tilde{X}$ into solid regions $\Omega^{\text{S}}$ and fluid regions $\Omega^{\text{F}}$. The fluid regions contain the liquid outer core $\Omega^{\text{OC}}$ and the oceans $\Omega^{\text{O}}$. The solid regions can be further subdivided into the crust and mantle $\Omega^{\text{CM}}$ and the inner core $\Omega^{\text{IC}}$. We use $\Sigma$ to represent the interfaces between these subregions. In summary,

[TABLE]

The interior interfaces can further be subdivided into three categories: interfaces between two fluid regions $\Sigma^{\text{FF}}$, interfaces between two solid regions $\Sigma^{\text{SS}}$, and interfaces between fluid and solid regions $\Sigma^{\text{FS}}$. We can subdivide $\Sigma^{\text{FS}}$ into two major interfaces: internal interfaces $\Sigma^{\text{FS}}_{\text{int}}$ and the bottom interface $\Sigma^{\text{FS}}_{\text{O}}$ of the oceans. The internal interfaces include the interfaces between the lower mantle and the outer core $\Sigma^{\text{CMB}}$, which is known as the Core-Mantle Boundary (CMB); the interface between the outer core and the inner core is denoted as $\Sigma^{\text{ICB}}$, which is known as the Inner-Core Boundary (ICB). Thus,

[TABLE]

In Fig. 1, we illustrate the concepts of the main mathematical symbols for the geometry used in this work. Since a general terrestrial planet may contain multiple complex discontinuities associated with different geological and geodynamical features, utilization of a flexible, fully unstructured tetrahedral mesh would be natural. We discretize the major discontinuities using triangulated surfaces that are generated via distmesh persson2004simple and then build up the Earth model using an unstructured tetrahedral mesh via TetGen si2015tetgen . In Fig. 2, we illustrate the interfaces and meshes with one hundred thousand and one million elements. These techniques show great flexibility and can provide models with multiple resolutions. In Figs. 3, we illustrate a three-dimensional Earth model built on a tetrahedral mesh. In Fig. 3 (a), we show the Moho interface that is constructed using an unstructured triangular mesh. The color shows the depth and the black lines are the edges of the triangles. In Fig. 3 (b), we illustrate the three-dimensional $V_{P}$ model based on MIT’s mantle tomographic results burdick2017model and crust 1.0 laske2013update . The core model is based on the Preliminary Reference Earth Model (PREM) dziewonski1981preliminary .

We also use a Mars model as an example to illustrate our construction of a terrestrial planet. The topography of Mars was measured by the Mars Orbiter Laser Altimeter (MOLA) zuber1992mars ; smith1999global with high accuracy. The thickness and density of the Martian crust were constructed with the help of the works of belleguic2005constraints ; goossens2017evidence . In Fig. 4 (a), we illustrate the topography of Mars using data from MOLA smith1999global ; in Fig. 4 (b), we show the crust-mantle interface of Mars using data provided by goossens2017evidence . In Figs. 5 (a)–(c), we illustrate $V_{P}$, $V_{S}$ and $\rho^{0}$ of Mars integrating a radial model khan2016single with a three-dimensional crust as shown in Fig. 4. In Figs. 6 (a) and (b), we illustrate the axial spin mode, $\Omega\times x$, and the centrifugal acceleration, $-\nabla\psi$, of the Mars model, respectively.

2.2 The basic equations

Given the reference density $\rho^{0}$ and the gravitational constant $G$, we let $\Phi^{0}$ denote the gravitational potential which satisfies,

[TABLE]

and $S(u)$ denote the Eulerian perturbation of the Newtonian potential associated with the displacement $u$,

[TABLE]

To include the centrifugal force, we introduce the centrifugal potential

[TABLE]

where $\Omega\in\mathbb{R}^{3}$ is the angular velocity of rotation. We form the gradient,

[TABLE]

where the reference gravitational field

[TABLE]

The initial stress $T^{0}$ satisfies the mechanical equilibrium given by the static momentum equations,

[TABLE]

The elastic-gravitational system of a rotating non-hydrostatic terrestrial planet has the form

[TABLE]

where $\omega$ denotes the angular frequency; $R_{\Omega}u=\Omega\,\times\,u$; $\rho^{\mathrm{E1}}=-\nabla\cdot(\rho^{0}u)$ denotes the first-order Eulerian density perturbation and $T^{\mathrm{L1}}=\Upsilon^{T^{0}}\colon\nabla u$ denotes the incremental Lagrangian Cauchy stress. The elasticity tensor, $\Upsilon^{T^{0}}_{ijkl}$, attains the form,

[TABLE]

where $c$ denotes the elastic stiffness tensor. In fact, (6) does not determine the entire tensor $T^{0}$. It is common practice to invoke the hydrostatic assumption when $T_{ij}^{0}=-p^{0}\delta_{ij}$; then $\Upsilon^{T^{0}}_{ijkl}$ reduces to $c_{ijkl}$. Under the hydrostatic assumption, we reduce (7) into

[TABLE]

The boundary conditions for the system (8) governing a hydrostatic planet are summarized in Table 1.

2.3 The weak formulation

We let $u^{s}$ denote displacement in the solid regions and $u^{f}$ denote displacement in the fluid regions. We treat the solid and fluid parts differently and then deal with $S(u)$ globally. We use $v$ to denote test functions and denote $v^{s}$ and $v^{f}$ for the solid and fluid test displacements, respectively. The mass term from the first and the second term of (8) take the form

[TABLE]

and

[TABLE]

respectively. We note that the coercivity of the original weak form of the right-hand side of (8), identified as $a_{\text{original}}(u,v)$ in (de2015system, , (3.5)), is not apparent. The early work by Valette valette1989spectre , which is written in French, analyzed this problem in a proper mathematical space while the details can be found in a preprint of a book chapter de2015system . In the work of de2015system , it is revisited and a proper form, $a_{2}(u,v)$, for the weak formulation is introduced. The coercivity of $a_{2}(u,v)$ is established in (de2015system, , Sections 5.2 and 6). The equivalence, that is, $a_{2}(u,v)=a_{\text{original}}(u,v)$ under the boundary conditions (cf. (dahlen1998theoretical, , Table 3.1)), is proven in (de2015system, , Lemma 4.1).

In this work, we will study $a_{2}(u,v)$ under the hydrostatic assumption. The right hand side of (8) can be written in the form

[TABLE]

where $N^{2}=(\nabla\rho^{0}/\rho^{0}-g^{\prime}\rho^{0}/\kappa)\cdot g^{\prime}$ signifies the square of the Brunt-Väisälä frequency; $\nu^{s\rightarrow f}$ denotes the normal vector at the fluid-solid boundary pointing from the solid to the fluid side; the symmetrization operation $\mathfrak{S}$ is defined as $\mathfrak{S}\{L(u,\overline{v})\}:=\frac{1}{2}(L(u,\overline{v})+L(\overline{v},u))$, for any bilinear form $L(u,v)$. The first integral over $\Omega^{\text{F}}$ is responsible for the inertial or gravity modes, and the second integral over $\Omega^{\text{F}}$ yields the acoustic modes. The integral over $\Sigma^{\text{FF}}$ generates Kelvin modes that occur at boundaries with density jumps. To solve the basic equation (8), we combine (9) with (11) and obtain the system

[TABLE]

However, it is computationally infeasible to obtain the accurate normal modes from the direct discretization of (12) due to the existence of spurious oscillations kiefling1976fluid . We discuss various approaches in Subsection 2.3.1 and note that the solution needs to be restricted to the space associated with the seismic point spectrum. In Subsections 2.3.2, 2.3.3, 2.3.4 and 2.3.5, we present our scheme to deal with the fluid-solid and fluid surface boundary conditions, fluid regions, solid regions and perturbation of the gravitational potential and field, respectively. In the Subsection 2.4, we introduce the mathematical spaces associated with the seismic point and essential spectra and their separation using a polynomial filtering eigensolver.

2.3.1 Choice of physical variables for fluid regions without rotation

To study planetary normal modes, we include the linear elasticity, compressible fluids, and the fluid-solid and free-surface boundary conditions. Discretization of the standard formulation leads to computational difficulties, since the non-seismic modes from the compressible fluid may pollute the computation of the point spectrum. In this paper, we use a displacement-pressure formulation and later substitute the pressure term using an equivalent formula.

Here, we review different approaches pertaining to the above-mentioned separation of the essential spectrum for non-rotating bodies and then include the rotation. The natural displacement formulation for a non-rotating body will result in a symmetric eigenvalue problem. However, the drawback is the existence of spurious oscillations kiefling1976fluid . Several finite-element methods have been developed for modeling the fluid regions with fluid-solid interaction: a displacement formulation hamdi1978displacement , a pressure formulation zienkiewicz1969coupled ; craggs1971transient , a displacement-pressure formulation wang1997displacement , and a velocity potential formulation everstine1981symmetric ; olson1985analysis . However, the pressure formulation leads to a non-symmetric eigenvalue problem zienkiewicz1969coupled ; craggs1971transient , and the velocity potential formulation everstine1981symmetric ; olson1985analysis leads to a quadratic eigenvalue problem.

In the engineering community, several approaches have been designed to resolve this issue. A penalty method hamdi1978displacement has been applied by imposing an irrotational constraint. However, the study by olson1983study has shown that this penalty method has issues dealing with a solid vibrating in the fluid cavity, which is the case in this paper. A four-node element with a reduced integration using a mass matrix projection technique chen1990vibration has been designed to eliminate the spurious modes. A method using different elements for solid and fluid regions was proposed for two-dimensional bermudez1994finite and three-dimensional cases bermudez1999finite when non-physical spurious modes appear bermudez1995finite . The displacement/pressure formulation wang1997displacement has been developed via introducing mixed elements; still, the fluid-solid coupling needs additional consideration bermudez1994finite ; bermudez1999finite .

Compared with the above-mentioned engineering problems, we encounter a more complicated system (8) with different boundary conditions (cf. Table 1). Due to the presence of the reference gravitational field and the incremental gravitational field, the essential spectrum of the elastic-gravitational system is more complicated than the one of the elastic systems with fluid structures in the engineering problems. In the geophysical community, the pressure formulation komatitsch2002spectral1 ; komatitsch2002spectral2 ; nissen20082 has been commonly used, which is based on replacing the displacement by a scalar potential in the fluid regions. It results in non-symmetric stiffness and mass matrices for a non-rotating body. An alternative approach chaljub2003solving ; chaljub2004spectral ; chaljub2007spectral , using several additional variables to represent the fluid displacement, also leads to a non-symmetric system. To preserve the necessary symmetry and guarantee the correct orthonormality condition for the eigenfunctions or normal modes, we note that the fluid displacement must be kept in the formulation.

2.3.2 Fluid-solid and fluid surface boundary conditions

In this work, to deal with fluid-solid and fluid surface boundary conditions we applied a similar approach wang1997displacement with no any penalty terms by augmenting the system of equations (cf. (11)) and introducing an additional variable, $p$, according to

[TABLE]

Here, $\kappa$ signifies the compressibility of the fluid. Imposing the fluid-solid boundary condition $\left.\left[\nu^{f\rightarrow s}\cdot u^{f}-\nu^{f\rightarrow s}\cdot u^{s}\right]\right|_{\Sigma^{\text{FS}}}=0$ naturally with the introduction of the additional variable $p$, we obtain the weak form for (13),

[TABLE]

for all the test functions $v^{p}$, where $\nu^{f\rightarrow s}$ denotes the normal vector at the fluid-solid boundary pointing from the fluid to the solid side. Due to the hydrostatic equilibrium, we note that $\nu|_{\partial\tilde{X}^{\text{F}}}$ is parallel to $g^{\prime}$. Using the boundary condition,

[TABLE]

we have the relation

[TABLE]

We using (16) to rewrite (14)

[TABLE]

A short-hand notation $c_{g}([u,p],v^{p})$ in (17) is introduced for simplification. In this work, since we only consider planets with a solid surface, the integral over $\partial\tilde{X}^{\text{F}}$ will be omitted. But it will be needed while including the oceans, or dealing with gas giants, such as Saturn or Jupiter.

2.3.3 Fluid regions

We use (13) in (11) and obtain

[TABLE]

Since

[TABLE]

we include the right-hand side of (19) in the contributions from the solid regions. Thus, we obtain the contributions to $a_{2}(u,v)$ in (11) from the fluid regions,

[TABLE]

2.3.4 Solid regions

For the solid regions, we add the right-hand side of (19) to the terms related to the solid regions in (11) and obtain

[TABLE]

2.3.5 Perturbation of the gravitational potential and field

Here, we discuss the contribution of the perturbation of the gravitational potential $S(u)$. Since the test functions are divided into test functions on solid and fluid regions, we have

[TABLE]

where $[\rho^{0}]^{s}$ denotes the solid density along the fluid-solid boundary. One can set up $S(u)$ as an independent variable and apply the finite-infinite element method to approximate (2), but here we follow a different approach.

Making use of Green’s function (dahlen1998theoretical, , Chapter 3, (3.98)), we have

[TABLE]

Again, we separate the displacement $u$ into $u^{s}$ and $u^{f}$, and rewrite (23) as

[TABLE]

Although we impose $\nu^{s\rightarrow f}\cdot u^{f}=\nu^{s\rightarrow f}\cdot u^{s}$ along the fluid-solid boundaries, we keep the construction of the incremental gravitational potential $S(u)$ as described in (24). This is to preserve the symmetry of the bilinear form as we substitute (24) into (22).

Since the Green’s solution is known, we apply the FMM to evaluate $S(u)$ for a given displacement $u$ via (24). The utilization of this approach is computationally attractive, but requires that the eigensolver can solve for the interior eigenpairs via matrix-vector multiplications.

2.3.6 Summary

To restrict the system to the computational domain, we can rewrite (11) as

[TABLE]

We obtain the complete formula for the rotating hydrostatic planetary model (9), (10), (25) and (17):

[TABLE]

A matrix representation can be derived from (26). In practice, we replace $p$ in $a_{2}$ by $p(u^{f},u^{s}_{\Sigma^{\text{FS}}})$ via solving the constraint $c_{g}([u,p],v^{p})=0$ in (17) and obtain

[TABLE]

The corresponding orthonormality condition is that, for an eigenpair $(\omega_{(i)},u^{(i)})$, any other eigenpair $(\omega_{(j)},u^{(j)})$ satisfies

[TABLE]

which is consistent with (dahlen1998theoretical, , (4.82)).

2.4 Hilbert space for the elastic-gravitational system

We introduce the space for the displacement field (de2015system, , Definition 5.4)

[TABLE]

where

[TABLE]

$L^{2}(\tilde{X},\rho^{0}\,\mathrm{d}x)$ denotes a weighted $L^{2}$ Hilbert space with

[TABLE]

We write $H=L^{2}(\tilde{X},\rho^{0}\,\mathrm{d}x)$ subject to the constraint $\int_{\tilde{X}}u\rho^{0}\,\mathrm{d}x=0$ removing rigid-body translations; $E$ is densely embedded in $H$ de2019note .

To describe the essential spectrum, we introduce operator $T$ in (valette1989spectre, , Section 4) and de2019note ,

[TABLE]

The adjoint, $T^{*}$, of $T$ is given by

[TABLE]

where $\varphi$ has the interpretation of potential. A subspace, $H_{2}$, of $H$ associated with the essential spectrum is defined by the constraints

[TABLE]

In fact, $u^{f}$ can be decomposed according to $\text{Ran}(T^{*})\oplus\text{Ker}(T)$, following the decomposition

[TABLE]

where spaces $H_{1}$ and $H_{2}$ are associated with the point and essential spectrum, respectively. The space $H_{2}$ is designed precisely to extract, via projections, the “subseismic” approximations to the full system of governing equations for a contained rotating, compressible, inhomogeneous, self-gravitating fluid. The rigid boundary condition, $u|_{\Omega^{\text{F}}}\cdot\nu=0$ on $\Sigma^{\text{FF}}\cup\partial\tilde{X}^{\text{F}}$, is consistent with a rigid mantle and rigid inner core as $u|_{\Omega^{\text{S}}}=0$.

In fact, $\forall u\in H_{2}$, we obtain $p=0$ and for Cowling approximation, we have

[TABLE]

where $a_{2}^{s}$ and $a_{2}^{f}$ are defined in (21) and (20), respectively. For the incremental gravitational potential in (2), we have

[TABLE]

Combining (33) and (34), we note that (27) will be reduced to

[TABLE]

Thus, restricting $\forall u\in H_{2}$, the associated spectrum of (35) will essentially depend on to $N^{2}$ and the rotating rates.

In this work, we solve for the eigenvalues and eigenfunctions of (27) inside a target frequency interval $[f_{1},f_{2}]$, where

[TABLE]

where $N^{2}_{\sup}$ denotes the supremum of the square of the Brunt-Väisälä frequency. We note that inequality (36) holds true for most planets because the minimal seismic normal mode frequency is typically much larger than the upper bound of the associated spectrum of (35), which is the right hand side of (36). For instance, the maximum of the Brunt-Väisälä frequency of the Earth is around 50 $\mu$Hz and $|\Omega|$ is 7.3 $\mu$Hz while the minimal seismic normal mode frequency is around 0.3 mHz. A well-designed polynomial filter applied with the eigensolver, will have the effect of boosting up the eigenvalues inside the interval $[f_{1},f_{2}]$ while lessening the rest of the spectrum, including the part associated with $H_{2}$.

Remark 1

It is important to understand the need for polynomial filtering in this context. First note that eigensolvers like ARPACK lehoucq1998arpack or subspace iteration, e.g., saad2011numerical , compute eigenvalues of a matrix on one end of the spectrum. After discretizetion, the essential spectrum will give rise to a large number of eigenvalues near zero. Computing the (discrete) eigenvalues in the interval $[f_{1},f_{2}]$ will be numerically challenging unless the small eigenvalues associated with the essential spectrum are eliminated. In numerical linear algebra, this is termed an interior eigenvalue problem in that the target eigenvalues of the discretized problem are located well inside the spectrum. If we use a standard package like ARPACK lehoucq1998arpack we could compute these eigenvalues starting from the smallest ones until we reach the desired interval $[f_{1},f_{2}]$, which would be prohibitive because of the large cluster near zero caused by the essential spectrum. Alternatively, we could compute them from the largest ones down. This would also entail computing a large number of unwanted eigenpairs. Finally, we could also use a shift-and-invert strategy Parlett-book within ARPACK. This requires using a direct solver with a very large matrix and is impractical in our context due to the large memory requirement. The advantage of polynomial filtering is that it eliminates the unwanted eigenvalues and allows the eigensolver to focus on those that are amplified, namely those in $[f_{1},f_{2}]$.

In Section 3, we study the hydrostatic equilibrium of the liquid regions with rotation and derive a proper density distribution. In Section 4, we introduce the mixed FEM to construct the system without the perturbation of the gravitational field. In Section 5, we utilize FMM to compute the gravitational field and the perturbation of the gravitational field and then obtain the complete matrix formula for (27).

3 Hydrostatic equilibrium of the liquid core with rotation

In this section, we discuss the hydrostatic equilibrium with rotation and how it constrains the shape of the boundaries and the density distribution in planets. Rotating fluids have been extensively studied greenspan1968theory ; chandrasekhar2013hydrodynamic ; zhang2017theory . The outer core’s properties have been studied through seismic normal modes since the 1970s gilbert1975application ; dziewonski1975parametrically ; dziewonski1981preliminary , but also with body waves morelli1993body ; kennett1995constraints . Much more recently, an alternative radial outer core model has been proposed using the parametrization of the equation of state for liquid iron alloys at high pressures and temperatures, inferred from eigenfrequency observations irving2018seismically . Furthermore, we mention models for the core of Mars rivoldini2011geodesy ; khan2016single albeit ignoring rotation.

To reach the hydrostatic equilibrium, the prepressure $p^{0}$ satisfies

[TABLE]

where $g^{\prime}$ is defined in (4). Well-posedness requires that

[TABLE]

see (de2015system, , Lemma 2.1) for details about the functional properties of $\rho^{0}$, $p^{0}$ and $g^{\prime}$.

The derivation of Clairaut’s equation clairaut1743theorie , and Radau approximation are put in the context of a general scheme imposing (38) in (dahlen1998theoretical, , Chapter 14.1). The bulk parameters of Earth and Mars are listed in Table 2. While the hydrostatic assumption seems to apply to Earth with reasonable accuracy, the derivative of the ellipticity at $r_{e}$, $\dot{\epsilon}(r_{e})$, of Mars appears to be negative, whence this assumption fails to hold dollfus1972new ; bills1978mars .

To construct models of liquid planet interiors, such as Jupiter and Saturn, equations of state and theory of figures are commonly used for calculating a self-consistent shape and gravity field jeans1919problems . We refer to militzer2016understanding for a review on modelling Jupiter’s interior using equations of state and multiple mission data. Since Radau assumptions break down for fast rotating plants (wahltidal, , Fig.3), we refer to hubbard2013concentric ; militzer2019models for constructing Saturn’s interior using the concentric Maclaurin spheroid method to match the Cassini measurements. The condition (38) is satisfied along with other conditions.

4 The Continuous Galerkin mixed finite-element method

In this section, we employ the Continuous Galerkin mixed FEM zienkiewicz2005finite ; bathe2006finite ; hughes2012finite ; brezzi2012mixed ; ern2013theory , for discretizing our system without the perturbation of the gravitational field. We thus obtain a matrix representation for the corresponding weak forms. The incremental gravitational potential will be introduced in the discretization in Subsection 5.2.

4.1 The Continuous Galerkin mixed finite-element approximation

Given a shape regular finite-element partitioning $\mathcal{T}_{h}$ of the domain $\tilde{X}$, we denote an element of the mesh by $K_{k}\in\mathcal{T}_{h}$ and a boundary element by $E_{l}\subset\partial K_{k}$ and have

[TABLE]

where $N_{K}$ denotes the total number of volume elements and $N_{E}$ denotes the total number of interior and exterior boundary elements. Furthermore, we let $K_{k}^{\text{S}}$ and $K_{k}^{\text{F}}$ be elements in the solid and fluid regions, respectively. Similarly, $E_{l}^{\text{S}}$, $E_{l}^{\text{F}}$ and $E_{l}^{\text{FS}}$ denote boundary elements on the solid $\Sigma^{\text{SS}}\cup\partial\tilde{X}^{\text{S}}$, fluid $\Sigma^{\text{FF}}\cup\partial\tilde{X}^{\text{F}}$ and fluid-solid $\Sigma^{\text{FS}}$ discontinuities, respectively. We have

[TABLE]

with

[TABLE]

where $N_{K}^{\text{S}}$ and $N_{K}^{\text{F}}$ denote the total number of volume elements in the solid and fluid regions, respectively, and $N_{E}^{\text{S}}$, $N_{E}^{\text{F}}$ and $N_{E}^{\text{FS}}$ denote the total number of boundary elements on the (interior/exterior) solid, fluid and fluid-solid boundaries, respectively. In the above, $h$ signifies the maximum value of diameters of all the elements.

Since we separate out the fluid and solid regions, we divide the finite-element partitioning accordingly into

[TABLE]

where $\mathcal{T}_{h}^{\text{S}}$, $\mathcal{T}_{h}^{\text{F}}$ and $\Sigma^{\text{FS}}_{h}$ denote the partitioning of the domains $\Omega^{\text{S}}$, $\Omega^{\text{F}}$ and boundary $\Sigma^{\text{FS}}$, respectively. We then introduce $E_{h}$ as the finite-element space corresponding with the displacement space $E$ in (29),

[TABLE]

and $\mathbf{V}_{h}^{p}$ as the finite-element space for $p$,

[TABLE]

Here, $\mathcal{P}_{p^{s}}(K)$ and $\mathcal{P}_{p^{f}}(K)$ are the spaces of polynomials of degrees $p^{s}$ and $p^{f}$, respectively; $\mathcal{P}_{p^{p}}(K)$ is the space of polynomials of degree $p^{p}$. Though the $u_{h}^{f}$ is discretized as $u_{h}^{f}\in H^{1}(\Omega^{\text{F}})$, the constraint equation (13) restricts $u_{h}^{f}\in H(\text{Div},\Omega^{\text{F}},L^{2}(\partial\Omega^{\text{F}})$. By the Galerkin method, the finite-element solutions, $u_{h}$, and the test functions, $v_{h}$, both lie in $E_{h}$ and $\mathbf{V}_{h}^{p}$. We note that the polynomial degree $p^{p}$ does not need to be equal to $p^{f}$.

We apply non-conforming finite elements across the fluid-solid boundaries. The fluid-solid transmission condition in the definition of $E$ has been replaced by the condition $\int_{E^{{\text{FS}}}}[u_{h}\cdot\nu]_{-}^{+}v^{p}_{h}\,\mathrm{d}\Sigma=0$ in the definition of $E_{h}$. The fluid-solid transmission condition holds in the form of a boundary integration. For low-degree polyomials we show, in the next subsection, that these conditions are compatible through our formulation. Such a compatibility was analyzed and discussed by bermudez1994finite ; bermudez1995finite ; brezzi2012mixed . Several numerical studies kiefling1976fluid ; zienkiewicz1978fluid ; olson1985analysis ; chen1990vibration ; bermudez1999finite have been performed using similar non-conforming schemes along the fluid-solid boundaries. For the general theory and analysis of the mixed FEM, we refer to brezzi2012mixed .

4.2 Matrix formulae

We introduce nodal-based Lagrange polynomials, $\{\ell^{s}_{i}\}$, $\{\ell^{f}_{i}\}$, $\{\ell^{p}_{i}\}$, on the respective volume elements $K\in\mathcal{T}^{\text{S}}_{h}$, $\mathcal{T}^{\text{F}}_{h}$. We set $N_{p^{s}}=(p^{s}+1)(p^{s}+2)(p^{s}+3)/6$, where $N_{p^{s}}$ is the number of nodes on a tetrahedron for the $p^{s}$-th order polynomial approximation. We have similar expressions for $N_{p^{f}}$ and $N_{p^{p}}$. We write

[TABLE]

for $x\in K$; similar representations hold for $v^{s}_{h}$, $v^{f}_{h}$, $v^{p}_{h}$, respectively. We collect the values of $u^{s}_{h}$, $u^{f}_{h}$, $p_{h}$ and $v^{s}_{h}$, $v^{f}_{h}$, $v^{p}_{h}$ at all the nodes, $\{x_{i}\}$, in the vectors $\tilde{u}^{s}$, $\tilde{u}^{f}$, $\tilde{p}$ and $\tilde{v}^{s}$, $\tilde{v}^{f}$, $\tilde{v}^{p}$, respectively. We can then construct the corresponding submatrices, $A_{sg}$, $A_{f}$, $A_{p}$, $A_{\text{dg}}$, $A_{\text{dg}}^{\mathsf{T}}$, $E_{\text{FS}}$, $E_{\text{FS}}^{\mathsf{T}}$, $R_{s}$, $R_{f}$, $M_{s}$ and $M_{f}$, see Table 3, in a standard way summarized in Appendix A.

5 Self-gravitation as an N-body problem

Self-gravitation can be treated as the solution of an N-body problem. We discretize the entire planet into many elements and consider them as individual bodies. The gravitational potential and field are then computed through the interaction between these bodies. We note that FMM is an ideal candidate for solving an N-body problem. FMM reduces the complexity of the N-body problem from $O(N^{2})$ to $O(N\log N)$ or even $O(N)$ greengard1987fast . We apply the FMM greengard1997new ; gimbutas2011fmmlib3d to calculate the reference gravitational potential in Subsection 6.1. We employ ExaFMM yokota2013fmm , a massively parallel N-body problem solver, to solve for the perturbation of the gravitational potential.

5.1 Reference gravitational potential and gravitational field

For calculating the reference gravitational potential and field, we need to evaluate two integrals (dahlen1998theoretical, , (3.2) and (3.3)). The N-body problem of gravitation requires the evaluation of

[TABLE]

for the potential in (1) and

[TABLE]

for the field in (5). Here, $\mathbf{x}_{i}$ denotes the location of the target vertex and $\mathbf{r}_{k}$ denotes the barycenter of element $K_{k}$.

5.2 Incremental gravitational potential

For calculating the incremental gravitational potential, we need to evaluate (23) containing both the volume and boundary integral terms. Given the finite-element partitioning, $\mathcal{T}_{h}$, we approximate $S(u_{h})$ in (2) via

[TABLE]

and

[TABLE]

where $k_{1}$ and $k_{2}$ label the elements $K_{k_{1}}$ and $K_{k_{2}}$, $S_{k_{2}}(u_{h})$ is the incremental gravitational potential $S(u_{h})$ at the barycenter of $K_{k_{2}}$, $l_{1}$ and $l_{2}$ label the triangular elements $E_{l_{1}}$ and $E_{l_{2}}$, $\mathbf{r}_{l_{1}}$ and $\mathbf{r}_{l_{2}}$ denote the barycenters of $E_{l_{1}}$ and $E_{l_{2}}$. The first terms in (45) and (46) indicate the self-contribution.

Since the variation of $\nabla\cdot(\rho_{k_{2}}^{0}(x)u_{h}(x))$ is small on element $K_{k_{2}}$, we simplify the first term in (45) according to

[TABLE]

where $|K_{k_{2}}|$ denotes the volume of element $K_{k_{2}}$. We let

[TABLE]

and obtain

[TABLE]

Similarly, we simplify the first term in (46) according to

[TABLE]

with

[TABLE]

where $|E_{l_{2}}|$ denotes the area of the boundary element $E_{l_{2}}$. Note that $R_{k_{2}}$ in (47) and $R_{l_{2}}$ in (48) can be precomputed on each element and surface. The second and third terms in (45) and (46) may be evaluated via FMM.

5.2.1 Solid planets

For solid planets, we substitute (47) and (48) into (45) and (46), respectively. To evaluate (22) for a solid planet, we need to compute

[TABLE]

We add (49) into the matrix representation and obtain

[TABLE]

where $C_{s}\tilde{u}^{s}$ evaluates $S_{k_{2}}(u^{s}_{h})$ and $S_{l_{2}}(u^{s}_{h})$, $S_{s}$ solves the N-body problem for the solid planet, and $C_{s}^{\mathsf{T}}S_{s}C_{s}\tilde{u}^{s}$ evaluates (49); the submatrix $A_{sg}$ and its corresponding weak formula is shown in Table 3, and the submatrices $C_{s},C_{s}^{\mathsf{T}},S$ and their corresponding weak formulae are shown in Table 4. Here, of course, $A_{sg}$, $C_{s}$ and $C_{s}^{\mathsf{T}}$ do not include terms related the fluid-solid boundaries $\Sigma^{\text{FS}}$.

5.2.2 Planets with fluid regions

For a planet with fluid regions, we also substitute (47) and (48) into (45) and (46), respectively. To ensure the Hermitian property of the system, we carefully treat the fluid-solid boundary terms and evaluate the incremental gravitational potential $S(u_{h})$ via (24) and obtain the volume integral contributions

[TABLE]

and boundary integral contributions

[TABLE]

With (51) and (52), we have the full solution for the incremental gravitational potential. To evaluate (22) for a planet with fluid regions, we need to compute

[TABLE]

We derive the matrix representation with (53) and obtain

[TABLE]

with

[TABLE]

where $C\tilde{u}=C_{s}\tilde{u}^{s}+C_{f}\tilde{u}^{f}$ evaluates (51) and (52) to get $S_{k_{2}}(u_{h})$ and $S_{l_{2}}(u_{h})$, $S$ solves the N-body problem, and $C^{\mathsf{T}}SC\tilde{u}$ evaluates (53); the submatrices $A_{sg}$, $A_{f}$, $A_{p}$, $R_{s}$, $R_{f}$, $M_{s}$, $M_{f}$, $E_{\text{FS}}$, $A_{\text{dg}}$ and their corresponding weak formulae are shown in Table 3 and the submatrices $C_{s}$, $C_{s}^{\mathsf{T}}$, $C_{f}$, $C_{f}^{\mathsf{T}}$, $S$ and their corresponding weak formulae are shown in Table 4. The construction of submatrices $C_{s}$, $C_{s}^{\mathsf{T}}$, $C_{f}$, $C_{f}^{\mathsf{T}}$ can be found in A.4. We note that $A_{p}$ is always symmetric positive definite since $\kappa$ is always positive. We note that (54) is the discretization of (27).

6 Computational experiments for non-rotating planets

In this section, we first show the computational accuracy of our algorithm for the reference gravitational field using FMM in Subsection 6.1. We then illustrate computational experiments yielding planetary normal modes with or without perturbation of the gravitational potential. In this section and Section 7, two supercomputers, Stampede2 (an Intel cluster) at the Texas Advanced Computing Center and Abel (a Cray XC30 cluster) at Petroleum Geo-Services are utilized for the computational experiments.

6.1 Computational accuracy for the reference gravitational field

In this subsection, we illustrate the computational accuracy for the reference gravitational field using FMM. We begin with a simple constant-density ball. In Table 5, we show the FMM solution for a gravitational field of a constant density ball and a comparison with the closed-form solution. We note that FMM provides an accurate solution for this example.

We use PREM to build our Earth models on unstructured meshes with different sizes. In Table 6, we show the approximation errors of different three-layer models, which contain two major discontinuities (CMB and ICB) when compared with the semi-analytical solution. In Fig. 7, we show the comparison of the gravitational field computed via FMM with the semi-analytical solution in PREM.

In Table 7, we show the approximation errors of different seven-layer models which contain six major discontinuities (Moho, top of Low Velocity Zone (LVZ), bottom of LVZ, 660, CMB and ICB) with the semi-analytical solution.

6.2 Computational accuracy for non-rotating planets

In this subsection, we do not consider rotation and study the computational accuracy with existing algorithms for spherically-symmetric planets. Let the angular velocity of rotation $\Omega=0$, without loss of generality, we write (54) and its pure solid planet version (50) in the form of generalized eigenvalue problems:

[TABLE]

where $A$ represents $A_{sg}-C_{s}^{\mathsf{T}}S_{s}C_{s}$ in (50) or $A_{G}-E_{G}A_{p}^{-1}E_{G}^{\mathsf{T}}-C^{\mathsf{T}}SC$ in (54) and $\omega_{N}$ denotes the frequency for the non-rotation planets. Since the explicit formation of $A$ with self-gravitation requires excessive storage, it is necessary to solve (55) via a matrix-free scheme, where $A$, $M$ and $M^{-1}$ are only accessed through matrix-vector multiplications. We combine several efficient parallel approaches to solve (55) with a matrix-free scheme.

In this work, we utilize polynomial filtering techniques saad:filtered-cr06 ; Filtlan-paper ; spectrumslicing as these do not involve solving linear systems with the indefinite matrices. Here, the bulk of computations are carried out in the form of matrix-vector products. The polynomial filtering technique is ideally suited for solving large-scale three-dimentional interior eigenvalue problems because it significantly enhances the memory and computational efficiency without any loss of accuracy DBLP:conf/sc/ShiLXSH18 . In this paper, we adopt the polynomial filtering algorithms recently developed in spectrumslicing ; DBLP:conf/sc/ShiLXSH18 ; li2019evsl due to their simplicity and robustness on a prescribed interval $[f_{1},f_{2}]$ mHz. The details about our parallel algorithms and their performance can be found in DBLP:conf/sc/ShiLXSH18 .

We show the convergence of our numerical formulation and approach for constant elastic balls and PREM. The constant balls have a radius of 6,371 km, density $\rho^{0}=5.51\times 10^{3}$ kg/$m^{3}$, P-wave speed $V_{P}$ = 10.0 km/s and S-wave speed $V_{S}$ = 5.7735 km/s. The PREM used in our tests is modified in an isotropic model without attenuation, with $V_{P}=(V_{PV}+V_{PH})/2$ and $V_{S}=(V_{SV}+V_{SH})/2$. The ocean layer in PREM is replaced by crust. In the work of Matchette-Downes2021 , a good agreement of the one-dimensional solution based on the classical approach MINEOS masters2011mineos and a radial FEM jingchen2018revisiting is demonstrated. The discretization of the radial FEM code is described in Appendix B, In this work, we show our three-dimensional results are in a good agreement with the one-dimensional solutions.

6.2.1 Solid models with self-gravitation

We present our results for purely solid models with self-gravitation. In Tables 8 and 10, we list the number of elements ‘#elm.’ as well as the problem sizes (labeled as ‘size of $A$’ for the solid cases and ‘size of $A_{G}$’ and ‘size of $A_{p}$’ for the Earth examples), the number of surfaces ‘#surf.’, the size of $S_{s}$ or $S$, and the target frequency interval in milliHertz (labeled as $[f_{1},f_{2}]$ (mHz)), the degree of the polynomial filter ‘deg’, the number of the Lanczos iterations required ‘#it’, and the number of the normal modes computed ‘#eigs’.

Since the pure solid models do not generate any essential spectra, we can directly compute the lowest-frequency normal modes. We note that the length ($\lambda_{\max}-\lambda_{\min}$) of the spectrum grows with the size of the problem determined by the discretization.

In Table 9, we show the convergence results for different solid models using P1 elements, that is, the finite-element polynomial orders $p^{s}=p^{f}=p^{p}=1$ are used throughout this work. Through comparison with 1D results, we observe that our computational results do converge. We accept relative differences of about 0.1%.

In Table 10, we list test cases for different solid models using P2 elements, that is, the finite-element polynomial orders $p^{s}=p^{f}=p^{p}=2$ are used throughout this work. From experiments C1p2 to C5p2, we double the number of elements and obtain proper convergence results in Table 11. We show that even with about 330,000 elements, we are able to achieve four-digit agreement.

6.2.2 PREM with self-gravitation

Here, we include a liquid outer core using PREM and the presence of the essential spectrum. In Table 12, we show test cases for PREM. We roughly double the number of elements from E1p1 to E7p1. In Table 13, we argue convergence by comparing with 1D results. For PREM with self-gravitation, we accept relative differences that are less than 0.1%.

6.3 Fully heterogeneous models

Here, we study the effects of heterogeneity on the normal modes. In Subsection 6.3.2 and Subsection 6.3.1, we study the effects of the crust and upper mantle, and shape of the CMB, respectively.

6.3.1 Shape of the CMB

Here, we study the effects of the CMB. Long-wavelength topography of the CMB was proposed by creager1986aspherical ; morelli1987topography . Many studies bataille1988inhomogeneities ; doornbos1989models ; pulliam1993bumps ; rodgers1993inference ; obayashi1997p ; earle1997observations ; earle1998observations ; garcia2000amplitude ; sze2003core ; lassak2010core ; tanaka2010constraints ; colombi2014seismic ; schlaphorst2015investigation were later performed to model the topography of the CMB.

In Fig. 8, we show the topography of the CMB from the result by tanaka2010constraints . We use a triangular mesh to model the shape with ellipticity combined. In Table 14, we show the information of the experiment CMB8, which indicates a PREM-like model with the mentioned CMB embedded. In Fig. 9, we illustrate the splittings of modes ${}_{1}S_{7}$ and ${}_{1}S_{8}$ due to the non-spherically symmetric CMB. Since the modes ${}_{1}S_{7}$ and ${}_{1}S_{8}$ are sensitive to the change of the CMB, the splittings of these modes are quite clear.

6.3.2 Heterogeneity of the crust and upper mantle

Self-gravitation is important for the normal modes with frequencies lower than 5.0 mHz or so kennett1998density . However, in this subsection, we restrict ourselves to models without perturbation of the gravitational potential for computational efficiency. We reduce the full generalized eigenvalue problem (55) into Cowling approximation

[TABLE]

where $\omega_{C}$ is the frequency for Cowling approximation.

In Table 15, we show three different Earth models using the Cowling approximation. We construct two three-dimensional Earth models using MIT’s mantle tomographic results burdick2017model and crust 1.0 laske2013update . The core model is based on PREM. The mantle seismic reference wave speeds are based on AK135 kennett1995constraints . One model is obtained by combining MIT’s mantle tomographic model and PREM for the core and density. The other one replaces PREM’s crust by crust 1.0, which is shown in Fig. 3. In the first three rows of Table 15, we show the information of three different tests for these three different Earth models. Since with similar degrees of freedom, the largest eigenvalue of the MIT model with the three-dimensional crust is much larger than these of the other two models, we expect that significant mode coupling and splitting occur deuss2001theoretical ; romanowicz2008computation ; beghein2008signal ; irving2009normal ; koelemeijer2012normal ; nader2015normal ; yang2015synthetic ; akbarashrafi2017exact ; al2018hamilton .

We visualize different modes. The normal modes computed in the two MIT models are non-degenerate. In Fig. 10, we compare different modes computed in the three models in the frequency range $[2.0,2.5]$ mHz. Since the background models have only slight differences, some of the eigenfrequencies are similar amongst PREM and the MIT models. We illustrate most of the modes computed in PREM. In Fig. 10 (a), we observe that, even at low frequencies, weak mode splitting occurs for surface wave modes, including ${}_{2}S_{8}$, ${}_{0}S_{13}$, ${}_{0}T_{14}$ and ${}_{1}T_{7}$. We also report that no coupled modes are observed in $[2.0,2.18]$ mHz. In Figs. 10 (b-d), we show the different modes in $[2.18,2.28]$, $[2.28,2.38]$ and $[2.38,2.48]$ mHz, respectively. The splitting of most surface wave modes becomes larger with increasing frequency. However, since modes like ${}_{1}S_{10}$ (strong at the core-mantle boundary) in Fig. 10 (a), ${}_{0}c_{4}$ (an inner core toroidal mode) and ${}_{3}S_{5}$ (an ICB Stoneley mode) in Fig. 10 (c), are not sensitive to the crust and upper mantle structure, no clear splitting is observed. We observe coupled modes in Figs. 10 (b-d) computed in the MIT model with the three-dimensional crust. The eigenfunction of one mode in Fig. 10 (b) shows that ${}_{0}S_{14}$ and ${}_{2}T_{2}$ are coupled. The ${}_{0}T_{15}$ and ${}_{8}S_{1}$ near ${}_{0}S_{14}$ and ${}_{2}T_{2}$ are isolated multiplets. The eigenfunctions of the two modes in Fig. 10 (c) show that ${}_{1}S_{11}$ and ${}_{0}T_{16}$ are coupled. The ${}_{0}S_{15}$ near ${}_{1}S_{11}$ and ${}_{0}T_{16}$ is an isolated multiplet. These coupled modes are interesting because ${}_{1}S_{11}$ is clearly sensitive to the core-mantle boundary and the fundamental Love mode ${}_{0}T_{16}$ illustrated can be measured at the surface. The left mode in Fig. 10 (d) is a ${}_{0}S_{16}$ and ${}_{1}T_{9}$ coupled mode. The right mode in Fig. 10 (d) is a ${}_{6}S_{2}$ and ${}_{0}T_{17}$ coupled mode. This mode is also very interesting because ${}_{6}S_{2}$ illustrated is an inner core mode and the fundamental Love mode ${}_{0}T_{17}$ illustrated can be detected at the surface. Since the relative wave speed variations of the MIT tomographic model vary roughly from -1.4% to 1.4% in the upper mantle and the crust’s thickness is small, strong mode coupling occurs only to two modes. In this frequency range $[2.0,2.5]$ mHz, the width of each multiplet is small and no significant coupling between three or more modes is observed.

7 Computational experiments for rotating planets

In this section, we include the rotation and study its effects on normal modes. To simplify (54) and (50) without any loss of the generality, we extend (55) and derive a standard form for the QEP,

[TABLE]

We note that $\tilde{R}_{\Omega}=-\tilde{R}^{\mathsf{T}}_{\Omega}$, that is, $2\operatorname{i}\tilde{R}_{\Omega}$ is Hermitian. The eigenfrequencies are real and come in pairs $(\omega,-\omega)$.

To solve the QEP of the original form, the QEP is often projected onto a properly chosen low-dimensional subspace to facilitate the reduction to a QEP directly of lower dimension, such as in the Jacobi–Davidson method sleijpen1996jacobi ; sleijpen1996quadratic . The reduced QEP can then be solved by a standard dense matrix technique. Both Arnoldi- and Lanczos-type processes hoffnung2006krylov have been developed to build such projections of the QEP. A subspace approximation method holz2004subspace was presented via applying perturbation theory to the QEP. A second-order Arnoldi procedure bai2005soar was developed to generate an orthonormal basis for solving a large-scale QEP directly. We note that the above mentioned methods typically utilize a shift-and-invert scheme for solving the interior eigenpairs. These techniques become impractical for eigenvalue problems of the size of ours due to the high memory costs.

Instead, we can utilize extended Lanczos vectors from solving the generalized eigenvalue problem (55) through the polynomial filtering method. We then approximate the solution $\tilde{u}$ using the basis computed from

[TABLE]

where $X_{e}$ stands for the Ritz vectors of the linear system and $\Lambda_{e}$ denotes a diagonal matrix whose diagonal is a collection of $\omega_{N}^{2}$ in (55). We take $m_{e}$ eigenvectors spanning a subspace and let $\tilde{u}_{e}=X_{e}y_{e}$ approximate $\tilde{u}$ in (57), where $y_{e}$ is complex. We apply

[TABLE]

to an equivalent form of (57),

[TABLE]

Making use of $X_{e}^{\mathsf{T}}AX_{e}=\Lambda_{e}$, we obtain

[TABLE]

It is apparent that if $\tilde{R}_{\Omega}=0$, we have $\omega_{e}=\omega_{N}=\Lambda_{e}^{1/2}$. The system (59) can be solved with a standard eigensolver such as the one implemented in LAPACK anderson1999lapack . Here, we study the spectra of two models: Earth 1066A gilbert1975application and a Mars model khan2016single . We use 23.9345 hours allen1973astrophysical and 24.6229 hours lodders1998planetary as Earth’s and Mars’ rotation periods, respectively. With a large $m_{e}$ and a relatively small $\Omega$, the numerical solution $\omega_{e}$ is close to $\omega$ in (57). The numerical accuracy can further be improved via solving (57) exactly.

7.1 Computational accuracy

For small models, we are able to compute the full mode expansion associated with the point spectrum using (59). In Table 16, we list the numerical parameter values pertaining to the testing of computational accuracy and estimating the cost in different models: The number of elements (labeled as # of elm.), size of $A_{p}$, size of $A$, size of $S$ and the target frequency interval in milliHertz (labeled as $[f_{1},f_{2}]$ (mHz)).

In Figs. 11 (a)–(c), we illustrate the computational accuracy of tests in three different models, C3kp1, E3kp1 and M2kp1, respectively, on the lowest seismic eigenfrequencies using P1 elements. We compare the differences in the eigenfrequencies between the full mode expansion and a 200 mode expansion. The differences are about $5\times 10^{-6}$ mHz, which is two digits below the accuracy of common normal mode measurements.

In Figs. 12 (a) and (b), we show the computational accuracy of M8kp1 on [0.4, 1.14] mHz as well as the error distribution. In Fig. 12 (a), we show that even with a 100 mode expansion, the differences are as low as $1\times 10^{-5}$ mHz. In Fig. 12 (b), we show that with a 1000 mode expansion, the differences are further reduced to about $1\times 10^{-6}$ mHz.

7.2 Benchmark experiments for Earth models with rotation

Over the past two decades, a significant number of observational studies have been carried out to the rotation effects on the Earth’s normal modes zurn2000observation ; millot2003normal ; park2005earth ; roult2010observation ; nader2015normal ; schimmel2018low . Our computational approach can aid and complement such studies through accurate and consistent simulations generating even relatively high eigenfrequencies. Here, we perform a benchmark experiment of Earth model 1066A gilbert1975application against a perturbation calculation dahlen1979rotational . In the perturbation calculation, the eigenfrequency perturbations $\delta\omega_{m}$ have the following form

[TABLE]

where $\omega_{0}$ denotes the eigenfrequency of the unperturbed spherically symmetric model, $l$ denotes the angular order in the spherical harmonic expansion, and $a$, $b$ and $c$ are the relevant coefficients. The values of $a$, $b$ and $c$ for different radial modes can be found in (dahlen1998theoretical, , Table 14.1). In Table 17, we list the numerical parameters of the Earth models in the benchmark test. The models E1Mp1 and E2Mp2 used to compute $\omega_{0}$ represent spherically symmetric ones without rotation. Experiments EE1Mp1 and EE2Mp2 represent elliptic Earth models and are used to compute eigenfrequencies with our proposed method. The ellipticities of the Earth models are computed by solving Clairaut’s equation (cf. Section 3). Since the eigenfrequencies of the Slichter modes slichter1961fundamental are close to the upper bound of the essential spectrum and the convergence of the proposed algorithm is relatively slow, we set $f_{1}=$ 0.04 mHz and use experiments E1Mp1 and EE1Mp1 to compute the Slichter modes using P1 elements. Experiments E2Mp2 and EE2Mp2 are used to compute other modes using P2 elements. In Fig. 13, we show the comparison between the perturbation and our methods. The values of the computed eigenfrequenies of our method agree with the perturbation results in as much as that the relative differences are commonly less than 0.3 $\mu$Hz. The degree of agreement is, of course, model dependent. The eccentricity in the Earth model is so small that the second-order perturbation is accurate within the typical error of our numerical computations. Higher rotation rates would increase the eccentricity and let the second-order perturbation loose accuracy.

7.3 Mars models

Here, we present our computational results for Mars models. The interiors of the Mars models are based on mineral physics calculations khan2016single . In Table 18, we list three Mars models labeled as M2Mp2, EM2Mp2 and TM2Mp2 which represent a spherically symmetric Mars model without rotation, a spheroidal Mars model with rotation, and a spheroidal Mars model with a three-dimensional crust and rotation using P2 elements. The shape of the spheroidal Mars model’s core-mantle boundary is computed by solving Clairaut’s equation. Since Mars presumably is not hydrostatic as discussed in Section 3, its solid region is estimated via a linear interpolation using the ellipticities of the core-mantle boundary ($\varepsilon=4.19\times 10^{-3}$) and the surface ($\varepsilon=5.89\times 10^{-3}$). Model TM2Mp2 is illustrated in Fig. 5.

In Fig. 14, we show eigenfrequencies computed in different Mars models listed in Table 18. Symbols $\bullet$, $\circ$ and $\times$ represent the eigenfrequencies computed in Mars models M2Mp2, EM2Mp2 and TM2Mp2 (cf. Table 18). The horizontal dashed lines represent the eigenfrequencies of a spherically symmetric Mars model computed with a one-dimensional solver masters2011mineos ; jingchen2018revisiting . Mode splitting is apparent due to ellipticity, rotation and heterogeneity in three dimensions. The three-dimensional crust does not have a clear influence on the lowest eigenfrequencies associated with ${}_{0}S_{2}$, ${}_{0}T_{2}$, ${}_{1}S_{1}$, ${}_{0}S_{3}$, ${}_{0}T_{3}$, ${}_{1}S_{2}$ and ${}_{0}S_{4}$ in Fig. 14 (a). The three-dimensional crust has a noticeable effect on the surface wave modes, such as ${}_{0}T_{6}$, ${}_{0}T_{7}$, ${}_{0}T_{8}$, ${}_{0}S_{6}$, ${}_{0}S_{7}$ and ${}_{0}S_{8}$, as expected. In Fig. 15, we show the eigenfrequencies in a subinterval of the interval used in Fig. 14 (d). Here, we note the splitting of modes ${}_{2}S_{4}$, ${}_{0}S_{8}$ and ${}_{0}T_{8}$ and highlight the effects of the three-dimensional crust. The maximum difference between the eigenfrequencies in Fig. 15 is 5.2 $\mu$Hz, which, in principle, can be detected. There is no mode-coupling observed in these experiments.

In Fig. 16, we plot the branch ${}_{1}S_{l}$ as well as the corresponding incremental gravitational fields $\nabla S(u)$. The superconducting gravimeters are expected to contribute to normal mode seismology crossley1999network ; van1999measuring ; widmer2003can ; rosat2003search ; hafner2012signature . We anticipate that both the seismic and gravity measurements of these modes could help to estimate the size of the Martian core.

8 Conclusion

In this work, we propose a method to compute the normal modes of a fully heterogeneous rotating planet. We apply the mixed finite-element method to the elastic-gravitational system of a rotating planet and utilize the FMM to calculate the self-gravitation. We successfully separate out the essential spectrum by using a polynomial filtering eigensolver and thus, are able to compute the normal modes associated with seismic point spectrum. To solve the relevant QEP, we utilize extended Lanczos vectors computed in a non-rotating planet – with the shape of boundaries of a rotating planet and accounting for the centrifugal potential – spanning a subspace to reduce the dimension of an equivalent linear form of the QEP. The reduced system can be solved with a standard eigensolver. We demonstrate our ability to compute the seismic normal modes with and without rotation accurately. We then study the computational accuracy and use a standard Earth model to perform a benchmark test against a perturbation calculation. We carry out computational experiments on various Mars models and illustrate mode splitting due to rotation, ellipticity and heterogeneity of the crust. The use of modern supercomputers enables us to capture normal modes associated with the seismic point spectrum of a fully heterogeneous planet accurately. The computational efficiency can be further improved by using acceleration techniques. The extension to include viscoelastic relaxation (for a review, see romanowicz20071 ), in particular Maxwell and Burger models, leads to a nonlinear rational eigenvalue problem, which is tractable at current subject of research.

Acknowledgement

We would like to thank Bernard Valette for his thoughtful comments. J.S. would like to thank Petroleum Geo-Services for using their supercomputer Abel, and Danny Sorensen, Ruichao Ye, and Harry Matchette-Downes for helpful discussions.

Appendix A Construction of orthonormal bases and submatrices

Here, we introduce three-dimensional polynomial bases $\{\psi_{n}^{s}\}_{n=1}^{N_{p^{s}}}$, $\{\psi_{n}^{f}\}_{n=1}^{N_{p^{f}}}$ and $\{\psi_{n}^{p}\}_{n=1}^{N_{p^{p}}}$ while addressing the fact that the Lagrange polynomials are not orthogonal to one another. We suppress superscripts $s$, $f$, $p$ in the notation in the remainder of this subsection. To simplify the computations, we introduce reference volume and boundary elements. That is, we introduce a mapping that connects any element $K$ to the reference tetrahedron defined by

[TABLE]

Likewise, we introduce a mapping that connects any boundary element $E$ to the reference triangle defined by

[TABLE]

We note that any two tetrahedra are connected through an affine transformation, $x\to r$, with a constant Jacobian, $J$, which is the determinant of $(\partial_{r}x)$. For the local approximation on the reference element $\mathbf{I}$, we have

[TABLE]

The vector fields are treated component-wise in our discretization. This yields the expression $\mathcal{V}\hat{u}_{j}=u_{j}$, where the generalized Vandermonde matrix takes the form of $\mathcal{V}_{in}=\psi_{n}(r_{i})$ with $i,n$ as indices of nodal points. Here, $\{\psi_{n}\}$ is a polynomial basis that is orthonormal on $\mathbf{I}$. We later introduce submatrices of $\mathcal{V}$. We then evaluate derivatives and mass matrices according to

[TABLE]

where $\mathcal{D}_{j}$ and $\mathcal{M}$ are the derivative matrix and the mass matrix on the reference tetrahedron. More details of the constructions of $J$, $\mathcal{V}$, $\mathcal{D}_{j}$ and $\mathcal{M}$ can be found in (hesthaven2007nodal, , Chapter 10.1). Thus, we introduce

[TABLE]

We employ the notation

[TABLE]

reflecting the mapping of the derivatives from the reference tetrahedron to the target element. We follow a similar approach for boundary elements and introduce

[TABLE]

where $\mathcal{M}_{s}^{2D}$ and $\mathcal{M}_{f}^{2D}$ are the mass matrices for solid and fluid boundary elements, respectively; $J^{2D}$ denotes the Jacobian, which is the determinant of $(\partial_{t}x)$ on the boundary element. The construction of the mass matrices $\mathcal{M}^{2D}_{s}$ and $\mathcal{M}^{2D}_{f}$ on the reference triangle $\mathbf{I}_{2D}$ is similar to the construction of $\mathcal{M}$ (hesthaven2007nodal, , Chapter 6.1).

A.1 Submatrices: $A_{sg}$, $A_{f}$, $A_{p}$, $M_{s}$, $M_{f}$, $R_{s}$ and $R_{f}$

We extract $\tilde{u}^{s}|_{K_{k}}$, $\tilde{u}^{f}|_{K_{k}}$ and $\tilde{p}|_{K_{k}}$ from $\tilde{u}^{s}$, $\tilde{u}^{f}$ and $\tilde{p}$, respectively, by restricting the nodes to the ones of element $K_{k}$. In a similar fashion, we extract $\tilde{v}^{s}|_{K_{k}}$, $\tilde{v}^{f}|_{K_{k}}$ and $\tilde{v}^{p}|_{K_{k}}$ on any element $K_{k}$. For the evaluation of matrix $A_{sg}$ in Table 3 we need to evaluate the submatrices on element $K_{k}$ through

[TABLE]

where $c^{k}_{ijmn}$, $\rho^{0}_{k}$ and $J_{k}$ denote the stiffness tensor, density and the Jacobian on element $K_{k}$, respectively; $D_{g^{\prime}_{i}}$ and $D_{\partial_{x_{i}}g^{\prime}_{j}}$ denote the diagonal matrices whose diagonal entries are $g^{\prime}_{i}$ and $\partial_{x_{i}}g^{\prime}_{j}$, respectively. For the evaluation of the boundary integration in $A_{sg}$, we need to evaluate the submatrix on element $E_{l}^{\text{FS}}$ through

[TABLE]

where $\rho^{0}_{l}$ and $\nu^{s\rightarrow f}_{j}|_{E_{l}}$ denote the density and normal vector on the boundary element $E_{l}^{\text{FS}}$, respectively, upon extracting $\tilde{v}^{s}_{i}|_{E_{l}}$ and $\tilde{u}^{s}_{i}|_{E_{l}}$. We can deal with the integral over $\Sigma^{\text{FF}}$ similarly.

We then evaluate the submatrices for $A_{f}$, $A_{p}$, $M_{s}$, $M_{f}$ in Table 3 and obtain

[TABLE]

where $D_{g^{\prime}_{j}/\|g^{\prime}\|}$ denotes a diagonal matrix whose diagonal entries are $g^{\prime}_{j}/\|g^{\prime}\|$ and $N_{k}^{2}$ denotes the square of the Brunt-Väisälä frequency on element $K_{k}$. We also obtain the rotation components $R_{s}$ and $R_{f}$,

[TABLE]

where $\epsilon_{ilm}$ denotes the Levi-Civita symbol.

A.2 Submatrices: $A_{\text{dg}}$ and $A_{\text{dg}}^{\mathsf{T}}$

Here, we discuss the integration between the different variables. For the inner products between $u^{f}_{h}$ and $p_{h}$ for $A_{\text{dg}}$ and $A_{\text{dg}}^{\mathsf{T}}$ in Table 3, we evaluate the mass matrices $\mathcal{M}_{pf}$ and $\mathcal{M}_{fp}$,

[TABLE]

where we refine the notation to indicate submatrices of $\mathcal{V}$; $\mathcal{V}(I)$ denotes the submatrix of $\mathcal{V}$ formed by columns indexed by $I\subseteq\{1,\ldots,N_{p}\}$. The selection of submatrices is based on the polynomial construction (hesthaven2007nodal, , (10.6)). For instance, if the polynomial orders used for both $u^{f}_{h}$ and $p_{h}$ are the same, i.e., $p^{f}=p^{p}$, $I_{f}=I_{p}=\{1,\ldots,N_{p^{f}}\}$; if $p^{p}=1$ and $p^{f}=2$, we have $N_{p^{p}}=4$, $N_{p^{f}}=10$ and $I_{f}=\{1,2,3,4\}$, $I_{p}=\{1,2,4,7\}$. It is apparent that $\mathcal{M}_{pf}=\mathcal{M}_{fp}^{T}$.

Evaluating $A_{\text{dg}}$ in Table 3 requires the evaluation of the submatrices on element $K_{k}$ through

[TABLE]

where $\kappa^{-1}_{k}$ denotes the inverse of the bulk modulus on element $K_{k}$. To evaluate $A_{\text{dg}}^{\mathsf{T}}$ in Table 3, we also need to evaluate the submatrices on element $K_{k}$ through

[TABLE]

A.3 Submatrices: $E_{\text{FS}}$ and $E_{\text{FS}}^{\mathsf{T}}$

For $E_{\text{FS}}$ and $E_{\text{FS}}^{\mathsf{T}}$ , similar to Section A.2, we introduce two new indices to construct $\mathcal{M}_{ps}^{2D}$ and $\mathcal{M}_{sp}^{2D}$ on the boundary elements associated with the fluid-solid boundary. The selection of the submatrix is based on (hesthaven2007nodal, , Chapter 6). $\mathcal{M}_{ps}^{2D}={\mathcal{M}_{sp}^{2D}}^{\mathsf{T}}$ holds true as well. To evaluate $E_{\text{FS}}^{\mathsf{T}}$ in Table 3, we need to compute the submatrix on boundary element $E^{\text{FS}}_{l}$ through

[TABLE]

upon extracting $\tilde{p}|_{E_{l}}$ on boundary element $E^{\text{FS}}_{l}$. To evaluate $E_{\text{FS}}$ in Table 3, we need to evaluate the submatrix on boundary element $E^{\text{FS}}_{l}$ through

[TABLE]

upon extracting $\tilde{v}^{p}|_{E_{l}}$ on $E^{\text{FS}}_{l}$.

We are now able to build all the submatrices for the evaluation of the integrals in Table 3. We then assemble the global matrices from all these submatrices using standard techniques similar to those in bathe2006finite ; hughes2012finite .

A.4 Construction of the submatrices for the perturbation of the gravitational potential

Similar to the previous subsections, we construct the submatrices in $C_{s}$ in Table 4,

[TABLE]

and the submatrices in $C_{s}^{\mathsf{T}}$,

[TABLE]

where $\mathbf{1}$ denotes a vector of all ones. The construction of the submatrices in $C_{f}$ and $C_{f}^{\mathsf{T}}$ is the same. We are now able to build all the submatrices for the evaluation of the integrals in Table 4.

Appendix B Full mode coupling

Concerning the Galerkin approximation, we can use different, nonlocal bases of functions in the appropriate energy space, for example, the spectral-Galerkin method shen1994efficient . In this appendix, we consider the use of the eigenfunctions of a spherically symmetric, non-rotating, perfectly elastic and isotropic (SNREI) reference model as a basis in this method. This has been implemented by woodhouse1978effect ; woodhouse1980coupling ; deuss2001theoretical ; deuss2004iteration , and named the full mode coupling approach. An immediate drawback of using this basis, however, is that the fluid-solid boundaries need to be spherically symmetric, as these are encoded in these basis functions.

We let $u_{km}$ represent the eigenfunctions associated with eigenfrequencies, $\omega_{k}$, in terms of spherical harmonics, $Y_{l}^{m}$, that is,

[TABLE]

where $k$ is the multi-index for the eigenfrequency; $m=-l,-l+1\ldots,l-1,l$ is the index corresponding with the degeneracy with $l$ denoting the spherical harmonic degree; $U_{km},V_{km}$ and $W_{km}$ are the three components of eigenfunctions and are functions of the radial coordinate; $\mathbf{P}_{lm}$, $\mathbf{B}_{lm}$ and $\mathbf{C}_{lm}$ are the vector spherical harmonics, see (dahlen1998theoretical, , (8.36)) for their definition. In addition, $p_{km}$ needs to be introduced to constrain the solution, cf. (13) (de2019note, , Subsection 3.3). Since $\nabla\cdot u_{km}(x)$ can be expanded using $Y_{l}^{m}(x)$ (dahlen1998theoretical, , (8.38)) and $u_{km}(x)\cdot g_{(r)}$ can also be expanded using $Y_{l}^{m}(x)$ for the radial models, we let $p_{km}=P_{km}Y_{l}^{m}$ with

[TABLE]

where $\rho_{(r)}^{0}$, $\kappa_{(r)}$ and $g_{(r)}$ denote the radial profiles of the density, bulk modulus and reference gravitational field of a radial model, respectively. Similarly, the incremental gravitational potential of the radial models takes the form, $s_{km}=S_{km}Y_{l}^{m}$, where $S_{km}$ is also a function in the radial coordinate. In the following, $l$ and $m$ are fixed.

In a SNREI model, for the computation of the toroidal modes, we only need to consider a solid annulus comprising the mantle and the crust. We exemplify the computations with the spheroidal modes and let $U^{\prime}_{km}$, $P^{\prime}_{km}$ and $S^{\prime}_{km}$ be test functions for $U_{km}$, $P_{km}$ and $S_{km}$ following the Galerkin method. We let the $\tilde{X}_{(r)}$ be the 1D interval of the radial planet and have $\tilde{X}_{(r)}=\Omega_{(r)}^{\text{S}}\cup\Omega_{(r)}^{\text{F}}$, where $\Omega_{(r)}^{\text{S}}$ and $\Omega_{(r)}^{\text{F}}$ denote the 1D intervals for the solid and fluid regions, respectively. Given a regular finite-element partitioning $\mathcal{T}_{h}^{(r)}$ of the interval $\tilde{X}_{(r)}$, we denote an element of the mesh by $L_{q}\in\mathcal{T}_{h}^{(r)}$ and have $\tilde{X}_{(r)}=\bigcup_{q=1}^{N_{L}}L_{q}$, where $N_{L}$ denotes the total number of 1D elements. Furthermore, we let $L_{q}^{\text{S}}$ and $L_{q}^{\text{F}}$ specifically be elements in the solid and fluid regions and have

[TABLE]

where $N_{L}^{\text{S}}$ and $N_{L}^{\text{F}}$ denote the numbers of 1D elements in the solid and fluid regions, respectively. We let $\Sigma^{\text{FS}}_{(r)}$ denote the fluid-solid boundary points in the radial interval. We introduce the finite-element solutions, $U_{km;h}^{s}$, $U_{km;h}^{f}$, $V_{km;h}^{s}$, $V_{km;h}^{f}$, $P_{km;h}$ and $S_{km;h}$, and test functions, $U_{km;h}^{s^{\prime}}$, $U_{km;h}^{f^{\prime}}$, $V_{km;h}^{s^{\prime}}$, $V_{km;h}^{f^{\prime}}$, $P^{\prime}_{km;h}$ and $S^{\prime}_{km;h}$. We set $N_{p^{U}}=(p^{U}+1)/2$, where $N_{p^{U}}$ is the number of nodes on a 1D element for the $p^{U}$-th order polynomial approximation. We have likewise expressions for $N_{p^{V}}$, $N_{p^{P}}$ and $N_{p^{S}}$. As in Subsection 4.2, we introduce nodal-based Lagrange polynomials, $\ell_{i}^{U}$, $\ell_{i}^{V}$, $\ell_{i}^{P}$, $\ell_{i}^{S}$, on the respective 1D elements $L\in\mathcal{T}_{h}^{(r)}$, and write

[TABLE]

for $x\in L^{\text{S}}$ and $x\in L^{\text{F}}$, respectively; similar representations hold for $U_{km;h}^{s^{\prime}}$, $U_{km;h}^{f^{\prime}}$, $V_{km;h}^{s^{\prime}}$, $V_{km;h}^{f^{\prime}}$, $P^{\prime}_{km;h}$ and $S^{\prime}_{km;h}$, respectively. We note that the fluid-solid boundary points coincide with nodes.

As in Subsection 4 and Section 5, we collect the “values” of $U_{km;h}^{s}$, $U_{km;h}^{f}$, $V_{km;h}^{s}$, $V_{km;h}^{f}$, $P_{km;h}$ and $S_{km;h}$ at all the nodes, in vectors $\tilde{U}_{km}^{s}$, $\tilde{U}_{km}^{f}$, $\tilde{V}_{km}^{s}$, $\tilde{V}_{km}^{f}$, $\tilde{V}_{km}$ and $\tilde{S}_{km}$, respectively, and collect the values of $U_{km;h}^{s^{\prime}}$, $U_{km;h}^{f^{\prime}}$, $V_{km;h}^{s^{\prime}}$, $V_{km;h}^{f^{\prime}}$, $P^{\prime}_{km;h}$ and $S^{\prime}_{km;h}$ at all the nodes, in “vectors” $\tilde{U}_{km}^{s^{\prime}}$, $\tilde{U}_{km}^{f^{\prime}}$, $\tilde{V}_{km}^{s^{\prime}}$, $\tilde{V}_{km}^{f^{\prime}}$, $\tilde{P}^{\prime}_{km}$ and $\tilde{S}^{\prime}_{km}$, respectively. We let

[TABLE]

and obtain the resulting eigenvalue problem (cf. (54))

[TABLE]

where

[TABLE]

in which $A_{sg}^{(r)}$, $A_{f}^{(r)}$, $A_{p}^{(r)}$, $E_{\text{FS}}^{(r)}$, ${E_{\text{FS}}^{(r)}}^{\mathsf{T}}$, $A_{\text{dg}}^{(r)}$, ${A_{\text{dg}}^{(r)}}^{\mathsf{T}}$, $M_{s}^{(r)}$, $M_{f}^{(r)}$, ${C_{s}^{(r)}}^{\mathsf{T}}$, ${C_{f}^{(r)}}^{\mathsf{T}}$, $S^{(r)}$, $C_{s}^{(r)}$ and $C_{f}^{(r)}$, are given in Tables 19 and 20. We note that the matrices in (87) are obtained using separation of variables with spherical harmonics in (54). We substitute

[TABLE]

upon solving (17) and

[TABLE]

upon solving (2). We only need to invoke a finite-element basis in the radial coordinate. We note that the resulting system can be solved via a standard eigensolver, such as LAPACK anderson1999lapack .

As mentioned above, we may consider the finite-element solution denoted as $\{u_{km;h}\}$ as an alternative basis. Since $\{u_{km;h}\}$ is a global basis for the general problem, we have no separation in the solid and fluid components and no longer have the fluid-solid boundary terms in the system. Following the Galerkin method, we then consider an expansion for the general solution $u_{c}=\sum_{km}y_{km}u_{km;h}$ and the corresponding test functions $v_{c}=\sum_{k^{\prime}m^{\prime}}y^{\prime}_{k^{\prime}m^{\prime}}u_{k^{\prime}m^{\prime};h}$. We introduce $s_{c}$ and its corresponding test functions $v^{s_{c}}$ for self-gravitation. We have $s_{c}=\sum_{km}z_{km}S_{km;h}$ and $v^{s_{c}}=\sum_{k^{\prime}m^{\prime}}z^{\prime}_{k^{\prime}m^{\prime}}S_{k^{\prime}m^{\prime};h}$. Assuming that all the discontinuities in a fully heterogeneous model coincide with the ones in the reference radial model and the fluid outer core, the eigenfuncions represented by the mentioned expansions lie in $H_{1}\subset E$ (cf. (32)) for the fully heterogeneous problem while the constraint equation disappears. We let $y$, $y^{\prime}$, $z$ and $z^{\prime}$ be the “vectors” with components $y_{km}$, $y^{\prime}_{k^{\prime}m^{\prime}}$, $z_{km}$ and $z^{\prime}_{k^{\prime}m^{\prime}}$, respectively, and obtain

[TABLE]

as the counterpart of (54). Here, $A_{G}^{(c)}$, $M^{(c)}$, ${C^{(c)}}^{\mathsf{T}}$, $S^{(c)}$ and $C^{(c)}$, obtained via substituting the above-mentioned expansion of $u_{c}$ in (54), are given in Tables 21 and 22.

If all the discontinuities in a fully heterogeneous model with a fixed fluid outer core coincide with the reference radial model, we note that the matrix elements in (88), Tables 21 and 22 are similar to (woodhouse1980coupling, , (A1)), which describe mode coupling in non-radial models. However, Woodhouse (woodhouse1980coupling, , (A1)) includes additional terms accounting for changes in the fluid-solid boundaries while in the previous work (woodhouse1978effect, , (42)), perturbation theory is used to compute the eigenfrequency changes in terms of the unperturbed eigenfunctions; both calculations violate the condition that normal modes need to remain in $E$ and in $H_{1}$.