Series Of Spherical Harmonics Research Articles

Divergence beneath the Brillouin sphere and the phenomenology of prediction error in spherical harmonic series approximations of the gravitational field

The Brillouin sphere is defined as the smallest sphere, centered at the origin of the geocentric coordinate system, that incorporates all the condensed matter composing the planet. The Brillouin sphere touches the Earth at a single point, and the radial line that begins at the origin and passes through that point is called the singular radial line. For about 60 years there has been a persistent anxiety about whether or not a spherical harmonic (SH) expansion of the external gravitational potential, V, will converge beneath the Brillouin sphere. Recently, it was proven that the probability of such convergence is zero. One of these proofs provided an asymptotic relation, called Costin’s formula, for the upper bound, EN , on the absolute value of the prediction error, eN , of a SH series model, VN(θ,λ,r) , truncated at some maximum degree, N=nmax . When the SH series is restricted to (or projected onto) a particular radial line, it reduces to a Taylor series (TS) in 1/r . Costin’s formula is EN≃BN−b(R/r)N , where R is the radius of the Brillouin sphere. This formula depends on two positive parameters: b, which controls the decay of error amplitude as a function of N when r is fixed, and a scale factor B. We show here that Costin’s formula derives from a similar asymptotic relation for the upper bound, An on the absolute value of the TS coefficients, an , for the same radial line. This formula, An≃Kn−k , depends on degree, n, and two positive parameters, k and K, that are analogous to b and B. We use synthetic planets, for which we can compute the potential, V, and also the radial component of gravitational acceleration, gr=∂V/∂r , to hundreds of significant digits, to validate both of these asymptotic formulas. Let superscript V refer to asymptotic parameters associated with the coefficients and prediction errors for gravitational potential, and superscript g to the coefficients and predictions errors associated with gr . For polyhedral planets of uniform density we show that bV=kV=7/2 and bg=kg=5/2 almost everywhere. We show that the frequency of oscillation (around zero) of the TS coefficients and the series prediction errors, for a given radial line, is controlled by the geocentric angle, α, between that radial line and the singular radial line. We also derive useful identities connecting KV,BV,Kg , and Bg . These identities are expressed in terms of quotients of the various scale factors. The only other quantities involved in these identities are α and R. The phenomenology of ‘series divergence’ and prediction error (when r < R) can be described as a function of the truncation degree, N, or the depth, d, beneath the Brillouin sphere. For a fixed r⩽R , as N increases from very low values, the upper error bound EN shrinks until it reaches its minimum (best) value when N reaches some particular or optimum value, Nopt . When N>Nopt , prediction error grows as N continues to increase. Eventually, when N≫Nopt , prediction errors increase exponentially with rising N. If we fix the value of N and allow R/r to vary, then we find that prediction error in free space beneath the Brillouin sphere increases exponentially with depth, d, beneath the Brillouin sphere. Because bg=bV−1 everywhere, divergence driven prediction error intensifies more rapidly for gr than for V, both in terms of its dependence on N and d. If we fix both N and d, and focus on the ‘lateral’ variations in prediction error, we observe that divergence and prediction error tend to increase (as does B) as we approach high-amplitude topography.

Optimizing the parameters of spherical harmonic reconstruction of the 3D morphology of manufactured aggregate with different particle sizes

The optimal reconstruction parameters of spherical harmonic (SH) series for the aggregates with different morphologies and sizes were not clear. In this paper, manufactured aggregates with different sizes (1.18–2.36 mm, 2.36–4.75 mm, 4.75–9.50 mm), and natural aggregate with size >1.18 mm were scanned using the digital light projection (DLP), and reconstructed using SH series. The effect of maximum number of SH functions (N) used in the reconstruction series and the number of terms in the Gaussian quadrature (G) used for solving the SH expansion coefficients on the reconstruction errors were studied. The optimal reconstruction parameters for different particle sizes and morphologies were proposed, via a derived relationship between N and G. The reconstruction errors decreased when the N and G values increased. The relationship between N and G was derived to be G ≥ 1.84 N + 3.60, which is probably caused by the issue of the algebraic precision of the Gaussian integral.

Breaking down walls: Continuous potential models for internal motions in NMR spin relaxation

Simple physical models for restricted diffusion in a potential, which provide important insights for NMR spin relaxation, usually are based on free diffusion within rigid boundaries or diffusion in relatively simple continuous potential energy surfaces. The diffusion-in-a-cone model is an example of the former and diffusion in an N-fold cosine potential is an example of the latter. The present work models restricted diffusion for arbitrary potential energy functions on the surface of a cone or a sphere, by expanding the potentials in Fourier or spherical harmonic series, respectively. The results exhibit simple relationships between generalized order parameters and effective correlation times, critical for analysis of experimental spin relaxation data, and illustrate the transition from diffusive-like to jump-like behavior in multi-well potentials.

A G-Modified Helmholtz Equation with New Expansions for the Earth’s Disturbing Gravitational Potential, Its Functionals and the Study of Isogravitational Surfaces

The G-modified Helmholtz equation is a partial differential equation that enables us to express gravity intensity g as a series of spherical harmonics having radial distance r in irrational powers. The Laplace equation in three-dimensional space (in Cartesian coordinates, is the sum of the second-order partial derivatives of the unknown quantity equal to zero) is used to express the Earth’s gravity potential (disturbing and normal potential) in order to represent other useful quantities—which are also known as functionals of the disturbing potential—such as gravity disturbance, gravity anomaly, and geoid undulation as a series of spherical harmonics. We demonstrate that by using the G-modified Helmholtz equation, not only gravity intensity but also disturbing potential and its functionals can be expressed as a series of spherical harmonics. Having gravity intensity represented as a series of spherical harmonics allows us to create new Global Gravity Models. Furthermore, a more detailed examination of the Earth’s isogravitational surfaces is conducted. Finally, we tabulate our results, which makes it clear that new Global Gravity Models for gravity intensity g will be very useful for many geophysical and geodetic applications.

Efficient computation of the geopotential gradient in graphic processing units

Efficient computation of the geopotential gradient is essential for numerical propagators, particularly in scenarios involving low Earth orbits. Conventional geopotential calculations are based on spherical harmonics series, which become computationally demanding as the degree/order increases. This computational burden can be mitigated by means of parallelized algorithms. Additionally, certain situations lend themselves to high parallelization, such as the propagation of space debris catalogs, satellite mega-constellations, or the dispersion of particles resulting from a space collision event. This paper introduces an optimized Graphics Processing Unit (GPU) implementation designed to facilitate extensive parallelization in the geopotential gradient calculation. The formulation developed in this study is not specific to any GPU. However, to illustrate the low-level optimizations necessary for an efficient implementation, we selected the Compute Unified Device Architecture (CUDA) as the dominant and de facto standard in parallel computing. Nevertheless, most of the concepts and optimizations presented in this paper are also valid for other GPU architectures. Built upon the spherical harmonic expansion using the Cunningham formulation, which is well-suited for GPU computations, our implementation offers several variants with different tradeoffs between speed and accuracy. Besides GPU double precision, we introduced a mixed precision arithmetic –a hybrid between single and double precision– that exploits GPU capabilities with a low penalty in accuracy. The proposed algorithm was implemented as a software reusable module, and its performance was evaluated against GMAT, GODOT, and Orekit astrodynamic codes. The algorithm’s accuracy in double precision is comparable to such codes. The mixed precision version showed enough accuracy for LEO satellite propagation, with around 1 m difference in four days. Testing across different CUDA architectures revealed very high speed-up factors compared to a single CPU, reaching a speed-up of 645 for the mixed precision variant and 450 for the double precision one in the propagation of about 3200 objects with a geopotential of degree/order 126 × 126 using an A100 GPU device.

Using spherical harmonics to solve the Boltzmann equation: an operator-based approach

ABSTRACT The transport of charged particles or photons in a scattering medium can be modelled with a Boltzmann equation. The mathematical treatment for scattering in such scenarios is often simplified if evaluated in a frame where the scattering centres are, on average, at rest. It is common therefore to use a mixed coordinate system, wherein space and time are measured in a fixed inertial frame, while momenta are measured in a ‘co-moving’ frame. To facilitate analytical and numerical solutions, the momentum dependence of the phase-space density may be expanded as a series of spherical harmonics, typically truncated at low order. A method for deriving the system of equations for the expansion coefficients of the spherical harmonics to arbitrary order is presented in the limit of isotropic, small-angle scattering. The method of derivation takes advantage of operators acting on the space of spherical harmonics. The matrix representations of these operators are employed to compute the system of equations. The computation of matrix representations is detailed and subsequently simplified with the aid of rotations of the coordinate system. The eigenvalues and eigenvectors of the matrix representations are investigated to prepare the application of standard numerical techniques, e.g. the finite volume method or the discontinuous Galerkin method, to solve the system.

Theoretical and numerical simulation to control the wideband three-dimensional interfering noise by the broadband constant beam pattern method for a spherical array

The broadband constant beam pattern (CBP) spherical array synthesis theory is applied to suppress or to cancel wideband three-dimensional interfering noise sources by implanting zeros in the array's original synthesized shading function. This modified array angular shading function can be expanded by a series of spherical harmonics that will converge to the beam pattern in the far-field such that the created nulls are in the wideband interfering noise source directions per the CBP theory where the ratio of the spherical array radius to the operating frequency wavelength is large. The simulated numerical examples given for this wideband noise source suppression method demonstrate that a broadband CBP performance is maintained for the spherical array with a Legendre polynomial shading function, the classic Dolph–Chebyshev shading function, or a combination. With the CBP design, one set of the real shading functions works for all frequencies in the array's operating band to cancel or to suppress three-dimensional wideband interfering noise sources.

A New Modified Helmholtz Equation for the Expression of the Gravity Gradient and the Intensity of an Electrostatic Field in Spherical Harmonics

In this work, it is shown that the geometry of a gravity field generated by a spheroid with low eccentricity can be described with the help of a newly modified Helmholtz equation. To distinguish this equation from the modified Helmholtz equation, we refer to it as the G-modified Helmholtz equation. The use of this new equation to study the spheroid’s gravity field is advantageous in expressing the gravity vector as a vector series of spherical harmonics. The solution of the G-modified Helmholtz equation involves both the gravity intensity g (or simply gravity g) and the intensity E of an electrostatic field as shown in sequel. An electrostatic field generated by an oblate spheroid charged with l electrons (uniform ellipsoidal charge distribution) is demonstrated to be a special case. Both gravity intensity g and intensity E are governed by the same law and can be expressed as a series of spherical harmonics, and thus the G-modified Helmholtz equation is useful for describing the gravity and electrostatic fields.

Vlasov-Fokker-Planck-Maxwell simulations for plasmas in inertial confinement fusion

This work presents a code to simulate collisions and laser-plasma interactions of multiple-species plasmas by solving 1D3V Vlasov-Fokker-Planck-Maxwell equations. The distribution functions of all species are expanded into spherical harmonic series of arbitrary order. Each electron or ion species possesses its distinct and suitable velocity mesh and the number of harmonics. Collisions between different ion species are considered and treated with the Fokker-Planck collision term. The code employs the total variation diminishing (TVD) scheme to avoid numerical instabilities in the region of steep density, rarefied plasma-vacuum interface, and non-periodic boundary conditions, which makes it possible to simulate cases that are close to actual physical scenes. The number of particles is globally conserved for each species, the energy for each species is conserved for collisions between the same species, and total energy is conserved for collisions between different species.

Comparison between the geological features of Venus and Earth based on gravity aspects

We probe the gravitational properties of two neighboring planets, Earth and Venus. To justify a comparison between gravity models of the two planets, spherical harmonic series were considered up to a degree and order of 100. The topography and gravity aspects, including {Gamma }_{zz} (vertical derivative of the vertical component of the gravity field), strike alignment (SA), comb factor (CF), and I2 invariant derived from the Marussi tensor, were calculated for the two planets at specifically selected zones that provided sufficient resolution. From Γzz we discovered that the N-NW edge of Lakshmi Planum does not show any subduction-like features. Its Γzz signature resembles passive continental margins on Earth, like those surrounding the Indian Peninsula. Moreover, according to SA and CF, the Pacific and Philippine-North American Contact Zone on Earth indicates significantly higher level of deformation due to convergent motion of the plates, whereas the deformation level on Venus is significantly smaller and local, when considering an equatorial rifting zone (ERZ) of Venus (between Atla-Beta Regios) as diverging boundaries. The strain mode on the East African Rift system is smaller in comparison with ERZ as its Venusian analog. The topography-I2 analysis suggests a complicated nature of the topographic rise on Beta Regio. We show that specific regions in this volcanic rise are in incipient stages of upward motion, with denser mantle material approaching the surface and thinning the crust, whereas some risen districts show molten and less dense underlying crustal materials. Other elevated districts appear to be due to mantle plumes and local volcanic activities with large density of underlying material.

Spherical Harmonics for the 1D Radiative Transfer Equation. I. Reflected Light

A significant challenge in radiative transfer theory for atmospheres of exoplanets and brown dwarfs is the derivation of computationally efficient methods that have adequate fidelity to more precise, numerically demanding solutions. In this work, we extend the capability of the first open-source radiative transfer model for computing the reflected light of exoplanets at any phase geometry, PICASO (Planetary Intensity Code for Atmospheric Spectroscopy Observations). Until now, PICASO has implemented two-stream approaches to the solving the radiative transfer equation for reflected light, in particular following the derivations of Toon et al. In order to improve the model accuracy, we have considered higher-order approximations of the phase functions; namely, we have increased the order of approximation from two to four, using spherical harmonics. The spherical harmonics approximation decouples spatial and directional dependencies by expanding the intensity and phase function into a series of spherical harmonics, or Legendre polynomials, allowing for analytical solutions for low-order approximations to optimize computational efficiency. We rigorously derive the spherical harmonics method for reflected light and benchmark the four-term method (SH4) against Toon et al. and two independent and higher-fidelity methods (CDISORT and doubling method). On average, the SH4 method provides an order-of-magnitude increase in accuracy, compared to Toon et al. Finally, we implement SH4 within PICASO and observe only a modest increase in computational time, compared to two-stream methods (20% increase).

ID: 216730 Comparing Tractography Methods for Preoperative Deep Brain Stimulation Planning in Essential Tremor Patients

Accurate lead placement is essential to successful treatment with deep brain stimulation (DBS). The Ventralis Intermediate Nucleus (VIM) of the thalamus has been identified as a surgical target for DBS neuromodulation of motor input with successful outcomes [1-3]. To aid in this, prospective pre-operative tractography can be implemented for targeting the VIM in ET patients [4,5]. For this, FSL’s ball-and-stick model based probabilistic tracking” (Probtrackx) has been validated to be an effective tractography method (TM) [5]. However, with many other tractography widely available and embedded in commercial surgical platforms, there is benefit in evaluating the relative performance of tractography methods (TMs). Here, we assessed the relative performance of five TMs from the MRTRIX toolbox: three of these techniques were based on fiber orientation distribution functions (ODFs) represented by spherical harmonic series: iFOD1, iFOD2, iFOD-rk4, and two tensor-based models: Tensor-Deterministic, and Tensor-Probabilistic.

The Gravity Force Generated by a Non-Rotating Level Ellipsoid of Revolution with Low Eccentricity as a Series of Spherical Harmonics

The gravity force of a gravity field generated by a non-rotating level ellipsoid of revolution enclosing mass M is given as a solution of a partial differential equation along with a boundary condition of Dirichlet type. The partial differential equation is formulated herein on the basis of the behavior of spherical gravity fields. A classical solution to this equation is represented on the basis of spherical harmonics. The series representation of the solution is exploited in order to conduct a rigorous asymptotic analysis with respect to eccentricity. Finally, the Dirichlet boundary problem is solved for the case of an ellipsoid of revolution (spheroid) with low eccentricity. This has been accomplished on the basis of asymptotic analysis, which resulted in the determination of the coefficients participating in the spherical harmonics expansion. The limiting case of this series expresses the gravity force of a non-rotating sphere.

Effect mechanism of block convexity on the shear behaviors of soil-rock mixtures by the developed 3D spherical harmonics-based modeling approach

The effect of block convexity on the shear behaviors of soil-rock mixtures (S-RMs) can be investigated using accurate 3D models constructed with the discrete element method (DEM). Based on the improved spherical harmonics series, a new block geometry modeling approach was proposed for generating 3D blocks with different convexities but the same form and angularity. By introducing tracer particles, a DEM modeling approach was also proposed for generating S-RM models with the same block distribution characteristics and possible block breakage. Numerical direct shear tests were performed on the generated stiff and soft S-RM models with different block convexities. The meso- and the macro-shear behaviors of the S-RMs were analyzed in detail to reveal the effect mechanism of block convexity. The results showed that the breakage degree of soft blocks increased linearly with the block convexity degree. The interlocking degrees of blocks with higher convexities are higher, leading to higher contributions of the blocks to the shear strength and greater block rotation. As the effect of block breakage is more significant, soft blocks with higher convexity degrees contribute smaller to the shear resistance. For S-RMs with higher block convexity degrees, the shear strength of stiff S-RMs is higher owing to their higher block interlocking degrees, while the shear strength of soft S-RMs under high normal stress is smaller due to their higher block breakage degrees. The dilation magnitude is greater for stiff S-RMs while lower for soft S-RMs with higher convexity degrees, owing to the greater effect of block breakage than the rotation behavior. The friction angles of stiff S-RMs with higher block convexity degrees are larger due to their higher block interlocking degrees, while the friction angles of soft S-RMs with lower block convexity degrees are smaller under high normal stress due to their higher block breakage degrees.

Elastic wave propagation in anisotropic polycrystals: inferring physical properties of glacier ice

An optimization problem is proposed for inferring physical properties of polycrystals given ultrasonic (elastic) wave velocity measurements, made across multiple sample orientations. The feasibility of the method is demonstrated by inferring both the effective grain elastic parameters and the grain c -axis orientation distribution function (ODF) of ice-core samples from Priestley glacier, Antarctica. The method relies on expanding the ODF in terms of a spherical harmonic series, which allows for a non-parametric estimation of the sample ODF. Moreover, any linear combination of the Voigt (strain) and Reuss (stress) homogenization scheme is allowed, although for glacier ice, the exact choice is found to matter little for bulk elastic behaviour, and thus for inferred physical properties, too. Finally, the accuracy of the inferred grain elastic parameters is discussed, including the well-posedness and shortcomings of the inverse problem, relevant for future adoptions in glaciology, geology and elsewhere.

Assessment of recent global geopotential models based on the Auvergne test area data

The gravitational field is important for many natural phenomena related to earth dynamics, especially mass transport. Its precise determination is essential for earth sciences such as geodesy, geophysics and oceanography. Determining the earth’s gravitational field is the same as determining the earth’s potential. Since this potential is a harmonic function outside the earth, spherical harmonic series are often used to represent the gravity field. Global Geopotential Models (GGMs) are sets of spherical harmonic coefficients representing the earth’s gravity field at different wavelengths. GGMs developed by scientists are published by the International Centre for Gravity Earth Models (ICGEM). When the structure of the GGMs is examined, it is seen that they consist of different degrees and various data groups. The accuracy and resolution provided by each GGM vary depending on the degree of the model and the data used for developing the GGM. Also, geodetic quantities such as potential, geoid undulation, deflection of the vertical components, gravity, and anomaly values can be derived from GGMs within the framework of mathematical principles. In this study, gravity, geoid heights, and free-air gravity anomaly values at test points using different GGMs produced in recent years have been investigated. The study area contains 98000 test points chosen from the Auvergne test area in France. Within the scope of the study, the geoid undulations, gravity values, and free-air gravity anomalies for all points derived from seven recent GGMs have been compared with ground-truth data and the statistical results have been obtained.

On-grid compressive sampling for spherical field measurements in acoustics.

We derive a compressive sampling method for acoustic field reconstruction using field measurements on a pre-defined spherical grid that has theoretically guaranteed relations between signal sparsity, measurement number, and reconstruction accuracy. This method can be used to reconstruct band limited spherical harmonic or Wigner D-function series (spherical harmonic series are a special case) with sparse coefficients. Contrasting typical compressive sampling methods for Wigner D-function series that use arbitrary random measurements, the new method samples randomly on an equiangular grid, a practical and commonly used sampling pattern. Using its periodic extension, we transform the reconstruction of a Wigner D-function series into a multi-dimensional Fourier domain reconstruction problem. We establish that this transformation has a bounded effect on sparsity level and provide numerical studies of this effect. We also compare the reconstruction performance of the new approach to classical Nyquist sampling and existing compressive sampling methods. In our tests, the new compressive sampling approach performs comparably to other guaranteed compressive sampling approaches and needs a fraction of the measurements dictated by the Nyquist sampling theorem. Moreover, using one-third of the measurements or less, the new compressive sampling method can provide over 20 dB better de-noising capability than oversampling with classical Fourier theory.

Continuum solutions to the two–center Coulomb problem in prolate spheroidal coordinates

We present a computational and theoretical framework for solving the Schrödinger equation (SE) for the two–center Coulomb problem in prolate spheroidal coordinates when the energy of the SE is positive. A general and robust computer code has been produced that calculates the separation constants, spheroidal harmonic expansion coefficients, regular quasi–radial two–center Coulomb wave functions, and two–center Coulomb phase shifts. These quantities can be calculated over a range of internuclear separations, angular momentum projections, and continuum electron momenta. A representative set of results are presented and compared with previous calculations, excellent agreement is found in many cases while significant disagreements are found in others. Program summaryProgram title:spheroidal-cwCPC Library link to program files:https://doi.org/10.17632/nf5gw7vjnh.1Code Ocean capsule:https://codeocean.com/capsule/7638744Licensing provisions: MITProgramming language: Fortran 90Nature of problem: A robust and accurate computer code for calculating Coulomb wave functions in prolate spheroidal coordinate is necessary for the calculation of cross sections for photonionization, and electron and positron scattering on charged diatomic targets.Solution method: The angular solution is expanded in a series of spherical harmonics leading to the diagonalization of a pentadiagonal matrix. The eigenvectors contain the expansion coefficients and the eigenvalues are the separation constants which, unlike the SE in spherical coordinates, are also functions of the energy. The radial two–center Coulomb wave functions are started using a power series solution and then propagated using a linear multistep method. An asymptotic expansion is used to calculate the two–center Coulomb phase shift, and to normalize the radial wave function.Additional comments including restrictions and unusual features: Two–center Coulomb radial wave functions and phase shifts can be accurately calculated for internuclear separations from 0.001 to 60.0 a0, angular momentum projection |m| from 0 to 44, momentum k from 0.01 to 100.0 a.u., and a combined nuclear charge Z+ from 1 to 12.

A New Normalization Form Based on a General Unified Normalization Computing for the Gravitational Potential Tensor of Arbitrary Orders

Taking a general unified recursive algorithm for computing the normalized gravitational potential based on Cunningham's formula into account, the numerical stability of both the recursion for spherical harmonic series and the pre-processing for the spherical harmonic coefficients are discussed. Then a new method of normalizing the associated Legendre polynomial is proposed.

Convergence Acceleration of Infinite Series Involving the Product of Riccati–Bessel Function and Its Application for the Electromagnetic Field: Using the Continued Fraction Expansion Method

A summation technique has been developed based on the continuous fractional expansion to accelerate the convergence of infinite series involving the product of Riccati–Bessel functions, which are common to electromagnetic applications. The series is transformed into a new and faster convergent sequence with a continued fraction form, and then the continued fraction approximation is used to accelerate the calculation. The well-known addition theorem formula for spherical wave function is used to verify the correctness of the algorithm. Then, some fundamental aspects of the practical application of continuous fractional expansion for Mie scattering theory and electromagnetic exploration are considered. The results of different models show that this new technique can be applied reliably, especially in the electromagnetic field excited by the vertical electric dipole (VED) source in the “earth-ionospheric” cavity. The comparison among the new technology, the Watson-transform, and the spherical harmonic series summation algorithm shows that this new technology only needs less than 120 series items which is already enough to obtain a small relative error, which greatly improves the convergence speed, and provides a new way to solve the problem.