Exponential shapelets: basis functions for data analysis of isolated features

TL;DR

Exponential shapelets are a new set of basis functions derived from quantum mechanics that efficiently model isolated features in data, outperforming traditional shapelets in data compression and analysis.

Contribution

The paper introduces exponential shapelets, a novel basis function set based on hydrogen atom eigenfunctions, extending shapelet methods to 2D and non-circular features with improved data compression.

Findings

Exponential shapelets outperform Gauss-Hermite/Gauss-Laguerre shapelets in data compression.

They inherit mathematical properties under Fourier transform, facilitating convolution operations.

Applications demonstrated in astronomy, physics, and geodesy.

Abstract

We introduce one- and two-dimensional `exponential shapelets': orthonormal basis functions that efficiently model isolated features in data. They are built from eigenfunctions of the quantum mechanical hydrogen atom, and inherit mathematics with elegant properties under Fourier transform, and hence (de)convolution. For a wide variety of data, exponential shapelets compress information better than Gauss-Hermite/Gauss-Laguerre (`shapelet') decomposition, and generalise previous attempts that were limited to 1D or circularly symmetric basis functions. We discuss example applications in astronomy, fundamental physics and space geodesy.

Figures

Equations

$${\Psi }_n^{+}(x;\beta )=\frac{(-1{)}^{n-1}}{\sqrt{n^3\beta }}\frac{2x}{n\beta }{L}_{n-1}^1\left(\frac{2x}{n\beta }\right)\exp \left(-\frac{x}{n\beta }\right)\forall x⩾0,$$

$$f(x)=\sum _{n=1}^{\infty }{f}_n{\Psi }_n(x;\beta ),$$

$$f_n={\int }_0^{\infty }f(x){\Psi }_n(x;\beta )\mathrm{d}x.$$

$$\sum _{n=1}^{\infty }|f_n{|}^2⩽||f|{|}^2,$$

$${\theta }_{\max }\approx \sqrt{\frac{\int x^2{\Psi }_{n_{\max }}(x;\beta )\mathrm{d}x}{\int {\Psi }_{n_{\max }}(x;\beta )\mathrm{d}x}}=\beta \sqrt{2n_{\max }(2n_{\max }+1)}.$$

$${\theta }_{\min }\approx \beta .$$

$$\begin{array}[]{ccc}\tilde{f}(k)&=&(2\pi)^{-\frac{1}{2}}\int_{-\infty}^{\infty}f(x){\rm e}^{ikx}\,{\rm d}x,\\ f(x)&=&(2\pi)^{-\frac{1}{2}}\int_{-\infty}^{\infty}\tilde{f}(k){\rm e}^{-ikx}\,{\rm d}k\end{array}$$

$${\tilde{\Psi }}_n(k;\beta )=(-1{)}^n\sqrt{\frac{2n\beta }{\pi }}\frac{(n\beta k-i{)}^{2n}}{[(n\beta k{)}^2+1{]}^{n+1}}.$$

$$\left|{\tilde{\Psi }}_n(k;\beta )\right|=\sqrt{\frac{2\pi }{n\beta }}\frac{1}{\pi }\frac{\Gamma /2}{k^2+(\Gamma /2{)}^2}.$$

$$\mathcal{L}({\Psi }_n(x;\beta ))(s)=(-1{)}^{n-1}2\sqrt{n\beta }\frac{(n\beta s-1{)}^{n-1}}{(n\beta s+1{)}^{n+1}}.$$

$$h(x)\equiv (f\star g)(x)\equiv {\int }_{-\infty }^{\infty }\mathrm{d}{x}^{\prime }f(x-x^{\prime })g(x^{\prime }).$$

$$h_n=\sum _{m,l=1}^{\infty }{C}_{nml}{f}_m{g}_l,$$

$$\begin{gathered} C_{nml}(\gamma ,\alpha ,\beta )=\sqrt{\frac{8\alpha \beta \gamma nml}{\pi }}\sum _{u=0}^{2m}\sum _{v=0}^{2l}\sum _{w=0}^{2n}{i}^{u+v+w}\times \\ \left(\genfrac{}{}{0ex}{}{2m}{u}\right)\left(\genfrac{}{}{0ex}{}{2l}{v}\right)\left(\genfrac{}{}{0ex}{}{2n}{w}\right)(m\alpha {)}^u(l\beta {)}^v(n\gamma {)}^n{\mathcal{I}}_{u,v,w}^{m,l,n} \end{gathered}$$

$${\mathcal{I}}_{u,v,w}^{m,l,n}\equiv {\int }_{-\infty }^{\infty }\frac{(-1{)}^w{k}^{u+v+w}\mathrm{d}k}{{\left[(m\alpha k{)}^2+1\right]}^{m+1}{\left[(l\beta k{)}^2+1\right]}^{l+1}{\left[(n\gamma k{)}^2+1\right]}^{n+1}}$$

$${\int }_0^{\infty }{\Psi }_n(x;\beta )\mathrm{d}x=2\sqrt{n\beta }.$$

$$f_{n;\beta }^{\prime }=a^{-1/2}{f}_{n;a\beta },$$

$$F\equiv {\int }_{-\infty }^{\infty }f(x)\mathrm{d}x=2\sqrt{\beta }\sum _{n=1}^{\infty }\sqrt{n}{f}_n,$$

$$x_c\equiv \frac{1}{F}{\int }_{-\infty }^{\infty }xf(x)\mathrm{d}x=\frac{4\sqrt{2}{\beta }^{3/2}}{F}\sum _{n=1}^{\infty }{n}^2{f}_n,$$

$$r_{\mathrm{c}}^2\equiv \frac{1}{F}{\int }_{-\infty }^{\infty }{x}^{2}f(x)\mathrm{d}x=\frac{4{\beta }^{5/2}}{F}\sum _{n=1}^{\infty }{n}^3(2n^2+1)f_n.$$

$$\begin{gathered} {\Psi }_{n,m}(r,\varphi ;\beta )=(-1{)}^n\sqrt{\frac{2}{\beta \pi (2n+1{)}^3}\frac{(n-|m|)!}{(n+|m|)!}} \\ \times {\left(\frac{2r}{\beta (2n+1)}\right)}^{|m|}{L}_{n-|m|}^{2|m|}\left(\frac{2r}{\beta (2n+1)}\right) \\ \times \exp \left(-\frac{r}{\beta (2n+1)}\right)\exp (-im\varphi ), \end{gathered}$$

$$f(r,\varphi )=\sum _{n=0}^{\infty }\sum _{m=-n}^n{f}_{n,m}{\Psi }_{n,m}(r,\varphi ;\beta ),$$

$$f_{n,m}={\int }_0^{2\pi }{\int }_0^{\infty }f(r,\varphi ){\Psi }_{n,m}(r,\varphi ;\beta )r\mathrm{d}r\mathrm{d}\varphi .$$

$${\theta }_{\max }$$

$${\theta }_{\min }\approx \beta .$$

$${\tilde{\Psi }}_{n,m}(\vec{k};\beta )=2\pi i^m{\mathrm{e}}^{-im{\varphi }_k}{\left[F_m(k)\right]}_{n,\beta },$$

$$\bar{f}(r)\equiv \frac{1}{2\pi }{\int }_0^{2\pi }f(r,\varphi )\mathrm{d}\varphi .$$

$$\bar{f}(r)=\sum _{n=0}^{\infty }{f}_{n,0}{\Psi }_{n,0}(r;\beta ),$$

$$\begin{gathered} {\Psi }_{n,0}(r;\beta )=(-1{)}^n\frac{2}{\beta \sqrt{2\pi }}(2n+1{)}^{-3/2} \\ \times \exp \left(-\frac{r}{\beta (2n+1)}\right)L_n\left(\frac{r}{\beta (n+\frac{1}{2})}\right). \end{gathered}$$

$$F_R\equiv {\int }_0^{2\pi }{\int }_0^{R}f(r,\varphi )r\mathrm{d}r\mathrm{d}\varphi .$$

$$L_n^{\alpha }(x)=\sum _{k=0}^n(-1{)}^k(\genfrac{}{}{0ex}{}{n+\alpha }{n-k})\frac{x^k}{k!}.$$

Full text

1 Introduction

A frequent task in data analysis is to categorise and quantify the shapes of localised objects – such as transient events in a (one-dimensional) time-series, or regions of interest in a (two-dimensional) image. It is such a universal challenge that methods developed for one field frequently turn out to be useful in others. For example, astrophysicists measure the shapes of distant galaxies by decomposing them into orthogonal basis functions, such as CHEFs (Jiménez-Teja & Benítez, 2012) or (Gaussian) shapelets (Bernstein & Jarvis 2002; Refregier 2003; Refregier & Bacon 2003; Massey & Refregier 2005). Shapelets have been used to analyse data in other branches of astrophysics, modelling extrasolar planets (Hoekstra et al., 2005; Amara & Quanz, 2012), the distribution of dark matter (Birrer et al., 2015; Tagore & Jackson, 2016), or flashes of pulsars (Lentati et al., 2015; Desvignes et al., 2016; Ellis & Cornish, 2015). They have also been used in medical imaging (Weissman, Hancewicz & Kaplan, 2004; Apostolopoulos et al, 2017), pattern recognition in human vision (Sharpee & Victor, 2009) or artificial vision (Sabzmeydani & Mori, 2007), data compression (Holbrey, 2006), and the manufacture of nanoscale thin films (Suderman et al., 2015; Akdeniz et al., 2018).

Gaussian shapelets are based on eigenfunctions of the quantum mechanics harmonic oscillator. In one-dimensional form, they are Gauss-Hermite functions, which seem to be well adapted in several time-series analyses, especially when transients (whose shape is close to damped oscillations) must be detected, characterised and/or corrected for. This is the case for instance in fundamental physics for MICROSCOPE ‘crackles’ (e.g. Baghi et al. 2015; Bergé et al. 2015), space geodesy for GRACE ‘twangs’ (Flury et al., 2008; Peterseim et al., 2010, 2014) or ‘glitches’ in gravitational waves searches with LIGO (Cornish & Littenberg 2015; Powell et al. 2015, 2017; Principe & Pinto 2017 and references therein). In two-dimensional form, they can be expressed equivalently as either Gauss-Hermite (Cartesian; Refregier 2003) or Gauss-Laguerre (polar; Bernstein & Jarvis 2002; Massey & Refregier 2005) functions. Owing to their quantum mechanical origin, they have elegant properties (they make a complete set) under Fourier transform, and are hence efficient for operations involving convolution or deconvolution of two images.

The main limitation of shapelets is that they are perturbations around a Gaussian, which is flat near its peak. They inefficiently parameterise peaky features, including the distant galaxies for which they were originally suggested (Melchior et al., 2010). Galaxies have a two-dimensional surface density that decreases approximately exponentially with distance from the centre. Capturing the steep gradient near the centre requires a weighted sum of many Gaussian shapelet basis functions, which then overfit noise in the extended wings. Attempting to overcome this limitation, Ngan et al. (2009) developed ‘Sersiclets’, although they were forced to be circularly symmetric, and have not seen wider applications.

In this paper, we extend the quantum mechanical framework underpinning Gaussian shapelets, and define 1D and 2D exponential shapelets based on wavefunctions of the hydrogen atom. These functions are perturbations around a decreasing exponential, and should efficiently model any arbitrarily-shaped but centrally-peaked regions of interest. The perturbations themselves are Laguerre polynomials.

We introduce 1D exponential shapelets in Section 2, and 2D exponential shapelets in Section 3. In both cases, we describe their main properties, compare them to Gaussian shapelets, and provide example uses. We conclude in Section 4.

2 1D exponential shapelets

2.1 Definition

The one-dimensional hydrogen atom is the solution to the motion of a particle in a 1D Coulomb potential $1/|x|$. In this paper, we will neither dwell on its rich history, nor the debates about the finitude of its ground state, and about the existence of even wavefunctions (see Nieto 1979; Palma & Raff 2006; Núñez-Yépez et al. 2011, 2014 and references therein). Instead, we simply exploit the normalized 1D hydrogen atom wavefunctions as given by Palma & Raff (2006) to define the 1D exponential shapelet basis functions

[TABLE]

for $n\geqslant 1$, where $L_{n-1}^{1}$ are the generalized Laguerre polynomials (see e.g. Massey & Refregier 2005). The characteristic scale size $\beta$ corresponds to the Bohr radius, and $n$ is the eigenfunction energy level.

Note that, contrary to normal procedure in quantum physics, we restrict 1D exponential shapelets on positive $x$, and refrain from defining the negative-$x$ part $\Psi_{n}^{-}(x;\beta)=-\Psi_{n}^{+}(-x;\beta)$ (for $x<0$). Hence, we do not follow Palma & Raff (2006) and do not define the even and odd wavefunctions $\Psi_{n}^{+}(x;\beta)\pm\Psi_{n}^{-}(x;\beta)$. Events in time series can often be adequately described by their behaviour after a certain moment (in this case, $x=0$). Henceforth, we shall therefore drop the + superscript, to write more simply $\Psi_{n}(x;\beta)\equiv\Psi_{n}^{+}(x;\beta)$. These functions are continuously differentiable, smoothly departing from zero at $x=0$ and tending back to zero as $x\rightarrow\infty$.

Because the basis functions (1) are eigenfunctions of the one-dimensional hydrogen atom’s Hamiltonian, they form an orthogonal basis of the square integrable functions $L^{2}([0,\infty[,\langle\cdot,\cdot\rangle)$ Hilbert space equiped with the inner product $\langle f(x),g(x)\rangle=\int_{0}^{\infty}f(x)g^{*}(x){\rm d}x$ and an asterisk denotes complex conjugation. They are orthonormal, in the sense that $\int_{0}^{\infty}\Psi_{n}(x;\beta)\Psi_{m}(x;\beta)\,{\rm d}x=\delta_{nm}$ where $\delta_{nm}$ is the Kronecker symbol, and complete, in the sense that $\sum_{n=1}^{\infty}\Psi_{n}(x;\beta)\Psi_{n}(x^{\prime};\beta)=\delta(x-x^{\prime})$.

They thus form a basis, on which we can uniquely decompose a function as

[TABLE]

with coefficients $f_{n}$ given by an overlap integral

[TABLE]

Bessel’s inequality then assures us that for any function $f\in L^{2}([0,\infty[,\langle\cdot,\cdot\rangle)$,

[TABLE]

where $||.||$ is the L2 norm, so that the series (2) converges and coefficients $f_{n}$ must vanish as $n$ increases. The series can therefore be truncated for suitably localised functions, at some value $n\leqslant n_{\rm max}$.

2.2 Properties

2.2.1 Maximum and minimum effective scales

The first few 1D exponential shapelet basis functions $\Psi_{n}(x;\beta)$ are shown in Fig. 1. As the order $n$ increases (with constant $\beta$), the largest-scale oscillation dominates, and rapidly acquires a larger extent. However, the smallest oscillations always remain roughly the same size.

To act as a convenient figure-of-merit, we follow Refregier (2003) in defining the largest scale that can be described by an exponential shapelet model as

[TABLE]

Empirically, the smallest scale that can be described is

[TABLE]

As ${n_{\rm max}}\rightarrow\infty$, the range of scales modelled by a 1D exponential shapelet decomposition, $\theta_{\rm max}/\theta_{\rm min}\propto n_{\rm max}=n_{\rm coeffs}$. However, the resolution of an exponential shapelet model is greatest near the origin, and decreases away from it. This behaviour is very different from Gaussian shapelets, where the resolution is more spatially uniform: increasing $n$ increases the scale of the model slowly, while simultaneously adding smaller-scale oscillations (see Fig. 1 of Refregier 2003). The wide effective range of scales for exponential shapelets is what will allow them to efficiently describe spiky but long-duration events.

2.2.2 Fourier and Laplace transforms

To compute the Fourier transform of the 1D exponential shapelet basis functions, we adopt a convention where the Fourier transform is defined as

[TABLE]

for any function $f(x)$, where $i^{2}=-1$. The Fourier transform of the basis function $\Psi_{n}(x;\beta)$ is then (setting $f(x)=0$ for $x<0$)

[TABLE]

Note that the real part of the Fourier transform is even for all $n$, while its imaginary part is odd for all $n$. Its modulus is a Lorentzian function centered on $k=0$ whose width parameter $\Gamma\equiv 2/(n\beta)$ is inversely proportional to $\beta$ (as typical for exponentially suppressed oscillations)

[TABLE]

With conventions similar to those for the Fourier transform, the Laplace transform is given by

[TABLE]

Although Eqs. (8) and (10) are not as simple as for Gaussian shapelets (the Fourier transform of a Gaussian shapelet is itself a Gaussian shapelet), closed forms exist for exponential shapelets that are easy to implement. The Fourier and Laplace transforms of 1D exponential shapelets are also well localized in frequency space (Fig. 2).

2.2.3 Convolution

Convolution is an inevitable operation during signal acquisition, whereby an input (‘real’) signal undergoes the measurement apparatus’ transfer function. What is ultimately measured is the convolution of the input signal with the transfer function.

Let us consider the convolution of two functions $f(x)$ and $g(x)$, whose convolution product is

[TABLE]

Following the same arguments as Refregier (2003), the functions can all be decomposed into exponential shapelets, perhaps with different scale sizes $\alpha$, $\beta$ and $\gamma$. We can then relate the 1D exponential shapelet space coefficients $h_{n;\gamma}$ to $f_{m;\beta}$ and $g_{l;\alpha}$ via

[TABLE]

where the convolution tensor is given by

[TABLE]

and $\binom{n}{m}$ is the binomial coefficient. A proof of Eq. (13) is provided in Appendix A.1.

The integral

[TABLE]

appearing in Eq. (13) is zero if $u+v+w$ is odd (a proof is provided in Appendix A.2). Hence, $C_{nml}(\gamma,\alpha,\beta)$ is always wholly real. Under some conditions, it can be estimated analytically and expressed as a sum of converging series (see Eqs. 59, 63 and 70 in Appendix A.3). Those conditions are often not met, but the integrand is a well-behaved function, so the integration can also be performed numerically. We find that numerical integration is most convenient if the function is cut into three segments (see Eqs. 55, 61 and 68 in Appendix A.3).

2.2.4 Rescaling

1D exponential shapelets obey the integral property

[TABLE]

Using this, it can be shown that under a rescaling $x\rightarrow x^{\prime}=ax$, the coefficients $f_{n}$ of a model (3) transform as

[TABLE]

and under $f(x)\rightarrow f^{\prime}(x)=kf(x)$, the coefficients are themselves multiplied by $k$.

2.3 Shape characterization

Let $f(x)=\sum_{n}f_{n}\Psi_{n}(x;\beta)$ be a 1D object decomposed into 1D exponential shapelets. In terms of its shapelet coefficients, its integral (total ‘flux’) is

[TABLE]

which can be readily shown from the integral property (15).

Provided $F\neq 0$, its barycenter (centroid) is

[TABLE]

and it has a characteristic size

[TABLE]

This size can be used to determine the exponential decay rate of the object, for instance, when modelling damped oscillations.

Note that series (19) converges only if the amplitudes of the 1D exponential shapelet coefficients decrease faster than $n^{-6}$, which may not always be the case. Care must therefore be taken to check for convergence when characterising the shape of a feature using this technique.

2.4 Exponential shapelets modelling in practice

As shown above (Eq. 2), a 1D feature is straightforward to model for a given couple ($n_{\rm max}$, $\beta$), where $n_{\rm max}$ is the maximum order of the truncated sum (2). For example, a linear least-square method can be efficiently used for this purpose. Then the model depends non-linearly on the two parameters $n_{\rm max}$ and $\beta$ that can be optimised by iteratively minimizing the residuals between the observed feature and its model. This procedure was described at length, in the 2D case, by Massey & Refregier (2005).

2.5 Example applications

This section demonstrates three possible applications of 1D exponential shapelets. We start by modelling exponentially suppressed oscillations, which are measurements throughout experimental physics, including the response of accelerometers111http://www.onera.fr/en/dphy onboard the space missions MICROSCOPE (Touboul et al., 2017), GRACE (Tapley et al., 2004) and GOCE (Rummel et al., 2011). We then discuss a potential application to unmodelled bursts in the analysis of gravitational waves. Exponential shapelets may also be convenient to model charge transfer inefficiency trailing due to radiation damage in spacebourne imaging detectors (Massey et al., 2010), although we do not explore that further here.

2.5.1 Cleaning accelerometer time series data,                     and modelling an experiment’s response function

MICROSCOPE tested the Weak Equivalence Principle (WEP) by precisely measuring the differential acceleration experienced by two concentric cylindrical test-masses onboard a drag-free satellite in Low Earth Orbit222The mission ended on October 16, 2018; data analysis is still underway.. In theory, any non-zero difference at a well defined frequency $f_{\rm EP}$ (which depends on the orbit and attitude characteristics of the satellite) is a signature of the violation of the WEP.

In practice, many transients are apparent in MICROSCOPE data (the upper panel of Fig. 3 shows a typical high-signal-to-noise example). These transients are generally caused by crackles of the satellite’s coating (because of temperature variations), crackles of the satellite’s gas tanks (as their pressure decreases as the gas is consumed), or impacts with micro-meteorites. Such transients occupy frequencies higher than $f_{\rm EP}$, so they do not directly impact a possible WEP violation signal. However, it is necessary to detect and mask them in measured time-series, then either reconstruct the corresponding ‘missing data’ (Bergé et al., 2015; Pires et al., 2016) to allow for a least square fit of the expected WEP violation signal in the frequency domain (Touboul et al., 2017), or adapt the maximum likelihood technique to take missing data into account (Baghi et al., 2015, 2016). Existing techniques are suboptimal, as they may affect the noise characteristics. Moreover, transients could be considered as conveying useful information. They are created by an external impulse; if this is assumed to be instantaneous (Dirac), the shape and relaxation time of the observed signal is a measurement of MICROSCOPE’s transfer function.

1D exponential shapelets are well-matched to these transient signals. The upper panel of Fig. 3 shows a transient in the time domain: the observed data are shown in black, a model (fitted between $t=10$s and $t=25$s) and its 68% confidence interval are shown in red, and residuals consistent with noise are shown in green. This model uses $n_{\rm max}=15$, meaning that it has compressed 60 data points in the time domain description (15 seconds sampled at 4Hz) into 15 shapelet coefficients. Nonetheless, the model has rich, empirical freedom to capture multiple response modes of the instrument’s complex structure. The lower panel of Fig. 3 shows the convergence of the model as we increase $n_{\rm max}$, quantified as the square residuals between the observed transient and the model. Most of the information is contained in coefficients $3\leqslant n\leqslant 13$.

An extension to this process will be presented in a future paper. By fitting 1D exponential shapelet coefficients to many transients, it is possible to model temporal variations in MICROSCOPE’s instrument’s relaxation time via parameter $\beta$ or the exponential envelope scale $r_{\rm rms}$. Interpolating these models then yields a model of the transfer function at any time, to either gain insight into the instrument’s performance, or to correct (deconvolve) signals with an ‘inverse transfer function’ transform. Note that this procedure is similar to techniques developed for 2D astronomical image processing, where a Point Spread Function (PSF) is determined from the shape of stars then interpolated across the data (see e.g. Bergé et al. 2012; Gentile et al. 2013; Kilbinger 2015).

2.5.2 Characterisation of perturbing signals in space-borne geodesy missions

The GRACE mission (Tapley et al., 2004) revolutionised geodesy by measuring the Earth’s gravitational field with unprecedented precision. Two satellites followed each other on the same orbit, monitoring the distance between them via microwave ranging. In theory, any variation in their relative speed or distance can be ascribed to variations in the Earth’s geopotential.

In practice, the satellites were also subject to non-gravitational forces. These were measured by ultrasensitive accelerometers, for removal during post-processing (Flury et al., 2008; Peterseim et al., 2010). Peterseim (2014) and Peterseim et al. (2014) modelled transients (known as ‘spikes’) in accelerometer data using a piecewise function made from the derivative of a Gaussian, a 3rd order polynomial and a damped oscillation. Some were successfully classified as due to either the satellite’s heaters, or activation of its magnetic torquers – but no physical origin could be assigned to others, known as ‘twangs’.

Tentative correlations of twangs with the position of the satellite along its orbit hint at a possible geophysical origin. Both categories of twangs are compactly-represented as 1D exponential shapelets, so we will investigate in a future paper whether these provide a cleaner set of shape parameters to categorise and understand their origin.

2.5.3 Unmodelled bursts and glitches in gravitational waves data analysis

A wealth of methods have been developed to search for, characterise and classify unmodelled bursts and instrumental glitches in searches for gravitational waves (see e.g. Cornish & Littenberg, 2015; Powell et al., 2015, 2017; Principe & Pinto, 2017, and references therein). Glitches are often modelled using ‘Sine-Gaussian waveforms’ (Principe & Pinto, 2017). These have similar properties to Gaussian shapelets, although shapelets can encode more details. 1D exponential shapelets could further optimise the data compression of complex glitch shapes that often exhibit damped oscillations.

Shapelets might therefore improve the detection, characterisation and classification of instrumental glitches in gravitational wave detectors. Indeed, bursts and glitches are usually detected in the time-frequency domain, which is easily reproduced in Gaussian shapelets-time space. If 1D exponential shapelets more efficiently model the information in a glitch, exponential shapelet-time space would be even better. We will investigate in a future paper whether glitches can be detected by scanning a matched-filter, and correlating the measured signal with a 1D exponential shapelet model, while leaving $\beta$ (and possibly $n_{\rm max}$) free.

2.6 Comparison with Gaussian shapelets

The MICROSCOPE transient modelled as exponential shapelets ($n_{\rm max}=15$) in Fig. 3, is modelled as Gaussian shapelets ($n_{\rm max}=18$) in Fig. 4. Achieving a similar precision of fit requires more coefficients, and creates more small-scale artefacts: the data are over-fit near the centre of the model, and underfit at the extremes. We also find that its coefficients are highly covariant, with the good fit produced by large positive and negative basis functions almost precisely cancelling each other. This is much less apparent for the 1D exponential shapelet coefficients. Indeed, we find that estimating the model again on points different from those observed data fails to reproduce the overall shape of the transient, whereas the 1D exponential shapelet model is itself predictive.

In this example, it appears that 1D exponential shapelets outperform Gaussian shapelets. This is because their perturbations around an exponential decay are better matched to the underlying signal, and because their wider extent spreads information density more evenly. In general, the type of shapelet to choose should depend on the decay rate of the target function.

3 2D exponential shapelets

3.1 Definition

The quantum mechanics of a hydrogen atom restricted to two dimensions is well studied (see e.g. Zaslow & Zandler, 1967; Yang et al., 1991; Chaos-Cador & Ley-Koo, 2007). The bound states of the electron (i.e. eigenfunctions of the Hamiltonian) provide a natural set of basis functions to represent a bounded distribution. After renormalizing these, we define the 2D exponential shapelet basis functions

[TABLE]

for $n\geqslant 0$, where $L_{i}^{j}\left(x\right)$ are generalized Laguerre polynomials. The $(-1)^{n}$ term is used to ensure that the integral of each basis function is positive, but this is otherwise the form used in quantum theory. In terms of quantum mechanics, $n$ is the principal quantum number (the eigenfunction energy level) and $m$ the magnetic quantum number, which takes integer values between $-n$ and $n$. The characteristic scale $\beta$ is linked to the Bohr radius.

These functions form an orthonormal set of basis functions in the $L^{2}([0,\infty[\times[0,2\pi[,\langle\cdot,\cdot\rangle)$ space equipped with inner product $\langle\Psi_{n,m}(r,\phi),\Psi_{n^{\prime},m^{\prime}}(r,\phi)\rangle=\int_{0}^{2\pi}\int_{0}^{\infty}\Psi_{n,m}(r,\phi)\,\Psi^{*}_{n^{\prime},m^{\prime}}(r,\phi)\,r\,{\rm d}r\,{\rm d}\phi=\delta_{nn^{\prime}}\delta_{mm^{\prime}}$.

Any localized function $f(r,\phi)$ can be uniquely decomposed into a weighted sum of these basis functions

[TABLE]

where the 2D exponential shapelet coefficients are given by

[TABLE]

Using Bessel’s inequality like in the 1D case, we find that series (21) converges, and the coefficients $f_{n,m}$ must vanish when $n$ and $m$ increase. For a real function $f(r,\phi)\in\mathbb{R}$, $f_{n,m}=f_{n-m}^{*}$. In this case, truncating the series at $n\leqslant n_{\rm max}$ results in $n_{\rm coeffs}=(n_{\rm max}+1)^{2}$ independent coefficients.

It may also be possible to define elliptical 2D exponential shapelets by applying a shear transformation to the coordinate system, as Nakajima & Bernstein (2007) suggested for Gaussian shapelets. This preserves their orthonormality, and can increase their rate of data compression, at the cost of two additional non-linear parameters to specify the shear.

3.2 Properties

3.2.1 Maximum and minimum effective scales

The first few 2D exponential shapelet basis functions are illustrated in Fig. 5, and their radial component is shown in Fig. 6, defined via $\Psi_{n,m}(r,\phi;\beta)\equiv R_{n,m}(r;\beta)\exp(-im\phi)$. Their resemblance with Gaussian polar shapelets is striking (see Fig. 2 of Massey & Refregier, 2005). However, due to their exponential kernel, exponential shapelets are both more peaky and further spread out than Gaussian shapelets. The size difference between the lowest order basis function $\Psi_{00}$ and even a low order one like $\Psi_{40}$ is striking. It is this behaviour that will help to describe spatially extended features.

Generalising the 1D case (Section 2.2.1), we define the largest effective scale in a 2D exponential shapelet model as

[TABLE]

Again empirically, the smallest resolved scales are roughly constant,

[TABLE]

The range of scales included in a 2D exponential shapelet model as $n_{\rm max}\rightarrow\infty$ is $(\theta_{\rm max}/\theta_{\rm min})^{2}\propto n_{\rm max}^{4}\propto n_{\rm coeffs}^{2}$. The resolution of an exponential shapelet model is greatest near the origin, and information density is concentrated there. This behaviour is different from the Gaussian shapelets, where resolution is more spatially uniform, and information density is constant, with $(\theta_{\rm max}/\theta_{\rm min})^{2}=n_{\rm max}^{2}+1\propto n_{\rm coeffs}$.

3.2.2 Fourier transformation and convolution

Using the same convention as in the 1D case, it can be shown that the Fourier transform of 2D exponential shapelets is given by

[TABLE]

where $k=|\vec{k}|$, $\phi_{k}$ is the angle between the $\vec{k}$ direction and the $\phi=0$ direction in polar space, and $\left[F_{m}(k)\right]_{n,\beta}$ is the Hankel transform of the $R_{nm}(r,\beta)$ radial function. Consequently, the convolution tensor involved in the convolution of two objects modelled in 2D exponential shapelets is the integral of products of Hankel transforms of $R_{n,m}$ functions. We were not able to find an analytic form for this Hankel transform, so it must instead be computed numerically.

3.2.3 Coordinate transforms

It can be useful to know the response of the basis functions to 2D linear coordinate transforms: either to know how to mix the coefficients of a shapelet decomposition in order to perform that transform, or to form combinations of coefficients that are invariant under some transforms. This was used to construct estimators of the shear distortion applied to images of galaxies by the effect of weak gravitational lensing (Kuijken, 2006; Massey et al., 2007).

A convenient shortcut to calculating those transforms for Gaussian shapelets was provided by the ladder operators associated with the quantum mechanical harmonic oscillator (Refregier & Bacon, 2003). Unfortunately, we have not yet found a useful form of the ladder operators for the 2D hydrogen atom. If it becomes necessary to perform linear transformations on 2D exponential shapelets, it will be necessary to apply the transforms manually to the basis functions (a long and arduous task, but one that is guaranteed to yield a closed form, because the basis is complete).

3.3 Object shape measurement

Once a feature has been decomposed into 2D exponential shapelets, its coefficients can be used to construct characteristic measurements of its shape. In this section, we derive expressions for an object’s (azimuthally-averaged) radial profile, flux, centroid, unweighted quadrupoles, size and ellipticity, in terms of its shapelet coefficients. We have not attempted it here, but it should also be possible to develop an exponential-shapelet classification of galaxy morphologies via e.g. their concentration, asymmetry and symmetry (Conselice et al., 2000).

3.3.1 Radial profile

Azimuthally averaging an object’s signal yields its mean radial profile

[TABLE]

Noting that all $m\neq 0$ basis functions average to zero, the radial profile reduces to

[TABLE]

where the (rotationally invariant) $m=0$ basis functions are

[TABLE]

Eq. (27) is identical to the equivalent derivation for Gaussian shapelets (see Eq. (16) of Massey & Refregier, 2005), up to the fact that $m=0$ basis functions with odd $n$ do not exist in the Gaussian case.

3.3.2 Flux

Integrating the signal inside a circular aperture of radius $R$ yields its ‘flux’

[TABLE]

To evaluate this integral, it is useful to note that, once again, all $m\neq 0$ basis functions cancel out to zero during integration over $\phi$, and also a closed form for the generalized Laguerre polynomials,

[TABLE]

Using this, it can be shown that

[TABLE]

where $\gamma(y,x)$ is the lower incomplete gamma function. Extrapolating $F_{R}$ to large radius, and taking the limit $R\rightarrow\infty$, we obtain

[TABLE]

3.3.3 Centroid

Similarly, it can be shown that the unweighted centroid ($x_{c},y_{c}$), defined by

[TABLE]

is, in terms of 2D exponential shapelet coefficients,

[TABLE]

3.3.4 Quadrupole moments

The unweighted quadrupole moments of an object

[TABLE]

can be used to define quantities such as its rms size and ellipticity (see below). They are, in terms of 2D exponential shapelet coefficients,

[TABLE]

[TABLE]

[TABLE]

3.3.5 Size and ellipticity

Using Eqs. (36)-(38), expressions can be obtained to quantify a feature’s size

[TABLE]

and ellipticity

[TABLE]

In this complex notation with $\varepsilon\equiv|\varepsilon|\cos{(2\phi)}+i|\varepsilon|\sin{(2\phi)}$, $\varepsilon=0$ denotes rotational invariance, and a positive real (imaginary) component denotes elongation along (at $45^{\circ}$ to) the $x$-axis.

3.3.6 Convergence

Expressions for the shape estimators resemble those obtained for polar Gaussian shapelets in Section 6 of Massey & Refregier (2005). In particular, a feature’s radial profile, flux and size are encoded within the $n=0$ coefficients, while the centroid is described by the $n=1$ coefficients and ellipticity by the $n=2$ coefficients. The expressions also contain the same power of $\beta$ as those for Gaussian shapelets, because that encodes information about the units in which the data is expressed.

However, shape measures for 2D exponential shapelets contain higher powers of $n$ than those for 2D Gaussian shapelets, and will converge more slowly. Some of this is due to the normalisation of the basis functions against an inner product. In the Gaussian shapelet case, this happens to yield an integral of basis functions (32) that is independent of $n$, while in exponential shapelets it is $\propto n^{3/2}$. That power of $n$ could have been included in the basis functions and removed from the coefficients – although doing so would merely make them look more stable rather than actually altering convergence. In the limit that $R\rightarrow\infty$, convergence of the flux estimator requires $|f_{n,0}|$ to decrease faster than $n^{-5/2}$; convergence of centroid requires $|f_{n,1}|$ to decrease faster than $n^{-11/2}$; convergence of size and ellipticity requires $|f_{n,0}|$ and $|f_{n,2}|$ to decrease faster than $n^{-13/2}$. This may not be a problem, since we expect 2D exponential shapelets to converge faster than shapelets (i.e. to have a lower $n_{\rm max}$ for a given galaxy). However, if this does not hold true, for example due to measurement noise, care should be taken to restrict the decomposition and measurement of photometry inside a finite aperture, or to truncate the 2D exponential shapelet model at finite $n_{\rm max}$.

3.4 Example application

This section demonstrates the use of 2D exponential shapelets to model the shapes of galaxies seen outside our own Milky Way (they may also be convenient to model the shapes of gamma-ray events; S. Pires, private communication). We examine the convergence speed of 2D exponential shapelets, and compare their performance to that of 2D Gaussian shapelets.

3.4.1 Distant galaxies in astronomical imaging

Galaxies are collections of (hundreds of billions of) stars, with a combined surface brightness that peaks sharply near the centre, and decreases to large radii in a way that is often exponential. Four galaxies in the COSMOS survey (Scoville et al., 2007a, b) observed by the Hubble Space Telescope (Koekemoer et al., 2007), are shown in Fig. 7. As mentioned in Sect. 2.4 for the 1D case, we follow the algorithm described by Massey & Refregier (2005) to choose the best non-linear parameters of the models –$n_{\rm max}$, $\beta$ and centroid). Note that galaxies in this figure were modelled into shapelets with the “diamond” truncation scheme introduced in Massey et al. (2007), which is intended to reduce small scales oscillations by ignoring the higher $m$ terms of the decomposition. In this case, the number of coefficients is divided by 2, allowing also a better compression of the information.

In all four cases, 2D exponential shapelets require fewer coefficients (lower $n_{\rm max}$; here $n_{\rm coeffs}=85$) to provide a model whose residuals are consistent with the noise than Gaussian shapelets (here $n_{\rm coeffs}=121$). Even more interestingly, exponential shapelets tend to provide a (visually) cleaner model of the central region of the galaxies. The top left galaxy has an irregular morphology, and the 2D exponential shapelet model shows fewer artifacts than the polar shapelet model. For an elliptical (top right) or edge-on spiral galaxy (bottom left), the ringing common in Gaussian shapelets disappears with 2D exponential shapelets: they do not possess high-frequency ripples at large radii that need large coefficients alternating in sign for them to cancel. Even for highly eccentric galaxies (bottom right), 2D exponential shapelets provide cleaner models than polar shapelets, especially in the outskirts of the galaxies, as the exponential wings are fundamentally a better match to the underlying signal.

3.4.2 Convergence speed

To assess the ability of 2D exponential shapelets to model galaxies, we require noise-free images whose true properties are known. We simulate elliptical galaxies with a Sersic radial profile

[TABLE]

where $I(r)$ is the galaxy’s s profile, $n_{s}$ is the Sersic index, $b_{n}$ is a normalisation constant and $r_{e}$ is the effective radius (which contains 50% of the galaxy’s flux). We simulate three types of galaxies: exponential ($n_{s}=1$), intermediate ($n_{s}=2$), and de Vaucouleurs ($n_{s}=4$).

Fig. 8 compares the radial profiles fitted with 2D exponential shapelets (upper row) or Gaussian shapelets (lower row), for different maximum orders of decomposition $n_{\rm max}$, and coarsely optimised $\beta$ (this could be further optimised at low $n_{\rm max}$ in the presence of noise). The normalisation of the ordinate depends on $b_{n}$ and pixellisation, so is arbitrary; only its relative scaling is informative. For all three galaxy types, 2D exponential shapelets greatly outperform Gaussian shapelets. A decomposition to $n_{\rm max}=2$ or $4$ is sufficient to model an exponential galaxy. Depending on the signal-to-noise of real data, a decomposition to $n_{\rm max}=6$ may be sufficient for an intermediate $n_{s}=2$ galaxy. It is more challenging to model a de Vaucouleurs galaxy, and even exponential shapelets require $n_{\rm max}=12$ to capture its extended tails three orders of magnitude fainter than the peak. However, they do so without the oscillatory ‘ringing’ introduced by the cancellation of positive and negative basis functions in Gaussian shapelets.

To illustrate the convergence speed, Fig. 9 shows the mean square residual of the models of the same galaxies as in Fig. 8, as a function of the maximum shapelet order $n_{\rm max}$ (upper row), or the total number of coefficients $n_{\rm coeffs}$. The normalisation of the ordinate is arbitrary, and only its relative scaling is informative. Once again, it is clear that 2D exponential shapelets (solid lines) outperform Gaussian shapelets (dashed lines). For exponential and intermediate galaxies (left and middle columns), the residual of the 2D exponential shapelet model decreases consistently, and ends up several orders of magnitude better than that of Gaussian shapelets. The enhanced performance of exponential shapelets is even more prominent for a de Vaucouleurs galaxy. The limitation for exponential shapelets in this case is numerical precision of our algorithm to pixellate the basis functions. At high-$n$, this breaks down, causing spurious residuals at the very centre of the model, shown as an artificial upturn in the bottom-right panel of Fig. 9. This could be circumvented by more accurate pixellisation, or by imposing limits on $\theta_{\rm min}$.

4 Conclusion

In this paper, we introduced exponential shapelets, a family of orthonormal basis functions that can be efficiently used to model 1D and 2D objects. They borrow many concepts from the original Gaussian shapelets, but are perturbations around a decaying exponential rather than a Gaussian, and more efficiently compress the information in damped 1D oscillations or centrally-concentrated 2D features. In particular, they can simultaneously capture the central peak and the extended wings of galaxies in astronomical imaging, thereby solving the main criticism levelled at Gaussian shapelets in the field for which they were originally intended.

Modelling a feature using exponential shapelets first requires a choice of characteristic scale size $\beta$ and expansion order $n_{\rm max}$; these can be selected using an non-linear optimisation method developed for Gaussian shapelets (Massey & Refregier, 2005). The function is then decomposed into a weighted sum of exponential shapelet basis functions using linear regression. A least-squares fit is often sufficient, although in some cases, modeling a transient may be stabilised by regularisation techniques such as sparsity constraints.

Once a feature is described as a weighted sum of exponential shapelets, simple combinations of the weights yield expressions for its total flux, centroid, size and ellipticity. These are similar to those for Gaussian shapelets. However, convolution and deconvolution in exponential shapelets is significantly more difficult than the Gaussian case, with the convolution tensor being time-consuming to calculate via numerical integration.

We described possible applications of exponential shapelets in several fields of experimental and observational science. 1D exponential shapelets can be used to model (and subtract or understand) spurious transients in time series, such as measurements by the MICROSCOPE or GRACE satellites, or the LIGO search for gravitational waves. Thanks to their exponential decay, 2D exponential shapelets overcome the main criticism aimed at Gaussian shapelets (though at the cost of losing some simplicity) and can be used to measure the brightness, shape and size of distant galaxies in astronomical imaging. The characteristics of exponential shapelets should make them well-suited to data compression and analysis in a wide range of fields; their convenient mathematical properties should see them frequently adopted.

Acknowledgments

The concepts in this paper matured over many years, after a preliminary study to generalise Gaussian polar shapelets for Hubble Space Telescope data, for which we acknowledge support from NASA through grant HST-AR-11747 – and which was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. We are indebted to Barnaby Rowe for researching Kummer functions and the literature of the 2D hydrogen atom. We also acknowledge helpful and motivating discussions with Jason Rhodes, Richard Ellis, Alexandre Refregier, Sandrine Pires, James Nightingale, Bernard Foulon, Bruno Christophe, Gilles Métris, Manuel Rodrigues, Alain Robert and Nicolas Touquoy. This work makes use of technical data from the CNES-ESA-ONERA-CNRS-OCA Microscope mission, and has received financial support from ONERA and CNES. JB acknowledges the financial support of the UnivEarthS Labex program at Sorbonne Paris Cité (ANR-10-LABX-0023 and ANR-11-IDEX-0005-02) and of CNES through the APR programme (LISA project). RM is supported by a Royal Society University Research Fellowship and through STFC grant ST/N001494/1.

Appendix A Derivation of convolution kernel tensor $C_{nml}$ for 1D exponential shapelets

A.1 Proof of Eq. (13)

Let $f(x)$ and $g(x)$ two 1-dimensional functions, whose convolution product is $h(x)=f(x)\star g(x)$. In Fourier space, $\tilde{h}(k)=\tilde{f}(k)\tilde{g}(k)$, where

[TABLE]

and similarly for $\tilde{g}(k)$.

Decomposing $f(x)$ and $g(x)$ into exponential shapelets ($f(x)=\sum_{m=1}^{\infty}f_{m}\Psi_{m}(x,\alpha)$), substituting in Eq. (42) and re-arranging terms, we get

[TABLE]

where we recall that

[TABLE]

Using the Parseval-Plancherel theorem, we get

[TABLE]

where the bar denotes the complex conjugate.

Substituting Eq. (43) in Eq. (45), and re-arranging terms, we find

[TABLE]

which defines $C_{nml}$ such that

[TABLE]

The final expression for $C_{nml}$ (Eq. 13) is then found using the binomial decompositions of $(m\alpha k-i)^{2m}$, $(l\beta k-i)^{2l}$ and $(n\gamma k-i)^{2n}$.

It can noted that $C_{nml}$ is a complex number, unless $u+v+w$ is even, such that $i^{u+v+w}=-1$. If $u+v+w$ is odd, then the integral is 0 (see A.2). Hence, we do not need to impose any restriction on $u+v+w$ for $C_{nml}$ to be real. We can also note that since $u\leqslant 2m$, $v\leqslant 2l$ and $w\leqslant 2n$, the integral in Eq. (13) converges.

We will then aim to look for an analytic expression for the integral in Eq. (13). We introduce

[TABLE]

A.2 Properties of $\mathcal{I}_{u,v,w}^{m,l,n}$’s integrand

In this section, we derive some properties of the integrand that appears in the definition of $\mathcal{I}_{u,v,w}^{m,l,n}$. Let us name it $f(x)$ (rigorously, we should write $f_{u,v,w}^{m,l,n}(x)$, but we drop the $u,v,w,n,m,l$ indices to simplify the notation), such that

[TABLE]

where $0\leqslant u\leqslant 2m$, $0\leqslant v\leqslant 2l$ and $0\leqslant w\leqslant 2n$ (see Eq. 13).

A.2.1 Alternative definition

When developping binomial expressions in $f$, we get another form for the function, that appears as the inverse of a series of $x$

[TABLE]

A.2.2 Parity

It is obvious that $f(-x)=(-1)^{u+v+w}f(x)$, so that $f$ is even (resp. odd) when $u+v+w$ is even (resp. odd). Hence, $\mathcal{I}_{u,v,w}^{m,l,n}$ vanishes when $u+v+w$ is odd, in which case $C_{nml}=0$. As mentioned in the main text, $C_{nml}$ is real if $u+v+w$ is even (then, $C_{nml}$ is real for all combinations of $u$, $v$, $w$, $n$, $m$, $l$). In the remainder of this appendix, we restrict ourselves to $x\geqslant 0$.

A.2.3 Value in 0

Evidently, $f(0)=0$ if $u+v+w\neq 0$. If $u+v+w=0$, then $f(0)=1$.

A.2.4 Limit at $x\rightarrow\infty$

Since $0\leqslant u\leqslant 2m$, $0\leqslant v\leqslant 2l$ and $0\leqslant w\leqslant 2n$, it is easy to see that $\lim_{x\rightarrow\infty}=0$.

A.2.5 Dependence on $n$, $m$, $l$

For a given $x$, $f$ quickly decreases as $m^{2(m+1)}$ (and similarly for $n$ and $l$), showing that only the contributions of low $n$, $m$ and $l$ are significant in $C_{nml}$. Then, the series $h_{n}$ (Eq. 47) converges quickly.

The left panel of Fig. 10 shows $f(x)$ for some realistic combinations $(\alpha,\beta,\gamma,n,ml,u,v,w)$. It can be seen that $f(x)$ is well-behaved, and tends quickly towards 0 for large $x$. Its extend depends on $\alpha$, $\beta$ and $\gamma$. As mentioned above, for a given $x$, it quickly decreases when $n$, $m$, $l$ increase. The right panel of Fig. 10 shows $f(x)$ when $\gamma=\alpha=\beta)=0.1$, $n$ = 2, $m$ = 1 $l$ = 1, $w$ = 2, $u=1$ and $v=1$, together with the values $x=\min\left((m\alpha)^{-1},(l\beta)^{-1},(n\gamma)^{-1}\right)$ and $x=\max\left((m\alpha)^{-1},(l\beta)^{-1},(n\gamma)^{-1}\right)$ (borders of the grey area), which are key values in computing $\mathcal{I}_{u,v,w}^{m,l,n}$ (see below).

A.3 Computation of the integral in $C_{nml}$

We now turn to look for an analytic expression for $I_{u,v,w}^{m,l,n}$. We first note that, due to the parity property of $f_{u,v,w}^{m,l,n}(x)$,

[TABLE]

if $u+v+w$ is even, or $I_{u,v,w}^{m,l,n}=0$ otherwise.

We decompose $I_{u,v,w,>}^{m,l,n}$ as

[TABLE]

where, as we show below, we impose $\mu<\min\left((m\alpha)^{-1},(l\beta)^{-1},(n\gamma)^{-1}\right)$ and $M>\max\left((m\alpha)^{-1},(l\beta)^{-1},(n\gamma)^{-1}\right)$.

A.3.1 Computation of first integral in r.h.s. of Eq. (52)

We first introduce

[TABLE]

where, for clarity, we wrote $f$ with all its indices.

Changing variable such that $k=y\mu$, we get

[TABLE]

Another change of variable, $t=y^{2}$, gives

[TABLE]

where we can recognize, if $\mu<\min\left((m\alpha)^{-1},(l\beta)^{-1},(n\gamma)^{-1}\right)$, Lauricella’s function of the fourth kind (e.g. Bezrodnykh 2016; Hasanov & Srivastava 2007)

[TABLE]

with $r=3$, $a=\frac{u+v+w+1}{2}$, $c=\frac{u+v+w+3}{2}$, $b_{1}=m+1$, $b_{2}=l+1$, $c_{2}=n+1$, $x_{1}=-(m\alpha\mu)^{2}$, $x_{2}=-(l\beta\mu)^{2}$ and $x_{3}=-(n\gamma\mu)^{2}$, such that

[TABLE]

$F_{D}^{(3)}$ is defined (as long as $\max(|x_{i}|)<1$) as a converging series (Eq. (1.4) in Hasanov & Srivastava 2007), such that

[TABLE]

where $(a)_{\nu}$ is the Pochhammer symbol, such that (after simplifying Pochhammer symbols):

[TABLE]

A.3.2 Computation of second integral in r.h.s. of Eq. (52)

We then introduce

[TABLE]

that we will compute in a similar fashion as the previous integral.

With successive changes of variables $y=k^{2}$, $z=y-\mu^{2}$, and $t=z/(M^{2}-\mu^{2})$ (under our assumptions, $M\neq\mu$; if $M=\mu$, $J_{u,v,w}^{m,l,n}(\mu,M)=0$), we obtain

[TABLE]

If $M\geqslant\sqrt{2}\mu$, then $J_{u,v,w}^{m,l,n}(\mu,M)$ must be calculated numerically. This is the case in the example of Fig. 10. However, if $M<\sqrt{2}\mu$ (i.e.  $\max\left\{(m\alpha)^{-1},(l\beta)^{-1},(n\gamma)^{-1}\right\}$ $<$ $\sqrt{2}\min\left\{(m\alpha)^{-1},(l\beta)^{-1},(n\gamma)^{-1}\right\}$, then we recognize a Lauricella $F_{D}^{(4)}$ function, such that

[TABLE]

If $M<\sqrt{2}\mu$, $J_{u,v,w}^{m,l,n}(\mu,M)$ can finally be expressed as a converging series

[TABLE]

A.3.3 Computation of third integral in r.h.s of Eq. (52)

We now introduce

[TABLE]

that we will compute in a similar fashion as $I_{u,v,w}^{m,l,n}$ and $J_{u,v,w}^{m,l,n}$.

We first note that $\int_{M}^{\infty}f_{u,v,w}^{m,l,n}(k){\rm d}k=\lim_{\xi\rightarrow\infty}I_{\xi}(M)$, where for simplicity, we droped the indices from the $K_{\xi}$ definition, and $\xi$ is some positive cut-off, and

[TABLE]

$K_{\xi}$ can be computed using successive changes of variable. First, we set $y=1/k$, so that (after rearranging some terms)

[TABLE]

Additional changes of variables ($z=y^{2}$, $t^{\prime}=z-1/\xi^{2}$, $t=\frac{M^{2}\xi^{2}}{\xi^{2}-M^{2}}t$) then yield

[TABLE]

where $g=3+n+m+l-(3+u+v+w)/2$. Taking the limit $\xi\rightarrow\infty$, we obtain

[TABLE]

Under the assumption that $M>\max\left((m\alpha)^{-1},(l\beta)^{-1},(n\gamma)^{-1}\right)$, we recognize a Lauricella $F_{D}^{(3)}$ function, such that

[TABLE]

and we can finally express $K_{u,v,w}^{m,l,n}(M)$ as a converging series

[TABLE]