Smooth Polar B-Splines with High-Order Regularity at the Origin

TL;DR

This work tackles axis regularity in polar-coordinate discretizations on the unit disk by introducing smooth polar splines that enforce high-order regularity at the polar origin without altering the underlying geometry. The key idea is a Galerkin $L^2$ projection of harmonic polar functions $S_l^{\pm m}(r,\theta)=r^l h_m(\theta)$ onto a center subspace spanned by a few radial B-splines and angular harmonics, yielding center-spline functions $\widetilde{B}_l^m$ that preserve $r^l$ for $l\le p$ and decouple by angular mode $|m|$, with $C^\infty$-regularity recovered asymptotically as $Δ\theta\to0$. The resulting smooth polar-spline space $\widetilde{V}_h$ embeds into the original tensor-product B-spline space via an exact prolongation operator, enabling unchanged force/potential evaluations while solving the Galerkin system in the regularized subspace. Numerically, the method improves matrix conditioning, suppresses spurious eigenvalues, and reduces PIC-driven statistical noise near the axis, demonstrated in Poisson, Helmholtz, and TAE contexts, and implemented within a gyrokinetic PIC setting (EUTERPE). Overall, the approach provides a robust, high-order, and efficient adaptation of tensor-product B-splines to polar coordinates with rigorous axis regularity control and practical benefits for plasma-physics simulations.

Abstract

We introduce a smooth B-spline discretization in polar coordinates on the unit disc that corrects the loss of regularity present at the origin caused by the coordinate singularity in standard tensor-product B-spline formulations. The method constructs "smooth polar splines" via a Galerkin projection of harmonic polar functions $S_l^{-m}(r,θ) := r^l \sin(mθ)$ and $S_l^{m}(r,θ) := r^l \cos(mθ)$, derived from the polar representation of Cartesian monomials, onto the central tensor-product B-spline basis in the innermost radial region. The radial component reproduces $r^l$ exactly for $0 \le l \leq p$, where $p$ is the B-spline degree, satisfying the near-origin regularity condition. However, exact compatibility with $C^\infty$-regularity at the origin is recovered only in the limit $Δθ\to 0$, when the angular component resolves all angular harmonics accurately. The smooth polar splines are linear combinations of standard tensor-product B-splines and lie in the same function space, enabling mapping between the $C^\infty$-regular subspace and the original discretization space via an exact prolongation operator and a corresponding restriction operator acting on the discrete variables. They match standard tensor-product B-splines away from the origin, preserve orthogonality among the newly constructed origin-centered basis functions, and maintain local support and sparse matrices. This smoothness and locality improve the conditioning of mass and stiffness matrices, conserve charge, and reduce statistical errors in particle-in-cell simulations near the origin, while eliminating spurious eigenvalues in eigenvalue problems. The approach provides a robust, high-order, and efficient adaptation of tensor-product B-splines for polar coordinates in physics simulations.

Figures (28)

• Figure 1: Linear, quadratic, and cubic B-spline basis functions at the left boundary.
• Figure 2: Quadratic and cubic periodic B-spline basis functions near $\theta=0$.
• Figure 3: Radial parts $\widetilde{R}_l$ of the center-spline basis functions (colored lines) expressed as linear combinations of cubic B-spline basis functions, reproducing the monomials in the innermost radial interval $\tilde{r} \in [0,1]$. Dashed lines indicate the monomials extended to $\tilde{r} \in [1,4]$.
• Figure 4: Coefficient spaces $\mathbb{R}^N$ of the tensor-product B-splines (green) and $\mathbb{R}^{\widetilde{N}}$ of the smooth polar splines (blue). The $C^\infty$-regular subspace $\mathcal{R}(\mathit{\bm{\mathsf{P}}}) \subset \mathbb{R}^N$ (yellow) contains coefficient vectors $\mathit{\bm{\mathsf{\phi}}}_{\mathrm{reg}}$ that are $C^\infty$-regular at the polar origin. The $L^2$-orthogonal projection $\widetilde{\mathit{\bm{\mathsf{M}}}}^{-1} \mathit{\bm{\mathsf{P}}}^{\mathrm T} \mathit{\bm{\mathsf{M}}}$ (blue arrow) maps $\mathit{\bm{\mathsf{\phi}}} \in \mathbb{R}^N$ onto $\widetilde{\mathit{\bm{\mathsf{\phi}}}} \in \mathbb{R}^{\widetilde{N}}$ in the smooth polar-spline space. The prolongation $\mathit{\bm{\mathsf{P}}}$ (red arrow) maps $\widetilde{\mathit{\bm{\mathsf{\phi}}}}$ into $\mathit{\bm{\mathsf{\phi}}}_{\mathrm{reg}} \in \mathcal{R}(\mathit{\bm{\mathsf{P}}})$. The projector $\mathit{\bm{\mathsf{\Pi}}}_\perp$ given by Eq. \ref{['eq.Pi_perp_operator']} (green arrow) maps any coefficient vector $\mathit{\bm{\mathsf{\phi}}} \in \mathbb{R}^N$ to its $C^\infty$-regular (at the polar origin) part $\mathit{\bm{\mathsf{\phi}}}_{\mathrm{reg}}$.
• Figure 5: The load vector $\boldsymbol{f}$ is calculated in the original tensor-product B-spline basis and then transformed to the smooth polar-spline basis, in which the matrix equation is solved. The resulting coefficient vector $\widetilde{\boldsymbol{\phi}}$ is subsequently mapped back to the original basis.