# Sixth-order spherical quadrature rule based on the trapezoidal rule on the Cubed Sphere

Jean-Baptiste Bellet

► To cite this version:

Jean-Baptiste Bellet. Sixth-order spherical quadrature rule based on the trapezoidal rule on the Cubed Sphere. 2026. hal-05497324

HAL Id: hal-05497324

https://hal.science/hal-05497324v1

Preprint submitted on 6 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# SIXTH-ORDER SPHERICAL QUADRATURE RULE BASED ON THE TRAPEZOIDAL RULE ON THE CUBED SPHERE

JEAN-BAPTISTE BELLET†

ABSTRACT. The Cubed Sphere is a spherical grid based on cartesian grids on the faces of a cube, widely used since its advent in the early seventies. The derivation of quadrature rules on this grid is still an active subject. Among the simplest rules, a bivariate trapezoidal rule on each face of the cube can be used. This strategy converges to any order for smooth functions supported inside the faces of the cube; more generally, it is only second-order accurate. However, J.-P. Croisille introduced about fifteen years ago an ingenious modification of the weight on the eight corners, and observed the fourth-order accuracy; an accuracy which has not been theoretically proved so far. The first main result of this paper solves this open question. Next, a deeper analysis of the remainder permits to identify a further correction at the corners, yielding the second main result: a new quadrature rule which is sixth-order accurate. Numerical experiments are in agreement with this order of accuracy, and pave the way towards corrections of higher orders, based on values near the corners only.

## 1. INTRODUCTION

Cubed Spheres are spherical grids introduced more than fifty years ago [39]. They have been further studied, for instance in [6, 33, 35–37], and a wide variety of numerical methods have been developed on Cubed Spheres, e.g. in [4, 7, 10, 12, 15–18, 21, 25–27, 29–31, 34, 38, 41, 43, 45]. Among these methods, various quadrature rules on the equiangular Cubed Sphere have been introduced over the past fifteen years, e.g. in [9, 11, 14, 19, 20, 32]. Such rules enter in the framework of spherical quadrature rules, which is an old but still active topic [1–3, 5, 13, 22–24, 28, 40, 44].

On the equiangular Cubed Sphere, two main classes of quadrature rules have been considered. The first class is related to data approximation and the correct integration of suitable spherical harmonics (or polynomials) [9,11,32]. The associated rules are very efficient and accurate, but the computation of suitable weights is not an easy task in general, nor the proof of their positivity. In particular, it is difficult to consider very thin grids. The second class focuses rather on the order of convergence [14,20,32], in order to make easier the increase of the resolution. Up to a change of variable, the Cubed Sphere contains especially six square faces, on which usual quadrature methods with simple weights can be applied; for instance, one tensor product of trapezoidal rules per face.

This paper focuses on the accuracy of such a bivariate trapezoidal rule applied for the equiangular Cubed Sphere. The famous Euler-Maclaurin formula expresses the remainder term of the one-dimensional trapezoidal rule in term of the derivatives on the boundary; in particular, convergence to any order is achieved for any smooth periodic function [42]. In the case of the Cubed Sphere, we show that the remainder term depends especially on derivatives on the twelve edges and on cross-derivatives at the eight vertices. One question of interest concerns whether or not these terms could be compensated (analogously as periodicity compensates error terms in one dimension).

An ingenious correction of the trapezoidal rule has been introduced at the eight corners in [19, 20], further studied in [14,32]. The resulting method is fourth-order accurate in the numerical experiments, but proving theoretically this order is still an open question. This paper solves this open problem. Indeed, we analyze the remainder. We succeed in expressing the second-order term, and the fourth-order one, in term of the values at the corners. The second-order term indicates the correction to be applied to reach a fourth-order accurate rule; the corrected

Date: February 6, 2026.

Key words and phrases. Cubed sphere, spherical quadrature, trapezoidal rule.
rule coincides with the method originally introduced in [19, 20]. Next, the fourth-order term in the remainder permits to apply a further correction which is simple and new; it yields a novel quadrature rule, which is sixth-order accurate.

The paper is organized as follows. In Section 2, we set the notation of the Cubed Sphere. In Section 3, we state the main results of the paper. In Section 4, we prove these results. In Section 5, numerical experiments complete the theoretical part. In Section 6, we conclude the work with some perspectives. Appendix A presents a two-dimensional Euler-Maclaurin formula. Appendix B contains a Maple code exploring the fourth-order term. Appendix C contains a Matlab function implementing the proposed quadrature rule.

## 2. NOTATION OF THE CUBED SPHERE

In this section, we outline the background about the equiangular Cubed Sphere.

**2.1. Cubed Sphere.** The equiangular Cubed Sphere with parameter $N \ge 1$, denoted by $CS_N$, is a subset of the unit sphere $S^2 \subset \mathbb{R}^3$, defined by

$$CS_N = \{ \frac{1}{r}(\pm 1, u, v), \frac{1}{r}(u, \pm 1, v), \frac{1}{r}(u, v, \pm 1); \text{ with } r = (1 + u^2 + v^2)^{1/2}, \\ u = \tan(-\frac{\pi}{4} + ih), v = \tan(-\frac{\pi}{4} + jh), 0 \le i, j \le N \}; \quad h = \frac{\pi}{2N}. \quad (1)$$

As displayed in Figure 1, it turns out that

$$CS_N = x(C_N), \quad \text{with} \quad C_N = C \cap \{-\frac{\pi}{4} + ih, 0 \le i \le N\}^3, \quad (2)$$

where the mapping

$$x: \alpha \in C \mapsto x(\alpha) = \frac{(\tan \alpha_1, \tan \alpha_2, \tan \alpha_3)}{(\tan^2 \alpha_1 + \tan^2 \alpha_2 + \tan^2 \alpha_3)^{1/2}} \in S^2 \quad (3)$$

defines a one-to-one correspondence between the boundary

$$C = \partial[-\frac{\pi}{4}, \frac{\pi}{4}]^3 \quad (4)$$

of the cube $[-\frac{\pi}{4}, \frac{\pi}{4}]^3$ and the sphere $S^2$. In (2), the Cubed Sphere $CS_N$ is in correspondence with a set $C_N$, which contains a uniform cartesian grid with step $h$ on each face of $C$.

<figure><img src="image_2.png" /><figcaption>FIGURE 1. One-to-one correspondence $x$ in (3), between the cube $C = \partial[-\frac{\pi}{4}, \frac{\pi}{4}]^3$ and the sphere $S^2$. The Cubed Sphere $CS_N \subset S^2$ in (1) corresponds with the set $C_N \subset C$ in (2), based on uniform cartesian grids on the faces of $C$ (step $h = \frac{\pi}{2N}$). The subset $CS_1 \subset CS_N$ corresponds with the set $C_1 \subset C_N$ of the vertices of $C$, described in (5).</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Edges and vertices of the cube C in (4). The set E in (6) contains the twelve edges of C, the set $C_1$ in (5) the eight vertices. The unit vectors $\nu_1$, $\nu_2$ are normal to E at $\alpha \in E \setminus C_1$; they point towards the interior of the faces. The normal vector $w = \nu_1 + \nu_2$ points towards the interior of the cube.</figcaption></figure>

2.2. **Vertices and edges.** The subset $C_1 \subset C_N$ contains the vertices of the cube C, and corresponds with the subset $CS_1 \subset CS_N$; these subsets, displayed in Figure 1, satisfy

$$C_1 = \{-\frac{\pi}{4}, \frac{\pi}{4}\}^3, \quad CS_1 = \{-3^{-1/2}, 3^{1/2}\}^3; \qquad (5)$$

the eight elements of $CS_1$ are also called the *vertices* (or *corners*) of $CS_N$.

We introduce a subset $E \subset C$ containing the twelve edges of the cube C,

$$E = \{\alpha \in C : |\alpha_1| = |\alpha_2| = \frac{\pi}{4} \text{ or } |\alpha_1| = |\alpha_3| = \frac{\pi}{4} \text{ or } |\alpha_2| = |\alpha_3| = \frac{\pi}{4}\}, \qquad (6)$$

and we define some normal vectors to E as displayed in Figure 2. For any point in the “interior” of an edge, $\alpha \in E \setminus C_1$, there are two faces of C which contain $\alpha$; we introduce one unit normal vector per face, $\nu_1(\alpha)$, $\nu_2(\alpha)$. In the first face, resp. in the second face, the vector $\nu_1(\alpha)$, resp. $\nu_2(\alpha)$, is normal to E at $\alpha$ and points towards the interior of the face; hence, $\nu_i(\alpha)$ is proportional to the vector from the midpoint of the edge to the center of the i-th face, $i = 1, 2$. Also, we fix $w = \nu_1 + \nu_2$, which is proportional with the opposite of the midpoint of the edge.

2.3. **Properties of the correspondence.** The inverse of the mapping $x$ in (3) is given by

$$x^{-1} : x \in \mathbb{S}^2 \mapsto (\arctan \frac{x_1}{|x|_{\infty}}, \arctan \frac{x_2}{|x|_{\infty}}, \arctan \frac{x_3}{|x|_{\infty}}) \in C. \qquad (7)$$

Furthermore, the set

$$\Omega =] -\frac{\pi}{2}, \frac{\pi}{2}[^3 \setminus \{0\} \qquad (8)$$

provides an open neighborhood of C such that the extension

$$x : \alpha \in \Omega \mapsto x(\alpha) = \frac{(\tan \alpha_1, \tan \alpha_2, \tan \alpha_3)}{(\tan^2 \alpha_1 + \tan^2 \alpha_2 + \tan^2 \alpha_3)^{1/2}} \in \mathbb{S}^2 \qquad (9)$$

belongs to $\mathcal{C}^{\infty}(\Omega)$ (but such an extension is not bijective). It is straightforward to check that the jacobian matrix of this extension is given by

$$(\nabla_{\alpha} x(\alpha))^{\top} = \frac{1}{|t|^3} \begin{bmatrix} (1 + t_1^2)(t_2^2 + t_3^2) & -(1 + t_2^2)t_1t_2 & -(1 + t_3^2)t_1t_3 \\ -(1 + t_1^2)t_1t_2 & (1 + t_2^2)(t_1^2 + t_3^2) & -(1 + t_3^2)t_2t_3 \\ -(1 + t_1^2)t_1t_3 & -(1 + t_2^2)t_2t_3 & (1 + t_3^2)(t_1^2 + t_2^2) \end{bmatrix}, \quad \alpha \in \Omega, \qquad (10)$$

where

$$t = (t_1, t_2, t_3) = (\tan \alpha_1, \tan \alpha_2, \tan \alpha_3), \quad |t| = (t_1^2 + t_2^2 + t_3^2)^{1/2}. \qquad (11)$$

3. CORRECTED TRAPEZOIDAL RULE ON THE CUBED SPHERE

In this section, we state the main results of the paper. We present fourth-order and sixth-order quadrature rules on the Cubed Sphere, based on simple corrections of trapezoidal rules.
**3.1. Second-order trapezoidal rule.** Let $f: \mathbb{R}^3 \to \mathbb{R}$ be a sufficiently smooth function in an open neighborhood of $\mathbb{S}^2$. We approximate the spherical integral $\int_{\mathbb{S}^2} f(x)d\sigma$ from the values $f(x), x \in \text{CS}_N$. Taking benefit from the one-to-one correspondence (2), we transform the spherical integral into an integral on the cube $C$, using the change of variable $x = x(\alpha)$:

$$ \int_{\mathbb{S}^2} f(x)d\sigma = \int_C f(x(\alpha))\mu(\alpha)d\sigma, \qquad (12) $$

where the weight $\mu$ is given by

$$ \mu(\alpha) = \frac{(1 + \tan^2 \alpha_1)(1 + \tan^2 \alpha_2)(1 + \tan^2 \alpha_3)}{2(\tan^2 \alpha_1 + \tan^2 \alpha_2 + \tan^2 \alpha_3)^{3/2}}. \qquad (13) $$

Indeed, considering (for instance) the face $\alpha_1 = \frac{\pi}{4}$, it can be checked that

$$ \mu(\alpha) = \|\partial_{\alpha_2}x(\alpha) \wedge \partial_{\alpha_3}x(\alpha)\| = \frac{(1 + \tan^2 \alpha_2)(1 + \tan^2 \alpha_3)}{(1 + \tan^2 \alpha_2 + \tan^2 \alpha_3)^{3/2}}, \quad \alpha \in C \text{ with } \alpha_1 = \frac{\pi}{4}. $$

On the right hand side of (12), the integral can be decomposed as the sum of the six integrals on the six faces of the cube $C$. For each of those integrals, we apply a bivariate trapezoidal rule, for the function $F = (f \circ x) \cdot \mu$ restricted to the considered face, in the square $[a, b]^2 = [-\frac{\pi}{4}, \frac{\pi}{4}]^2$, and the step $h$. Summing over the six faces yields

$$ \int_C f(x(\alpha))\mu(\alpha)d\sigma = h^2 \sum_{\alpha \in C_N} f(x(\alpha))\mu(\alpha)(1 - \frac{1}{4}\mathbb{1}_{\alpha \in C_1}) + R_h(f), \qquad (14) $$

where $R_h(f)$ denotes the total remainder. Equivalently to (14), we obtain the following spherical quadrature rule,

$$ \int_{\mathbb{S}^2} f(x)d\sigma = h^2 \sum_{x \in \text{CS}_N} f(x)\omega(x)(1 - \frac{1}{4}\mathbb{1}_{x \in \text{CS}_1}) + R_h(f), \qquad (\text{T-CS}) $$

where $\omega(x) = \mu(x^{-1}(x))$. This rule is called a *trapezoidal rule on the Cubed Sphere* due to its derivation. Note that the factor $(1 - \frac{1}{4}\mathbb{1}_{x \in \text{CS}_1})$ takes the value $\frac{3}{4}$ on the eight corners $x \in \text{CS}_1$, and the value 1 otherwise. Also, the combination of (7) and (13) reveals that

$$ \omega(x) = \frac{1}{2}|x|_{\infty}^3 \left(1 + \frac{x_1^2}{|x|_{\infty}^2}\right) \left(1 + \frac{x_2^2}{|x|_{\infty}^2}\right) \left(1 + \frac{x_3^2}{|x|_{\infty}^2}\right), \quad x = (x_1, x_2, x_3) \in \text{CS}_N. \qquad (15) $$

A particular question of interest concerns the analysis of the remainder $R_h(f)$. As a consequence of the Euler-Maclaurin formula in Proposition 8, it is straightforward to prove that if $f \in \mathcal{C}^2(O)$ in an open neighborhood $O$ of $\mathbb{S}^2$, then $R_h(f)$ satisfies

$$ R_h(f) = \mathcal{O}(h^2). $$

**3.2. Fourth-order rule.** Numerical results in previous works suggest that very simple corrections of the rule in (T-CS) increase the order of accuracy. Indeed, in [14, 19, 20, 32], the special factor (3/4) occurring at the corners ($x \in \text{CS}_1$) has been ingeniously removed. This trick permits to integrate the constant function 1 with a fourth-order accuracy (proof in [14]), while a surprising fourth-order accuracy has been observed in the numerical experiments, for any smooth function $f$. This suggests that the following relation may be true,

$$ -\frac{h^2}{4} \sum_{x \in \text{CS}_1} f(x)\omega(x) + R_h(f) = \mathcal{O}(h^4). $$

This claim, due to the method introduced in [19, 20], has never been proved so far. The first main result of this paper is a proof of this claim. We prove indeed the following theorem.

**Theorem 1 (Fourth-order rule).** Let $\text{CS}_N$ be the Cubed Sphere with step $h$, defined in (1), and consider the weight $\omega$ in (15). For every $f \in \mathcal{C}^4(O)$ with $O \subset \mathbb{R}^3$ an open neighborhood of the sphere $\mathbb{S}^2$,

$$ \int_{\mathbb{S}^2} f(x)d\sigma = h^2 \sum_{x \in \text{CS}_N} \omega(x)f(x) + \mathcal{O}(h^4). \qquad (\text{T-CS4}) $$
Remark 2. On the cube C, we get a uniform quadrature rule on uniform cartesian grids,

$$ \int_C F(\alpha) d\sigma = h^2 \sum_{\alpha \in C_N} F(\alpha) + \mathcal{O}(h^4), \quad \text{with} \quad F(\alpha) = f(x(\alpha))\mu(\alpha); $$

the remarkable fourth-order of accuracy is due to the form of F.

**3.3. Sixth-order rule.** A further correction has been mentioned in [8], as follows. Consider the remainder $R_h(1)$ in (T-CS) with $f = 1$. Add $\frac{1}{8}R_h(1) \sum_{x \in CS_1} f(x)$ to the trapezoidal rule in (T-CS), in order to exactly integrate the constant function 1. In practice, the value $R_h(1)$ in (T-CS) is computed once for all in $\mathcal{O}(N^2)$ operations; then, the proposed correction is observed to achieve a sixth-order of accuracy for any smooth function $f$. This suggests that the following claim may be true,

$$ R_h(f) = \frac{1}{8} R_h(1) \sum_{x \in CS_1} f(x) + \mathcal{O}(h^6). \quad (16) $$

In particular, it should be possible to correct the fourth-order term of Theorem 1, using only the values on the corners, $f(x)$, $x \in CS_1$. The second main result of this paper tackles this problem. It provides a framework in which

$$ -\frac{h^2}{4} \sum_{x \in CS_1} f(x)\omega(x) + R_h(f) = -\frac{h^4}{15} \sum_{x \in CS_1} f(x)\omega(x) + \mathcal{O}(h^6). $$

Hence, the sixth-order accuracy is achieved using only a correction on the corners, as follows.

**Theorem 3 (Sixth-order rule).** Assume that $f: \mathbb{R}^3 \to \mathbb{R}$ can be expanded as a multivariate series converging absolutely in $[-\rho, \rho]^3$, for some $\rho > 1$. Then,

$$ \int_{S^2} f(x)d\sigma = Q_h(f) + \mathcal{O}(h^6), \quad \text{with, } Q_h(f) = h^2 \sum_{x \in CS_N} \omega(x)\left(1 - \frac{h^2}{15}\mathbb{I}_{x \in CS_1}\right)f(x), \quad (\text{T-CS6}) $$

$CS_N$ and $h$ in (1), $\omega$ in (15), and $CS_1$ in (5).

Remark 4. As a corollary, the claim (16) was true.

## 4. PROOF OF THE THEOREMS

In this section, we prove Theorems 1-3. In brief, we derive an Euler-Maclaurin formula on the Cubed Sphere; it quantifies how the remainder $R_h(f)$ is due to edges and vertices. In a second step, the obtained formula is expressed as a combination of the values on the corners $f(x)$, $x \in CS_1$, using algebraic manipulation, or symbolic computation.

**Lemma 5 (Euler-Maclaurin formula on the Cubed Sphere).** For all $k \ge 0$ and $f \in \mathcal{C}^{2k+2}(S^2)$, the remainder $R_h(f)$ in (T-CS) expands as follows,

$$ R_h(f) = \sum_{p=1}^{k} r_p(f)h^{2p} + \mathcal{O}(h^{2k+2}), \quad \text{with,} \quad (17) $$

$$ \begin{aligned} r_p(f) = & \frac{b_{2p}}{(2p)!} \int_E (\partial_{\nu_1}^{2p-1} + \partial_{\nu_2}^{2p-1}) [f(x(\alpha))\mu(\alpha)] d\ell \\ & - \sum_{\substack{i+j=p \\ 1 \le i,j \le p-1}} \frac{b_{2i} b_{2j}}{(2i)! (2j)!} \sum_{\alpha \in CS_1} \left[ \begin{aligned} & \operatorname{sg}(\alpha_1 \alpha_2) \partial_{\alpha_1}^{2i-1} \partial_{\alpha_2}^{2j-1} \\ & + \operatorname{sg}(\alpha_2 \alpha_3) \partial_{\alpha_2}^{2i-1} \partial_{\alpha_3}^{2j-1} \\ & + \operatorname{sg}(\alpha_1 \alpha_3) \partial_{\alpha_1}^{2i-1} \partial_{\alpha_3}^{2j-1} \end{aligned} \right] [f(x(\alpha))\mu(\alpha)]. \end{aligned} \quad (18) $$

Here, $b_{2i}$, $1 \le i \le p$, are Bernoulli numbers, $x \in \mathcal{C}^\infty(\Omega)$ is given in (9), $\mu \in \mathcal{C}^\infty(\Omega)$ is given in (13), the set of edges $E$, the set of vertices $C_1$, and the normal vectors $\nu_1, \nu_2$ are displayed in Figure 2.
*Proof.* We introduce $F : \alpha \in \Omega \mapsto f(x(\alpha))\mu(\alpha) \in \mathbb{R}$, with $\Omega$ in (8). This function belongs to $\mathcal{C}^{2k+2}(\Omega)$. On each face of $C$, we apply the bivariate trapezoidal rule with step $h$, on the square $[a, b]^2 = [-\frac{\pi}{4}, \frac{\pi}{4}]^2$, for the function $F$ restricted to the face. We estimate the remainder using the bi-dimensional Euler-Maclaurin formula in (37). Then, we sum over the six faces to get the total remainder. The contribution of the normal derivatives along the boundary of the faces forms an integral on $E$; along each edge, the contribution of two adjacent faces are added. On each corner, the contribution of the cross-derivatives from three faces are added. $\square$

*Proof of Theorem 1.* Due to Lemma 5 with $k = 2$,

$$R_h(f) = r_2(f)h^2 + \mathcal{O}(h^4), \quad \text{with} \quad r_2(f) = \frac{1}{2^6} \int_E (\partial_{\nu_1} + \partial_{\nu_2})[f(x(\alpha))\mu(\alpha)]d\ell. \quad (19)$$

We introduce the vector $w = \nu_1 + \nu_2$ (displayed in Figure 2), so that

$$\partial_{\nu_1} + \partial_{\nu_2} = \partial_w = w^\top \nabla$$

represents a derivative in the direction $w$ (which points towards the interior of $C$).

The integral in (19) can be decomposed as the sum over the twelve edges in $E$. Consider one of the edges, denoted by $E_0$, and compute

$$\int_{E_0} \partial_w [f(x(\alpha))\mu(\alpha)] d\ell. \quad (20)$$

Up to a permutation, along the edge $E_0$,

$$|\alpha_1| = |\alpha_2| = \frac{\pi}{4}, \quad \alpha_3 \in [-\frac{\pi}{4}, \frac{\pi}{4}], \quad w = [-\operatorname{sg}(\alpha_1), -\operatorname{sg}(\alpha_2), 0]^\top, \quad \alpha \in E_0; \quad (21)$$

(on the endpoints of $E_0$, we also define $w$ by this relation). Due to Lemma 6 stated hereafter the proof, the integrated function $\partial_w[f(x(\alpha))\mu(\alpha)]$ satisfies

$$\partial_w[f(x(\alpha))\mu(\alpha)] = \partial_{\alpha_3}[f(x(\alpha))2x_3(\alpha)(1 - x_3(\alpha)^2)], \quad \alpha \in E_0. \quad (\star)$$

The equality $(\star)$ is the key of the proof. The integrated function along $E_0$ in (20) is a derivative in the direction $w$, which enters inside the cube $C$. The relation $(\star)$ expresses this function as a derivative in the direction of the edge $E_0$. It enables the computation as follows,

$$\begin{align*}
\int_{E_0} \partial_w [f(x(\alpha))\mu(\alpha)] d\ell &= \int_{-\frac{\pi}{4}}^{\frac{\pi}{4}} \partial_{\alpha_3} [f(x(\alpha))2x_3(\alpha)(1-x_3(\alpha)^2)] d\alpha_3 \\
&= [f(x(\alpha))2x_3(\alpha)(1-x_3(\alpha)^2)]^{\alpha_3=\frac{\pi}{4}}_{\alpha_3=-\frac{\pi}{4}} \\
&= \frac{4}{3\sqrt{3}} [f(\operatorname{sg} \alpha_1 \frac{1}{\sqrt{3}}, \operatorname{sg} \alpha_2 \frac{1}{\sqrt{3}}, \frac{1}{\sqrt{3}}) + f(\operatorname{sg} \alpha_1 \frac{1}{\sqrt{3}}, \operatorname{sg} \alpha_2 \frac{1}{\sqrt{3}}, -\frac{1}{\sqrt{3}})],
\end{align*}$$

where the two nodes correspond to the endpoints of the edge $E_0$. Similar results apply for the other edges. Hence, summing over the twelve edges yields

$$r_2(f) = \frac{1}{3\sqrt{3}} \sum_{\alpha \in C_1} f(x(\alpha)) = \frac{1}{4} \sum_{x \in C_1} f(x)\omega(x). \quad (22)$$

The combination of this result and (19) achieves the proof. $\square$

**Lemma 6.** Let $f \in \mathcal{C}^1(O)$ for some open neighborhood $O$ of $\mathbb{S}^2$. Let $E_0$ be an edge of $C$ with direction $(0, 0, 1)$; along $E_0$, fix $w = \nu_1 + \nu_2$ as in Figure 2. Then,

$$\partial_w[f(x(\alpha))\mu(\alpha)] = \partial_{\alpha_3}[f(x(\alpha))\partial_w x_3(\alpha)], \quad \alpha \in E_0, \quad (23)$$

$$\partial_w x_3(\alpha) = 2x_3(\alpha)(1 - x_3(\alpha)^2), \quad \alpha \in E_0. \quad (24)$$

*Proof.* Along the edge $E_0$, $\alpha$ and $w$ satisfy (21). By derivation of products, we obtain that

$$\begin{align*}
\partial_w[f(x(\alpha))\mu(\alpha)] &= \mu(\alpha)\nabla_x f(x(\alpha))^{\top}\partial_w x(\alpha) + f(x(\alpha))\partial_w\mu(\alpha), \quad \alpha \in \Omega, \\
\partial_{\alpha_3}[f(x(\alpha))\partial_w x_3(\alpha)] &= \nabla_x f(x(\alpha))^{\top}(\partial_{\alpha_3}x(\alpha))\partial_w x_3(\alpha) + f(x(\alpha))\partial_{\alpha_3}\partial_w x_3(\alpha), \quad \alpha \in \Omega.
\end{align*}$$
We check that these two derivatives are equal if $\alpha \in E_0$. Due to the two relations above, it is in fact sufficient to check that

$$ \partial_w \mu = \partial_{\alpha_3} \partial_w x_3, \quad \alpha \in E_0, \qquad (25) $$

$$ \mu \partial_w x = (\partial_{\alpha_3} x)(\partial_w x_3), \quad \alpha \in E_0. \qquad (26) $$

For any $\alpha \in E_0$, the jacobian matrix $(\nabla_{\alpha} x(\alpha))^{\top}$ is given by (10) with $t = (\operatorname{sg} \alpha_1, \operatorname{sg} \alpha_2, \operatorname{tan} \alpha_3)$. Considering the last element of $(\nabla_{\alpha} x(\alpha))^{\top}$, we observe that

$$ \mu(\alpha) = \partial_{\alpha_3} x_3(\alpha), \quad \alpha \in E_0. \qquad (27) $$

Considering the first two columns of $(\nabla_{\alpha} x(\alpha))^{\top}$, we obtain

$$ \begin{aligned} \partial_w x(\alpha) = \nabla x(\alpha)^{\top} w &= -\operatorname{sg}(\alpha_1)\partial_{\alpha_1} x(\alpha) - \operatorname{sg}(\alpha_2)\partial_{\alpha_2} x(\alpha), \quad \alpha \in E_0, \\ &= \frac{2 \operatorname{tan} \alpha_3}{(2 + \operatorname{tan}^2 \alpha_3)^{3/2}} \begin{bmatrix} -\operatorname{sg} \alpha_1 \operatorname{tan} \alpha_3 \\ -\operatorname{sg} \alpha_2 \operatorname{tan} \alpha_3 \\ 2 \end{bmatrix} = \begin{bmatrix} -2x_1(\alpha)x_3(\alpha)^2 \\ -2x_2(\alpha)x_3(\alpha)^2 \\ 2x_3(\alpha)(1 - x_3(\alpha)^2) \end{bmatrix}, \end{aligned} $$

where the third component gives the relation (24).

Next, we differentiate the relation (24) with respect to $\alpha_3$ using the chain rule, and we obtain with (27) that

$$ \partial_{\alpha_3} \partial_w x_3(\alpha) = 2(1 - 3x_3(\alpha)^2)\partial_{\alpha_3} x_3(\alpha) = 2(1 - 3x_3(\alpha)^2)\mu(\alpha), \quad \alpha \in E_0. $$

Furthermore, the gradient of $\mu$ is given by

$$ \nabla_{\alpha} \mu(\alpha) = \frac{\mu(\alpha)}{|t|^2} \begin{bmatrix} t_1(2t_2^2 + 2t_3^2 - 3 - t_1^2) \\ t_2(2t_1^2 + 2t_3^2 - 3 - t_2^2) \\ t_3(2t_1^2 + 2t_2^2 - 3 - t_3^2) \end{bmatrix}, \quad \text{with } t \text{ in } (11), \quad \alpha \in \Omega. $$

For $\alpha \in E_0$, $t = (\operatorname{sg} \alpha_1, \operatorname{sg} \alpha_2, \operatorname{tan} \alpha_3)$, and we deduce

$$ \begin{aligned} \partial_w \mu(\alpha) = \nabla \mu(\alpha)^{\top} w &= -\operatorname{sg}(\alpha_1)\partial_{\alpha_1} \mu(\alpha) - \operatorname{sg}(\alpha_2)\partial_{\alpha_2} \mu(\alpha) \\ &= 2\mu(\alpha)(1 - 3x_3(\alpha)^2), \quad \alpha \in E_0. \end{aligned} $$

We recognize the expression above, which proves (25).

Lastly, with the third column of the jacobian matrix $(\nabla_{\alpha} x(\alpha))^{\top}$, we observe that

$$ \partial_{\alpha_3} x(\alpha) = \frac{1 + \tan^2 \alpha_3}{(2 + \tan^2 \alpha_3)^{3/2}} \begin{bmatrix} -\operatorname{sg} \alpha_1 \operatorname{tan} \alpha_3 \\ -\operatorname{sg} \alpha_2 \operatorname{tan} \alpha_3 \\ 2 \end{bmatrix} = \frac{1}{2} \mu(\alpha) \begin{bmatrix} -\operatorname{sg} \alpha_1 \operatorname{tan} \alpha_3 \\ -\operatorname{sg} \alpha_2 \operatorname{tan} \alpha_3 \\ 2 \end{bmatrix}, \quad \alpha \in E_0. $$

We recognize the vector $\partial_w x(\alpha)$ computed above, up to a proportionality coefficient. Considering the third line, the proportionality coefficient is observed to be $\frac{\partial_w x_3}{\mu}$, and we finally get (26). $\square$

*Proof of Theorem 3.* First, assume (only) that $f \in \mathscr{C}^6(\mathcal{O})$ for some open neighborhood $\mathcal{O}$ of $\mathbb{S}^2$. Due to Lemma 5 with $k = 4$, the remainder $R_h(f)$ in (T-CS) satisfies

$$ R_h(f) = r_2(f)h^2 + r_4(f)h^4 + \mathcal{O}(h^6), \qquad (28) $$

with $r_2(f)$ in (19) and

$$ r_4(f) = -\frac{1}{30 \cdot 4!} \int_E (\partial_{\nu_1}^3 + \partial_{\nu_2}^3) [f(x(\alpha))\mu(\alpha)] d\ell \\ - \frac{1}{(2 \cdot 6)^2} \sum_{\alpha \in V} [\operatorname{sg}(\alpha_1 \alpha_2)\partial_{\alpha_1} \partial_{\alpha_2} + \operatorname{sg}(\alpha_2 \alpha_3)\partial_{\alpha_2} \partial_{\alpha_3} + \operatorname{sg}(\alpha_1 \alpha_3)\partial_{\alpha_1} \partial_{\alpha_3}] [f(x(\alpha))\mu(\alpha)]. \quad (29) $$

For the second order term, Theorem 1 shows (22). The rest of the proof aims at expressing the fourth-order coefficient $r_4(f)$ in term of the values $f(x)$, $x \in \mathcal{CS}_1$.
For that purpose, we start with $f = 1$. We consider the edge $E_0 = \{(\frac{\pi}{4}, \frac{\pi}{4}, \alpha_3), \alpha_3 \in [-\frac{\pi}{4}, \frac{\pi}{4}]\}$, $\partial_{\nu_2} = -\partial_{\alpha_1}$ (normal derivative along $E_0$, for the face $\alpha_2 = \frac{\pi}{4}$), and the vertex $v_0 = (\frac{\pi}{4}, \frac{\pi}{4}, \frac{\pi}{4})$. A symbolic computation (with Maple) has shown that

$$ \int_{E_0} \partial_{\nu_2}^3 \mu(\alpha) d\ell = \frac{32\sqrt{3}}{27}, \quad \partial_{\alpha_1} \partial_{\alpha_3} \mu(v_0) = \frac{32\sqrt{3}}{27}. $$

Taking benefit from symmetry, it is expected that

$$ r_4(1) = -\frac{1}{30 \cdot 4!} 12 \times \left(2 \times \frac{32\sqrt{3}}{27}\right) - \frac{1}{(2 \cdot 6)^2} 8 \times \left(3 \times \frac{32\sqrt{3}}{27}\right) = -8 \frac{4\sqrt{3}}{135}; $$

this quantity can be distributed over the eight vertices as follows,

$$ r_4(1) = -\frac{4\sqrt{3}}{135} \sum_{x \in \text{CS}_1} 1 = -\frac{1}{15} \sum_{x \in \text{CS}_1} \omega(x). $$

Therefore, if there is a constant $C$ such that $r_4(f) = C \sum_{x \in \text{CS}_1} \omega(x)f(x)$, it is expected that $C = -\frac{1}{15}$.

This approach has enabled to guess that more generally, the term $r_4(f)$ should satisfy

$$ r_4(f) = -\frac{1}{15} \sum_{x \in \text{CS}_1} \omega(x) f(x), \qquad (30) $$

which is in fact the result to be proved to achieve the proof. Unfortunately, we do not have at disposal a direct proof of this relation for a general function $\mathcal{C}^6(O)$. That is the reason why we add the assumption that $f$ can be expanded as a series converging absolutely in $[-\rho, \rho]^3$ with $\rho > 1$. In this case, $f \in \mathcal{C}^\infty([-\rho, \rho]^3)$, so (28) is valid, and the relation (30) has been checked in Lemma 7, stated hereafter. $\square$

**Lemma 7.** Assume that $f : \mathbb{R}^3 \to \mathbb{R}$ can be expanded as a multivariate series converging absolutely in $[-\rho, \rho]^3$, for some $\rho > 1$. Then, the quantity $r_4(f)$ defined in (29) satisfies (30).

*Proof.* In a first step, assume that $f$ is the multivariate polynomial $f(x) = x^d = x_1^{d_1} x_2^{d_2} x_3^{d_3}$, with $d = (d_1, d_2, d_3) \in \mathbb{N}^3$; the total degree of such a polynom is given by $|d| = d_1 + d_2 + d_3$. In this case, the relation (30) can be checked using symbolic computation (Maple code available in Appendix B). This intermediate step implies by linearity that (30) is true for every polynomial.

Now, consider a multivariate series $f(x) = \sum_{d \in \mathbb{N}^3} \hat{f}_d x^d$, with $\hat{f}_d \in \mathbb{R}$, $d \in \mathbb{N}^3$. Assume that the series converges absolutely in $[-\rho, \rho]^3$, i.e. $\sum_{d \in \mathbb{N}^3} |\hat{f}_d| |\rho|^d| < \infty$, for some $\rho > 1$. Fix $r$ such that $1 < r < \rho$. Then $f \in \mathcal{C}^\infty([-\rho, \rho]^3)$. The partial series $f_D(x) = \sum_{|d| \le D} \hat{f}_d x^d$, $D \in \mathbb{N}$ converge uniformly to $f$ on $[-r, r]^3$ (when $D \to \infty$); and for any order of derivation $\beta = (\beta_1, \beta_2, \beta_3) \in \mathbb{N}^3$, $\partial^\beta f_D$ converges uniformly to $\partial^\beta f$ on $[-r, r]^3$. We apply the relation (30) for the polynomial $f_D$, and we obtain (30) when $D \to \infty$, due to the uniform convergence of $f_D$ and its derivatives. $\square$

## 5. NUMERICAL RESULTS

5.1. **Observation of the sixth-order.** We introduce a set of test functions in Table 1 with various smoothness properties. The function $f_{-1}$ is piecewise constant; it is discontinuous across the hyperplane $u^\top x = c$, with

$$ c = -0.324, \quad u = \frac{v}{\|v\|}, \quad v = (1, 0.2, -0.376). \qquad (31) $$

For every $0 \le k \le 7$, $f_k \in \mathcal{C}^k(\mathbb{R}^3)$; the $k$-th directional derivative $\partial_u^k f_k$ is proportional to $f_0$ which is not differentiable across $u^\top x = c$. Lastly, $f_\infty \in \mathcal{C}^\infty(\mathbb{R}^3)$, as a multivariate series converging absolutely in $\mathbb{R}^3$.

For each test function $f_k$, the exact value of $\int_{\mathbb{S}^2} f_k d\sigma$ is reported in Table 1. We apply the quadrature rule studied in Theorem 3,

$$ Q_h(f_k) = h^2 \sum_{x \in \text{CS}_N} \omega(x) \left(1 - \frac{h^2}{15} \mathbb{1}_{x \in \text{CS}_1}\right) f_k(x). $$
<table><caption>TABLE 1. Test functions with various smoothness properties (with c, u defined in (31)).</caption><thead><tr><th>k</th><th>f<sub>k</sub>(x)</th><th>∫<sub>S<sup>2</sup></sub> f<sub>k</sub> dσ</th><th>smoothness</th></tr></thead><tbody><tr><td>k = -1</td><td>1<sup>l</sup><sub>u<sup>l</sup>x≥c</sub></td><td>2π(1 - c)</td><td>piecewise constant</td></tr><tr><td>0 ≤ k ≤ 7</td><td>(u<sup>l</sup>x - c)<sup>k+1</sup>1<sup>l</sup><sub>u<sup>l</sup>x≥c</sub></td><td>2π(1 - c)<sup>k+2</sup></td><td>C<sup>k</sup>(R<sup>3</sup>)</td></tr><tr><td>k = ∞</td><td>exp(u<sup>l</sup>x)</td><td>2π(exp(1) - exp(-1))</td><td>series with infinite radius</td></tr></tbody></table>

<table><caption>TABLE 2. Quadrature errors of the test functions $f_k$ in Table 1. The relative quadrature errors $\epsilon_k(h)$ in (32) and the corresponding “numerical orders” of convergence $p_k(h)$ in (33) are reported. Convergence is always observed. A sixth order of convergence is observed if $f_k$ is at least $C^4$ ($4 \le k \le 7$ of $k = \infty$).</caption><thead><tr><th>N</th><th>1</th><th>2</th><th>4</th><th>8</th><th>16</th><th>32</th><th>64</th><th>128</th><th>256</th></tr></thead><tbody><tr><td>ε<sub>-1</sub>(h)</td><td>4.6e-02</td><td>6.3e-02</td><td>2.1e-02</td><td>1.4e-03</td><td>1.0e-03</td><td>1.2e-04</td><td>3.0e-04</td><td>3.6e-05</td><td>1.9e-05</td></tr><tr><td>p<sub>-1</sub>(h)</td><td>×</td><td>-0.4</td><td>1.6</td><td>3.9</td><td>0.4</td><td>3.2</td><td>-1.4</td><td>3.1</td><td>0.9</td></tr><tr><td>ε<sub>0</sub>(h)</td><td>1.4e-02</td><td>6.6e-03</td><td>7.8e-04</td><td>1.6e-04</td><td>1.9e-05</td><td>1.8e-06</td><td>1.8e-07</td><td>8.6e-08</td><td>2.7e-08</td></tr><tr><td>p<sub>0</sub>(h)</td><td>×</td><td>1.1</td><td>3.1</td><td>2.2</td><td>3.1</td><td>3.3</td><td>3.4</td><td>1.0</td><td>1.6</td></tr><tr><td>ε<sub>1</sub>(h)</td><td>2.5e-02</td><td>2.4e-03</td><td>1.7e-04</td><td>2.9e-06</td><td>2.3e-07</td><td>5.0e-08</td><td>5.9e-09</td><td>1.4e-10</td><td>4.5e-11</td></tr><tr><td>p<sub>1</sub>(h)</td><td>×</td><td>3.4</td><td>3.8</td><td>5.9</td><td>3.6</td><td>2.2</td><td>3.1</td><td>5.4</td><td>1.6</td></tr><tr><td>ε<sub>2</sub>(h)</td><td>4.5e-03</td><td>1.6e-03</td><td>3.6e-05</td><td>3.2e-07</td><td>2.1e-08</td><td>3.1e-10</td><td>1.4e-11</td><td>1.3e-12</td><td>1.2e-13</td></tr><tr><td>p<sub>2</sub>(h)</td><td>×</td><td>1.5</td><td>5.5</td><td>6.8</td><td>3.9</td><td>6.1</td><td>4.5</td><td>3.4</td><td>3.4</td></tr><tr><td>ε<sub>3</sub>(h)</td><td>5.2e-02</td><td>2.7e-04</td><td>4.6e-07</td><td>8.7e-08</td><td>2.1e-09</td><td>5.3e-11</td><td>1.7e-13</td><td>9.2e-15</td><td>5.2e-16</td></tr><tr><td>p<sub>3</sub>(h)</td><td>×</td><td>7.6</td><td>9.2</td><td>2.4</td><td>5.4</td><td>5.3</td><td>8.3</td><td>4.2</td><td>4.1</td></tr><tr><td>ε<sub>4</sub>(h)</td><td>1.0e-01</td><td>2.9e-04</td><td>8.1e-06</td><td>9.2e-08</td><td>1.3e-09</td><td>2.2e-11</td><td>3.4e-13</td><td>5.4e-15</td><td>0.0e+00</td></tr><tr><td>p<sub>4</sub>(h)</td><td>×</td><td>8.5</td><td>5.1</td><td>6.5</td><td>6.1</td><td>6.0</td><td>6.0</td><td>6.0</td><td>×</td></tr><tr><td>ε<sub>5</sub>(h)</td><td>1.6e-01</td><td>5.6e-04</td><td>2.0e-05</td><td>3.1e-07</td><td>4.6e-09</td><td>7.1e-11</td><td>1.1e-12</td><td>1.7e-14</td><td>1.4e-16</td></tr><tr><td>p<sub>5</sub>(h)</td><td>×</td><td>8.1</td><td>4.8</td><td>6.0</td><td>6.1</td><td>6.0</td><td>6.0</td><td>6.0</td><td>7.0</td></tr><tr><td>ε<sub>6</sub>(h)</td><td>2.1e-01</td><td>4.5e-04</td><td>3.3e-05</td><td>5.1e-07</td><td>7.8e-09</td><td>1.2e-10</td><td>1.9e-12</td><td>3.0e-14</td><td>3.6e-16</td></tr><tr><td>p<sub>6</sub>(h)</td><td>×</td><td>8.9</td><td>3.8</td><td>6.0</td><td>6.0</td><td>6.0</td><td>6.0</td><td>6.0</td><td>6.4</td></tr><tr><td>ε<sub>7</sub>(h)</td><td>2.6e-01</td><td>3.8e-03</td><td>4.6e-05</td><td>7.3e-07</td><td>1.1e-08</td><td>1.7e-10</td><td>2.7e-12</td><td>4.2e-14</td><td>8.1e-16</td></tr><tr><td>p<sub>7</sub>(h)</td><td>×</td><td>6.1</td><td>6.4</td><td>6.0</td><td>6.0</td><td>6.0</td><td>6.0</td><td>6.0</td><td>5.7</td></tr><tr><td>ε<sub>∞</sub>(h)</td><td>9.2e-03</td><td>9.1e-04</td><td>2.0e-05</td><td>3.4e-07</td><td>5.5e-09</td><td>8.6e-11</td><td>1.3e-12</td><td>2.1e-14</td><td>8.4e-16</td></tr><tr><td>p<sub>∞</sub>(h)</td><td>×</td><td>3.3</td><td>5.5</td><td>5.9</td><td>6.0</td><td>6.0</td><td>6.0</td><td>6.0</td><td>4.6</td></tr></tbody></table>

The associated relative errors

$$ \epsilon_k(h) = \frac{|\int_{S^2} f_k d\sigma - Q_h(f_k)|}{\int_{S^2} f_k d\sigma}, \quad \text{with} \quad h = \frac{\pi}{2N}, \quad N \in \{1, 2, 4, 8, 16, 32, 64, 128, 256\}, \quad (32) $$

have been computed in double precision accuracy in Matlab, and have been reported in Table 2. To study the behavior of the error, we have also reported the exponent $p_k(h)$ such that $\epsilon_k(h/2) = 2^{-p_k(h)}\epsilon_k(h)$, i.e.

$$ p_k(h) = \log_2 \epsilon_k(h) - \log_2 \epsilon_k(h/2). \quad (33) $$

This exponent represents the numerical order of convergence of $\epsilon_k(h)$ (if it converges).

For the function $f_\infty$, which satisfies the assumption of Theorem 3, the numerical order of convergence is six, as expected theoretically ($\epsilon_\infty(h) = O(h^6)$). A sixth-order of convergence is still observed for the functions $f_k \in C^k$, $4 \le k \le 7$, despite they do not satisfy the assumption of Theorem 3. Therefore, for these test functions, the $C^4$-smoothness is sufficient to reach the sixth-order accuracy. For the functions $f_k$, $-1 \le k \le 3$, ranging from discontinuous to $C^3$, some convergence is still observed; the convergence rate increases with the smoothness, but it is difficult to state a clear rule.
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Display of the function $|g| + 10^{-14}$ on $\mathbb{S}^2$ (log$_{10}$-scale), with $g \in \mathcal{C}^{\infty}$ in (34). The function $g$ and its derivatives are not trivial around the set $x(E)$ of the edges (displayed in bold gray), but any derivative of $g$ is null on the set $\text{CS}_1$ of the vertices.</figcaption></figure>

<table><caption>TABLE 3. Quadrature errors of the function $g$ in (34), whose derivatives on the set $\text{CS}_1$ of the vertices are null. The estimated relative errors $\eta(h)$ and the corresponding “numerical orders” of convergence $p(h)$ in (35) are reported. An ultra-fast convergence appears, despite $g$ is not trivial around the edges.</caption><thead><tr><th>N</th><th>1</th><th>2</th><th>4</th><th>8</th><th>16</th><th>32</th><th>64</th></tr></thead><tbody><tr><td>η(h)</td><td>1.0</td><td>6.2</td><td>7.3e-01</td><td>1.1e-03</td><td>2.6e-08</td><td>1.8e-14</td><td>1.4e-16</td></tr><tr><td>p(h)</td><td>×</td><td>-2.6</td><td>3.1</td><td>9.3</td><td>15.4</td><td>20.5</td><td>7.0</td></tr></tbody></table>

5.2. **Invisible contribution of the edges.** Due to Lemma 5, it is obvious that if $f \in \mathcal{C}^{\infty}$ “does not see” the edges, with $f(x(\alpha))$ supported in the interior of the faces of $C$, then the remainder $R_h(f)$ converges to any order. In this section, we study the case of a function which “does not see” the vertices, but “sees” the edges (and the faces).

Consider a function $g \in \mathcal{C}^{\infty}(\mathbb{R}^3)$ with null derivatives on the set $\text{CS}_1$ of the vertices, defined by

$$g(x) = \cos(9.2 u^\top x) \prod_{z \in \text{CS}_1} \mathbb{1}_{z^\top x < 1} \exp \left( \frac{1}{2} + \frac{1}{z^\top x - 1} \right), \quad u = \frac{v}{\|v\|}, \quad v = (1, 0.2, -0.376). \quad (34)$$

As displayed in Figure 3, the function $g$ is flattened around the vertices, but is not trivial around the set $x(E)$ of the edges.

In Table 3, we have reported estimated errors, and associated “numerical orders” of convergence,

$$\eta(h) = \frac{|Q_{\bar{h}}(g) - Q_h(g)|}{|Q_{\bar{h}}(g)|}, \quad p(h) = \log_2 \eta(h) - \log_2 \eta(h/2), \quad (35)$$

where $Q_{\bar{h}}(g)$, with $\bar{h} = \frac{\pi}{2 \cdot 128}$ (corresponding to $N = 128$), is considered as a reference value of $\int_{\mathbb{S}^2} g d\sigma$. The computation has been performed in double precision accuracy in Matlab.

An ultra-fast convergence is observed. This suggests that in the expression (17) of the remainder $R_h(g)$, many terms $r_p(g)$, with $r_p$ defined in (18), may be null. In the expression of $r_p(g)$, the cross-derivatives on $C_1$ are null; the fast convergence suggests that many integrals on the set $E$ of the edges may be null too. In brief, as soon as a smooth function “does not see” the vertices, an ultra fast-convergence may be expected, even if the function “sees” the edges.

## 6. CONCLUSION AND DISCUSSION

This paper has dealt with quadrature rules on the Cubed Sphere, based on the trapezoidal rule with suitable correction on the eight vertices. We have proved that the rule introduced in [19,20] has the fourth-order accuracy, and we have proposed a further correction which achieves the sixth-order accuracy. Numerical experiments confirm this order of accuracy for smooth-enough functions (at least $\mathcal{C}^4$). These two rules obviously satisfy the following further properties. They
contain positive weights. And they are invariant under the octahedral group, so they exactly integrate a large class of functions, including 15/16 of the real spherical harmonics [9].

The fourth-order of accuracy of (T-CS4) is due to the relation (23), and has been proved in $\mathcal{C}^4$. The sixth-order of accuracy of (T-CS6) has been proved for polynomials and then extended to multivariate series; in particular, it is achieved for any spherical harmonics. Giving a more direct proof for $\mathcal{C}^6$-functions, with relations analogous to (23), is still open.

Also, the numerical experiments suggest that the remainder term may be explained by the behavior around the vertices only. An ultra-fast convergence has been observed for a smooth function which has been flattened around the vertices. This observation is surprising due our expression of the remainder term in (17). So justifying this fast convergence is still open.

At the end, new perspectives are emerging from this work. It may be possible to increase the order of accuracy of the rule (T-CS6), using additional correction based on the values of $f$ in a neighborhood of the vertices.

## APPENDIX A. TWO-DIMENSIONAL EULER-MACLaurin FORMULA

In this section, we derive a two-dimensional Euler-Maclaurin formula in order to estimate the remainder of a bivariate trapezoidal rule. The classical one-dimensional Euler-Maclaurin formula estimates the remainder of the usual trapezoidal rule as follows. For all $a < b \in \mathbb{R}$, $f \in \mathcal{C}^{2k+2}([a, b])$ and $N \ge 1$, with $h = \frac{b-a}{N}$, there exists $\xi \in (a, b)$ such that

$$ \int_a^b f(x)dx = \frac{h}{2}(f(a) + f(b)) + h \sum_{m=1}^{N-1} f(a+mh) - \sum_{j=1}^k h^{2j} \frac{b_{2j}}{(2j)!} (f^{(2j-1)}(b) - f^{(2j-1)}(a)) - h^{2k+2}(b-a) \frac{b_{2k+2}}{(2k+2)!} f^{(2k+2)}(\xi). \quad (36) $$

Here, $b_{2j}, 1 \le j \le k + 1$, are Bernoulli numbers ($b_2 = \frac{1}{6}$, $b_4 = -\frac{1}{30}$, $b_6 = \frac{1}{42}$, $b_8 = -\frac{1}{30}, \dots$). To shorten the notation,

$$ \int_a^b f(x)dx = h \sum_{m=0}^{N} f(a+mh) + \sum_{j=1}^{k} h^{2j} \frac{b_{2j}}{(2j)!} [f^{(2j-1)}(x)]_{x=b}^{x=a} + \mathcal{O}(h^{2k+2}), $$

where the remainder is bounded from above by $h^{2k+2}(b-a) \frac{|b_{2k+2}|}{(2k+2)!} \max_{[a,b]} |f^{(2k+2)}|$. As a corollary of this formula, we estimate the remainder of a bivariate trapezoidal rule as follows.

**Proposition 8** (Two-dimensional Euler-Maclaurin formula). Let $k \ge 0$, $a < b \in \mathbb{R}$, $F \in \mathcal{C}^{2k+2}([a, b]^2)$, $N \ge 1$ and $h = \frac{b-a}{N}$. Then,

$$ \iint_{[a,b]^2} F(x,y)dxdy = h^2 \sum_{m=0}^{N} \sum_{n=0}^{N} F(a+mh, a+nh) + \sum_{j=1}^{k} h^{2j} \frac{b_{2j}}{(2j)!} \int_{\partial[a,b]^2} \partial_{\nu}^{2j-1} F(x,y)d\ell \\ - \sum_{\substack{1 \le i,j \le k \\ i+j \le k}} h^{2(i+j)} \frac{b_{2i}}{(2i)!} \frac{b_{2j}}{(2j)!} \sum_{(x,y) \in \{a,b\}^2} \epsilon(x,y) \partial_x^{2i-1} \partial_y^{2j-1} F(x,y) + \mathcal{O}(h^{2k+2}), \quad (37) $$

where $\ell$ denotes the length measure on $\partial[a, b]^2$, $\nu$ denotes the interior unit normal vector to $\partial[a, b]^2$, and $\epsilon(a, a) = \epsilon(b, b) = 1$, $\epsilon(a, b) = \epsilon(b, a) = -1$. More precisely, the remainder is bounded from above by

$$ h^{2k+2} C \max_{2k+1 \le i+j \le 2k+2} \| \partial_x^i \partial_y^j F(x, y) \|_{\infty}, $$

for some constant $C > 0$ depending only on $b-a$ and $k$.

Note that on the right hand-side in (37), the first term coincides with the one from a bivariate trapezoidal rule on a uniform cartesian grid on $[a, b]^2$, with nodes $(a+mh, a+nh)$, $0 \le m, n \le N$. The other terms relate the error of this quadrature rule with normal derivatives on the boundary $\partial[a, b]^2$, and cross derivatives on the corners $(x, y) \in \{a, b\}^2$.
*Proof of Proposition 8.* Fix $y \in [a, b]$, and apply the relation (36) with $f = F(\cdot, y)$ (and $\xi$ depending on $y$),

$$ \int_a^b F(x, y)dx = h \sum_{m=0}^N' F(a + mh, y) + \sum_{j=1}^k h^{2j} \frac{b_{2j}}{(2j)!} [\partial_x^{(2j-1)} F(x, y)]_{x=b}^{x=a} + \mathcal{O}(h^{2k+2}); $$

here the remainder term is dominated by $h^{2k+2}(b-a) \frac{|b_{2k+2}|}{(2k+2)!} \|\partial_x^{(2k+2)} F\|_{\infty}$, where for any $G \in \mathcal{C}^0([a, b]^2)$, $\|G\|_{\infty} = \max_{[a,b]^2} |G|$. Integrate the relation above for $y \in [a, b]$, and apply again the formula (36), for each integral $\int_a^b F(a + mh, y)dy$ ($0 \le m \le N$); we get

$$ \iint_{[a,b]^2} F(x,y)dxdy = h^2 \sum_{m=0}^{N} \sum_{n=0}^{N}' F(a+mh, a+nh) \\ + \sum_{j=1}^{k} h^{2j} \frac{b_{2j}}{(2j)!} \left( \int_{a}^{b} [\partial_{x}^{2j-1} F(x,y)]_{x=b}^{x=a} dy + h \sum_{m=0}^{N}' [\partial_{y}^{2j-1} F(a+mh,y)]_{y=b}^{y=a} \right) + \mathcal{O}(h^{2k+2}), \quad (38) $$

where the remainder is dominated by $h^{2k+2}(b-a)^2 \frac{|b_{2k+2}|}{(2k+2)!} (\|\partial_x^{2k+2} F\|_{\infty} + \|\partial_y^{2k+2} F\|_{\infty})$. Next, apply again (36) with $f(x) = [\partial_y^{2j-1} F(x,y)]_{y=b}^{y=a} \in \mathcal{C}^{2(k-j+1)}([a,b])$ and a remainder $\mathcal{O}(h^{2(k-j+1)})$,

$$ h \sum_{m=0}^{N}' [\partial_y^{2j-1} F(a+mh, y)]_{y=b}^{y=a} = \int_a^b [\partial_y^{2j-1} F(x, y)]_{y=b}^{y=a} dx - \sum_{i=1}^{k-j} h^{2i} \frac{b_{2i}}{(2i)!} [\partial_x^{2i-1} [\partial_y^{2j-1} F(x, y)]_{y=a}^{y=b}]_{x=a}^{x=b} \\ + \mathcal{O}(h^{2(k-j+1)}), $$

where the remainder is dominated by $h^{2(k-j+1)}(b-a)^2 \frac{|b_{2(k-j+1)}|}{(2(k-j+1))!} \|\partial_x^{2(k-j+1)} \partial_y^{2j-1} F\|_{\infty}$. Inserting this relation into (38) yields

$$ \iint_{[a,b]^2} F(x,y)dxdy = h^2 \sum_{m=0}^{N} \sum_{n=0}^{N}' F(a+mh, a+nh) \\ + \sum_{j=1}^{k} h^{2j} \frac{b_{2j}}{(2j)!} \left( \int_{a}^{b} [\partial_{x}^{2j-1} F(x,y)]_{x=b}^{x=a} dy + \int_{a}^{b} [\partial_{y}^{2j-1} F(x,y)]_{y=b}^{y=a} dx \right) \\ - \sum_{\substack{1 \le i,j \le k \\ i+j \le k}} h^{2(i+j)} \frac{b_{2i}}{(2i)!} \frac{b_{2j}}{(2j)!} \left[ [\partial_{x}^{2i-1} \partial_{y}^{2j-1} F(x,y)]_{y=b}^{y=a} \right]_{x=b}^{x=a} + \mathcal{O}(h^{2k+2}); $$

here, there are two constants $c_1 > 0$ and $c_2 > 0$ depending only $b-a$ and $k$ such that the remainder is dominated by

$$ h^{2k+2} \left( c_1 \|\partial_x^{2k+2} F\|_{\infty} + c_2 \|\partial_y^{2k+2} F\|_{\infty} + c_2 \max_{1 \le j \le k} \|\partial_x^{2(k-j+1)} \partial_y^{2j-1} F\|_{\infty} \right), $$

which concludes the proof. $\square$

<figure><img src="image_5.png" /></figure>

## APPENDIX B. SYMBOLIC COMPUTATION

In this section, we provide a Maple code which checks that $r_4(f)$ in (29) coincides with the given value in (30), for any multivariate polynomial $f(x) = x^d$, $d \in \mathbb{N}^3$.

restart;

\#Fix mu, depending on t=(t1,t2,t3), with ti=tan(alphai), i=1..3

mu:=(1+t1**2)*(1+t2**2)*(1+t3**2)/(2*sqrt(t1**2+t2**2+t3**2)**3);

\# Consider nonnegative integers, d1, d2, d3

assume(d1,int,d1>=0,d2,int,d2>=0,d3,int,d3>=0):

about(d1):about(d2):about(d3):
```c
# Fix f(x)=x^d=x1^d1 x2^d2 x3^d3, with x=t/|t|, (ti=tan(alphai), i=1..3)
# and add the weight mu
f:=t1**d1*t2**d2*t3**d3/sqrt(t1**2+t2**2+t3**2)**(d1+d2+d3);
fmu:=f*mu;
```

```c
#Compute the third derivatives d^3(fmu)/dalphai^3, with i=1..3
#Using the chain rule, with ti=tan(alphai)
#       dg/dalphai=dg/dti dti/dalphai=dg/dti (1+ti^2)
d3alpha1fmu:=simplify(diff(diff(diff(fmu,t1)
                                *(1+t1**2),t1)*(1+t1**2),t1)*(1+t1**2)):
d3alpha2fmu:=simplify(diff(diff(diff(fmu,t2)
                                *(1+t2**2),t2)*(1+t2**2),t2)*(1+t2**2)):
d3alpha3fmu:=simplify(diff(diff(diff(fmu,t3)
                                *(1+t3**2),t3)*(1+t3**2),t3)*(1+t3**2)):
```

```c
#Compute the integral (in r4f) with the third derivatives on the twelve edges,
#       IE=int_E (dnu1^3+dnu2^3)fmu
#Using the change of variable ti=tan(alphai),
#       int_{alphai=-pi/4..pi/4} g(alphai) dalphai
ti := int_{ti=-1..1} g(alphai)/(1+ti^2) dti
IE:=simplify(
    int((-subs(t2=1,subs(t1=1,d3alpha1fmu))-subs(t2=1,subs(t1=1,
                                                    d3alpha2fmu)))/(1+t3**2),t3=-1..1)
    +int((-subs(t2=-1,subs(t1=1,d3alpha2fmu))-subs(t2=-1,subs(t1=1,
                                                    d3alpha1fmu)))/(1+t3**2),t3=-1..1)
    +int((-subs(t3=1,subs(t1=1,d3alpha3fmu))-subs(t3=1,subs(t1=1,
                                                    d3alpha1fmu)))/(1+t2**2),t2=-1..1)
    +int((-subs(t3=-1,subs(t1=1,d3alpha3fmu))-subs(t3=-1,subs(t1=1,
                                                    d3alpha1fmu)))/(1+t2**2),t2=-1..1)
    +int((-subs(t2=1,subs(t1=-1,d3alpha1fmu))-subs(t2=1,subs(t1=-1,
                                                    d3alpha2fmu)))/(1+t3**2),t3=-1..1)
    +int((-subs(t2=-1,subs(t1=-1,d3alpha2fmu))+subs(t2=-1,subs(t1=-1,
                                                    d3alpha1fmu)))/(1+t3**2),t3=-1..1)
    +int((-subs(t3=1,subs(t1=-1,d3alpha3fmu))+subs(t3=1,subs(t1=-1,
                                                    d3alpha1fmu)))/(1+t2**2),t2=-1..1)
    +int((-subs(t3=-1,subs(t1=-1,d3alpha3fmu))+subs(t3=-1,subs(t1=-1,
                                                    d3alpha1fmu)))/(1+t2**2),t2=-1..1)
    +int((-subs(t2=1,subs(t3=1,d3alpha3fmu))-subs(t2=1,subs(t3=1,
                                                    d3alpha2fmu)))/(1+t1**2),t1=-1..1)
    +int((-subs(t2=1,subs(t3=-1,d3alpha3fmu))-subs(t2=1,subs(t3=-1,
                                                    d3alpha2fmu)))/(1+t1**2),t1=-1..1)
    +int((-subs(t2=-1,subs(t3=1,d3alpha3fmu))+subs(t2=-1,subs(t3=1,
                                                    d3alpha2fmu)))/(1+t1**2),t1=-1..1)
    +int((-subs(t2=-1,subs(t3=-1,d3alpha3fmu))+subs(t2=-1,subs(t3=-1,
                                                    d3alpha2fmu)))/(1+t1**2),t1=-1..1));
```

```c
#Compute the cross-derivatives d_alphai d_alphaj (fmu), i neq j, in 1..3
#       (with the chain rule d/dalphai=d/dti (1+ti^2))
der12:=diff(diff(fmu,t1)*(1+t1**2),t2)*(1+t2**2):
der13:=diff(diff(fmu,t1)*(1+t1**2),t3)*(1+t3**2):
der23:=diff(diff(fmu,t2)*(1+t2**2),t3)*(1+t3**2):

#Compute the term with cross-derivatives on the corners (in r4f)
```
```matlab
# (with ti=-1,1 for alphai=-pi/4, pi/4)
SV:=simplify(
    subs(t1=1,subs(t2=1,subs(t3=1,t1*t2*der12+t1*t3*der13+t2*t3*der23)))
    +subs(t1=-1,subs(t2=1,subs(t3=1,t1*t2*der12+t1*t3*der13+t2*t3*der23)))
    +subs(t1=1,subs(t2=-1,subs(t3=1,t1*t2*der12+t1*t3*der13+t2*t3*der23)))
    +subs(t1=1,subs(t2=1,subs(t3=-1,t1*t2*der12+t1*t3*der13+t2*t3*der23)))
    +subs(t1=-1,subs(t2=-1,subs(t3=1,t1*t2*der12+t1*t3*der13+t2*t3*der23)))
    +subs(t1=-1,subs(t2=1,subs(t3=-1,t1*t2*der12+t1*t3*der13+t2*t3*der23)))
    +subs(t1=1,subs(t2=-1,subs(t3=-1,t1*t2*der12+t1*t3*der13+t2*t3*der23)))
    +subs(t1=-1,subs(t2=-1,subs(t3=-1,t1*t2*der12+t1*t3*der13+t2*t3*der23)));
```

```matlab
#Compute the fourth order term r4f
c4:=-1/(30*2*3*4):c2:=1/12:r4f:=simplify(c4*IE-c2**2*SV);
```

```matlab
#Compute the proposed value for r4f
# (using the weight omega=mu(1,1,1)=4/(3^3/2) on any vertex)
omega:=4/(3**(3/2)):
```

```matlab
propValue:=-1/15*omega*simplify(
    subs(t1=1,subs(t2=1,subs(t3=1,f)))
    +subs(t1=-1,subs(t2=1,subs(t3=1,f)))
    +subs(t1=1,subs(t2=-1,subs(t3=1,f)))
    +subs(t1=1,subs(t2=1,subs(t3=-1,f)))
    +subs(t1=-1,subs(t2=-1,subs(t3=1,f)))
    +subs(t1=-1,subs(t2=1,subs(t3=-1,f)))
    +subs(t1=1,subs(t2=-1,subs(t3=-1,f)))
    +subs(t1=-1,subs(t2=-1,subs(t3=-1,f))));
```

```matlab
#Comparison of r4f and the proposed value
simplify(r4f-propValue);
```

## APPENDIX C. MATLAB CODE

The following Matlab function implements the Cubed Sphere CS<sub>N</sub> and the (corrected) quad-rature weight of the sixth-order rule (T-CS6).

```matlab
function [xyz,w]=CSw(N)
% [xyz,w]=CSw(N)
% Input:
% N>=1, parameter of the Cubed Sphere CS_N
% Output:
% xyz, 6*N^2+2 x 3, Cubed Sphere CS_N
% w, 6*N^2+2 x 1, weight of a 6th correction of the trapezoidal rule

%Cubed Sphere in gnomonic coordinates
t=tan(linspace(-pi/4,pi/4,N+1));
xyz=[ones((N+1)^2,1), kron(t,ones(N+1,1)),... 
    kron(ones(N+1,1),t(end:-1:1));...%Face I, with edges
    kron(t(end-1:-1:1),ones(N+1,1)), ones(N*(N+1),1),... 
    kron(ones(N,1),t(end:-1:1));...%Face II, without the edges in I
    -ones(N*(N+1),1), kron(t(end-1:-1:1),ones(N+1,1)),... 
    kron(ones(N,1),t(end:-1:1));...%Face III, without the edges II
    kron(t(2:end-1),ones(N+1,1)), -ones((N-1)*(N+1),1),... 
    kron(ones(N-1,1),t(end:-1:1));...%Face IV, without the edges in I,III
    kron(ones(N-1,1),t(2:end-1)), kron(t(2:end-1),ones(N-1,1)),... 
    ones((N-1)^2,1);...%interior of face V
```
kron(ones(N-1,1),t(end-1:-1:2)), kron(t(2:end-1),ones(N-1,1)),... -ones((N-1)^2,1)];%interior of face VI nxyz=sqrt(sum(xyz.^2,2));%norm %corrected trapezoidal weight h=pi/(2*N);w=h^2/2*prod(1+xyz.^2,2)/nxyz.^3; iV=[1 N+1 N*(N+1)+1 (N+1)^2 2*N*(N+1)+1 (N+1)*(2*N+1)... 3*N*(N+1)+1 (N+1)*(3*N+1)];%indices of the vertices w(iV)=(1-h^2/15)*w(iV); %radial projection xyz=xyz./nxyz; end

REFERENCES

[1] C. Ahrens and G. Beylkin. Rotationally invariant quadratures for the sphere. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 465(2110):3103–3125, 2009.

[2] C. An, X. Chen, I. H. Sloan, and R. S. Womersley. Well Conditioned Spherical Designs for Integration and Interpolation on the Two-Sphere. *SIAM journal on numerical analysis*, 48(6):2135–2157, 2010.

[3] C. An and H.-N. Wu. Spherical Configurations and Quadrature Methods for Integral Equations of the Second Kind. *SIAM Journal on Numerical Analysis*, 64(1):148–169, 2026.

[4] T. C. Andrews and C. Jablonowski. Stability analyses of divergence and vorticity damping on gnomonic cubed-sphere grids. Working document https://arxiv.org/abs/2505.05624, 2025.

[5] C. H. Beentjes. Quadrature on a spherical surface. Technical report, Oxford University - https://cbeentjes.github.io/notes/2015-Quadrature-Sphere, 2015.

[6] J.-B. Bellet. Symmetry group of the equiangular cubed sphere. *Quarterly of Applied Mathematics*, 80:69–86, 2022.

[7] J.-B. Bellet. A discrete Funk transform on the Cubed Sphere. *Journal of Computational and Applied Mathematics*, 429(115205), 2023.

[8] J.-B. Bellet. Mathematical and numerical methods for three-dimensional reflective tomography and for approximation on the sphere. Habilitation thesis, Université de Lorraine, 2023.

[9] J.-B. Bellet, M. Brachet, and J.-P. Croisille. Quadrature and symmetry on the Cubed Sphere. *Journal of Computational and Applied Mathematics*, 409(114142), 2022. Numerical tables available at https://hal.science/hal-03223150v1/file/xyzwCSN.zip.

[10] J.-B. Bellet, M. Brachet, and J.-P. Croisille. Interpolation on the Cubed Sphere with Spherical Harmonics. *Numerische Mathematik*, 153:249–278, 2023.

[11] J.-B. Bellet, M. Brachet, and J.-P. Croisille. Quadrature on the cubed sphere: the low resolution case. Working document https://hal.science/hal-04807672v1, 2024.

[12] J.-B. Bellet and J.-P. Croisille. Least Squares Spherical Harmonics Approximation on the Cubed Sphere. *Journal of Computational and Applied Mathematics*, 429(115213), 2023.

[13] L. Bourhrara, T. Collignon, B. Aylaj, and J.-F. Vidal. On the Integrals of Products of Multiple Real Spherical Harmonics over the Full Sphere and Hemispheres. *Journal of Scientific Computing*, 106(11), 2026.

[14] M. Brachet. *Schémas compacts hermitiens sur la Sphère: applications en climatologie et océanographie numérique*. PhD thesis, Université de Lorraine, 2018 (in French).

[15] M. Brachet and J.-P. Croisille. Spherical Shallow Water simulation by a cubed sphere finite difference solver. *Quarterly Journal of the Royal Meteorological Society*, 147(735):786–800, 2021.

[16] A. Chen and K. Robert. A Cubed Sphere Fast Multipole Method. Working document https://arxiv.org/abs/2508.13550, 2025.

[17] C. Chen and F. Xiao. Shallow water model on cubed-sphere by multi-moment finite volume method. *Journal of Computational Physics*, 227(10):5019–5044, 2008.

[18] S. Chevrot, R. Martin, and D. Komatitsch. Optimized discrete wavelet transforms in the cubed sphere with the lifting scheme-implications for global finite-frequency tomography. *Geophysical Journal International*, 191(3):1391–1402, 2012.

[19] J.-P. Croisille. Hermitian Compact Interpolation on the Cubed-Sphere Grid. *Journal of Scientific Computing*, 57:193–212, 2013.

[20] J.-P. Croisille. Hermitian approximation of the spherical divergence on the Cubed-Sphere. *Journal of Computational and Applied Mathematics*, 280:188–201, 2015.

[21] M. Faham and H. Nasir. Weakly Orthogonal Spherical Harmonics in a Non-Polar Spherical Coordinates and its Application to Functions on Cubed-Sphere. *Sultan Qaboos University Journal for Science*, 17(2):200–213, 2012.

[22] J. Fliege and U. Maier. The distribution of points on the sphere and corresponding cubature formulae. *IMA Journal of Numerical Analysis*, 19(2):317–334, 1999.
[23] B. Fornberg and J. M. Martel. On spherical harmonics based numerical quadrature over the surface of a sphere. *Advances in Computational Mathematics*, 40(5-6):1169–1184, 2014.

[24] K. Hesse, I. H. Sloan, and R. S. Womersley. Numerical integration on the sphere. In W. Freeden, Z. M. Nashed, and T. Sonar, editors, *Handbook of Geomathematics*. Springer, 2010.

[25] L. Ivan, H. De Sterck, S. A. Northrup, and C. P. T. Groth. Multi-dimensional finite-volume scheme for hyperbolic conservation laws on three-dimensional solution-adaptive cubed-sphere grids. *Journal of Computational Physics*, 255:205–227, 2013.

[26] B. A. Jones, G. H. Born, and G. Beylkin. Comparison of the Cubed-Sphere Gravity Model with the Spherical Harmonics. *Journal of Guidance, Control, and Dynamics*, 33(2):415–425, 2010.

[27] H.-G. Kang and H.-B. Cheong. An efficient implementation of a high-order filter for a cubed-sphere spectral element model. *Journal of Computational Physics*, 332:66–82, 2017.

[28] V. I. Lebedev. Quadratures on a sphere. *USSR Computational Mathematics and Mathematical Physics*, 16(2):10–24, 1976.

[29] D. Lee and A. Palha. A mixed mimetic spectral element model of the rotating shallow water equations on the cubed sphere. *Journal of Computational Physics*, 375:240–262, 2018.

[30] J. L. McGregor. Semi-Lagrangian Advection on Conformal-Cubic Grids. *Monthly Weather Review*, 124(6):1311–1322, 1996.

[31] R. D. Nair, S. J. Thomas, and R. D. Loft. A Discontinuous Galerkin Transport Scheme on the Cubed Sphere. *Monthly Weather Review*, 133(4):814–828, 2005.

[32] B. Portelenelle and J.-P. Croisille. An efficient quadrature rule on the Cubed Sphere. *Journal of Computational and Applied Mathematics*, 328:59–74, 2018.

[33] R. J. Purser and M. Tong. A minor modification of the gnomonic cubed-shaped sphere grid that offers advantages in the context of implementing moving hurricane nests. *Office note, National Centers for Environmental Prediction (U.S.)*, 486, 2017.

[34] W. M. Putman. *Development of the finite-volume dynamical core on the cubed-sphere*. PhD thesis, The Florida State University, 2007.

[35] M. Rančić, R. J. Purser, and F. Mesinger. A global shallow-water model using an expanded spherical cube: Gnomonic versus conformal coordinates. *Quarterly Journal of the Royal Meteorological Society*, 122(532):959–982, 1996.

[36] M. Rančić, R. J. Purser, D. Jović, R. Vasic, and T. Black. A Nonhydrostatic Multiscale Model on the Uniform Jacobian Cubed Sphere. *Monthly Weather Review*, 145(3):1083–1105, 2017.

[37] C. Ronchi, R. Iacono, and P. S. Paolucci. The “cubed sphere”: a new method for the solution of partial differential equations in spherical geometry. *Journal of Computational Physics*, 124(1):93–114, 1996.

[38] J. A. Rossmanith. A wave propagation method for hyperbolic systems on the sphere. *Journal of Computational Physics*, 213(2):629–658, 2006.

[39] R. Sadourny. Conservative finite-difference approximations of the primitive equations on quasi-uniform spherical grids. *Monthly Weather Review*, 100(2):136–144, 1972.

[40] S. L. Sobolev. Cubature formulas on the sphere invariant under finite groups of rotations. *Doklady Akademii Nauk SSSR*, 146(2):310–313, 1962.

[41] S. J. Thomas, J. M. Dennis, H. M. Tufo, and P. F. Fischer. A Schwarz preconditioner for the cubed-sphere. *SIAM Journal on Scientific Computing*, 25(2):442–453, 2003.

[42] L. N. Trefethen and J. A. C. Weideman. The Exponentially Convergent Trapezoidal Rule. *SIAM Review*, 56(3):385–458, 2014.

[43] P. A. Ullrich, P. H. Lauritzen, and C. Jablonowski. Some considerations for high-order 'incremental remap'-based transport schemes: edges, reconstructions, and area integration. *International Journal for Numerical Methods in Fluids*, 71(9):1131–1151, 2013.

[44] R. S. Womersley. Efficient spherical designs with good geometric properties. In *Contemporary computational mathematics-A celebration of the 80th birthday of Ian Sloan*, pages 1243–1285. Springer, 2018.

[45] M. Zerroukat and T. Allen. On the corners of the cubed-sphere grid. *Quarterly Journal of the Royal Meteorological Society*, 148(743):778–783, 2022.

<footnote>† UNIVERSITÉ DE POITIERS, CNRS, LMA, POITIERS, FRANCE</footnote>

<footnote>Email address: jean.baptiste.bellet@univ-poitiers.fr</footnote>