# Stein’s method and approximating the multidimensional quantum harmonic
oscillator

Ian W. McKeague and Yvik Swan

###### Abstract

Stein’s method is used to study discrete representations of multidimensional
distributions that arise as approximations of states of quantum harmonic
oscillators. These representations model how quantum effects result from the
interaction of finitely many classical “worlds,” with the role of sample size
played by the number of worlds. Each approximation arises as the ground state
of a Hamiltonian involving a particular interworld potential function. Such
approximations have previously been studied for one-dimensional quantum
harmonic oscillators, but the multidimensional case has remained unresolved.
Our approach, framed in terms of spherical coordinates, provides the rate of
convergence of the discrete approximation to the ground state in terms of
Wasserstein distance. The fastest rate of convergence to the ground state is
found to occur in three dimensions. This result is obtained using a discrete
density approach to Stein’s method applied to the radial component of the
ground state solution.

Mathematics Subject Classification (MSC2020): 60B10, 81Q05, 81Q65
Keywords: Coupling, Hamiltonian, point configurations, Stein’s method,
Wasserstein distance

Author affiliations:
Ian W. McKeague: Columbia University, New York, U.S.A.
Yvik Swan: Université libre de Bruxelles, Bruxelles, Belgium.

Corresponding author: Yvik Swan

Acknowledgments. The authors thank Louis Chen, Davy Paindaveine, Erol Peköz,
Gesine Reinert and Thomas Verdebout for helpful discussions.

Declarations:
Funding: The research of Ian McKeague was partially supported by National
Science Foundation Grant DMS 2112938 and National Institutes of Health Grant
AG062401. The research of Yvik Swan was supported by by grant CDR/OL J.0197.20
from FRS-FNRS.
Availability of data and material: provided in the Supplementary material.
Code availability: provided in the Supplementary material.

###### Contents

  1. 1 Introduction
  2. 2 Extending interworld potentials to higher dimensions
    1. 2.1 Two-dimensional case
    2. 2.2 Interworld potentials for $d\geq 3$
    3. 2.3 Interworld potentials for excited states
    4. 2.4 Wasserstein distance in spherical coordinates
  3. 3 Stein’s method for approximating radial distributions
  4. 4 Convergence of multidimensional ground states
  5. A Supplementary material
    1. A.1 Preliminary remarks on Stein operators
    2. A.2 Stein’s method for radial distributions
      1. A.2.1 Notations and background
      2. A.2.2 When the base distribution is standard Gaussian
      3. A.2.3 When $k=0,1$ and a proof of Proposition 3.3
      4. A.2.4 When $k\geq 2$
    3. A.3 Upper bounds on the rate of convergence via coupling
  6. B Codes
    1. B.1 Mathematica code
    2. B.2 R Code

##  1 Introduction

The many interacting worlds (MIW) approach to quantum mechanics due to Hall et
al. [9] posits a Hamiltonian for a one-dimensional harmonic oscillator of the
form

$H_{1}({\bf x},{\bf p})=E({\bf p})+V_{1}({\bf x})+U_{1}({\bf x}),$

where the locations of particles (having unit mass) in $N$ worlds are
specified by ${\bf x}=(x_{1},\ldots,x_{N})$ with $x_{1}>x_{2}>\ldots>x_{N}$,
and their momenta by ${\bf p}=(p_{1},\ldots,p_{N})$. Here $E({\bf
p})=\sum_{n=1}^{N}p_{n}^{2}/2$ is the kinetic energy, $V_{1}({\bf
x})=\sum_{n=1}^{N}x_{n}^{2}$ is the potential energy (for the parabolic trap),
and

$U_{1}({\bf x})=\sum_{n=1}^{N}\left({1\over x_{n+1}-x_{n}}-{1\over
x_{n}-x_{n-1}}\right)^{2}$

is called the “interworld” potential, where $x_{0}=\infty$ and
$x_{N+1}=-\infty$.

The ground state particle locations $x_{n}$ are symmetric (about the origin)
and satisfy [9]

$x_{n+1}=x_{n}-\frac{1}{x_{1}+\cdots+x_{n}}.$ (1)

McKeague and Levin [15] showed that the empirical distribution of the
resulting ground state solution $\\{x_{n},n=1,\ldots,N\\}$ tends to standard
Gaussian, and conjectured that the optimal rate of convergence in Wasserstein
distance is of order $\sqrt{\log N}/N$; this conjecture was recently proved by
Chen and Thành [3].

For the first excited state, McKeague, Peköz and Swan [16] showed that an
extension of the above interworld potential leads to the two-sided Maxwell
distribution as the limit (agreeing with the classical quantum harmonic
oscillator). Ghadimi et al. [7] studied non-locality by introducing other
extensions of the MIW interworld potential for the first excited state using
higher-order smoothing methods.

In this article we study the question of how the MIW approach can be extended
to general $d$-dimensional settings. Herrmann et al. [10] examined the case
$d=2$ and proposed using Delaunay triangulations and Voronoi tesselations of
point configurations to extend the notion of the interworld potential. They
developed a numerical algorithm to estimate the resulting ground state
configuration, but the question of whether it is possible to establish an
asymptotically valid approximation of such a solution to the classical quantum
harmonic oscillator ground state solution was not addressed.

Our contribution is two-fold: (1) we introduce new interworld potentials that
apply in general $d$-dimensional settings and that lead to tractable ground
state solutions (as well as some excited state solutions) for the
corresponding MIW Hamiltonian; (2) we provide a new version the discrete
density approach to Stein’s method that furnishes upper bounds on the
convergence rate of the radial components of these ground state solutions.
Evaluating these bounds numerically, leads to optimal rates of convergence (in
terms of Wasserstein distance) of the full MIW configurations that apply to
general multidimensional settings. In contrast to [15, 16], which rely on the
coupling version of Stein’s method, the present approach leads to optimal
rates. In the Supplementary material we explore the types of upper bounds that
can be obtained using the coupling approach in the multidimensional case.

The paper is organized as follows. Section 2 develops the proposed interworld
potentials, first for the two- and three-dimensional cases, then for the
general $d$-dimensional setting, including some that apply to excited states.
The resulting MIW ground state solutions are described in terms of their
spherical coordinates, and plots of the solutions are provided for $d=2$ and
$3$. An upper bound on the Wasserstein distance between two probability
measures that have independent radial and directional components is provided
at the end of Section 2. Section 3 restricts attention to the radial component
and develops the discrete density approach to Stein’s method mentioned above.
In Section 4 we wrap up by providing optimal rates for full $d$-dimensional
ground state solutions. The Supplementary material includes discussion of the
coupling approach mentioned above along with computer code.

##  2 Extending interworld potentials to higher dimensions

In this section we introduce interworld potentials that lead to satisfactory
MIW approximations to the the ground states and excited states of the
$d$-dimensional isotropic quantum harmonic oscillator. First we consider the
two-dimensional case.

###  2.1 Two-dimensional case

When $d=2$, the idea is to couch the problem in terms of polar coordinates,
which allows the interworld potential to adapt to the basic geometry of the
problem, specifically via a separation of its angular (directional) and radial
components. To that end, a configuration ${\bf x}$ of points in
$\mathbb{R}^{2}$ is specified in terms of (signed) polar coordinates as

$\\{(r_{nj},\theta_{j})\colon\ n=1,\ldots,N_{j},\ j=1,\ldots,M\\},$

where there are $M=M({\bf x})\geq 1$ (distinct) angular coordinate values
satisfying $0\leq\theta_{1}<\theta_{2}<\ldots<\theta_{M}<\pi$. The number of
points $N_{j}$ in any direction $\theta_{j}$ is assumed to satisfy $N_{j}\geq
2$, to avoid the possibility of a “degenerate” radial contribution from a
single point at the origin. The radial coordinates $r_{nj}$ in direction
$\theta_{j}$ in this representation are signed (can be negative as well as
positive) and satisfy $r_{1j}>r_{2j}>\ldots>r_{N_{j}j}$; a point $(r,\theta)$
with $r<0$ is understood to correspond to the point with polar coordinates
$(|r|,\theta+\pi)$. The total number of points is $N=N_{1}+\ldots+N_{M}$, with
the understanding that points at the origin arising from different directions
are distinct elements of the configuration.

We propose the following ansatz for the interworld potential:

$\displaystyle U_{2}({\bf x})$
$\displaystyle=4\sum_{j=1}^{M}\sum_{n=1}^{N_{j}}\left[{1\over
R_{2}(r_{n+1,j})-R_{2}(r_{nj})}-{1\over
R_{2}(r_{nj})-R_{2}(r_{n-1,j})}\right]^{2}r_{nj}^{2}$ $\displaystyle\ \ \
+\pi\sum_{j\sim k}{1\over|\theta_{j}-\theta_{k}|_{\pi}}+{N^{2}\over
M^{2}}L\left(\sum_{j\sim k}|N_{j}-N_{k}|\right),$ (2)

where the notation $R_{d}(r)=|r|^{d}\,{\rm sign}(r)$ for general integers
$d\geq 2$, and $L(x)=\max(1,x/2)$, will be used frequently in the sequel. By
convention we set $R_{d}(r_{0j})=\infty$ and $R_{d}(r_{N_{j}+1,j})=-\infty$.
The summation in the second and third terms is over all pairs of angular
values that are neighbors (mod $\pi$), and $|\cdot|_{\pi}$ denotes absolute
value mod $\pi$; in the sequel we also use such notation when $\pi$ is
replaced by other positive real numbers, depending on the context.

The intuition behind the first term in this ansatz comes from the well-known
representation of the standard Gaussian density in two-dimensions in terms of
polar coordinates. That is, the angular coordinate (in the sense we defined
above) can be taken as uniformly distributed on $[0,\pi]$, independent of the
(signed) radial coordinate, which has the (two-sided) Rayleigh density
$p(r)=b(r)\varphi(r)$, where $\varphi$ is a standard normal density and
$b(r)=\sqrt{\pi/2}|r|$, $r\in\mathbb{R}$.

The first term in $U_{2}({\bf x})$ is based on a derivation given in [16], in
which an interworld potential is proposed for densities of the general form
$p(x)=b(x)\varphi(x)$, $x\in\mathbb{R}$. The focus in that paper was on the
case of the two-sided Maxwell distribution for which $b(x)=x^{2}$. The general
ansatz for the interworld potential of a one-dimensional $N$-point
configuration $x_{1}>x_{2}>\ldots>x_{N}$ was proposed to be

$U_{b}({\bf x})=\sum_{n=1}^{N}\left[{1\over B(x_{n+1})-B(x_{n})}-{1\over
B(x_{n})-B(x_{n-1})}\right]^{2}b(x_{n})^{2},$ (3)

where $B(x)=\int_{0}^{x}b(t)\,dt$, $B(x_{0})=\infty$ and $B(x_{N+1})=-\infty$.
When $b(x)$ is proportional to $|x|^{k}$ for some positive integer $k$, as for
the Rayleigh density ($k=1$), it can be shown that the minimizer of the
Hamiltonian $U_{b}({\bf x})+V_{1}({\bf x})$ is a symmetric solution of the
recursion

$B(x_{n+1})=B(x_{n})-\left(\sum_{i=1}^{n}\frac{x_{i}}{b(x_{i})}\right)^{-1}.$
(4)

For the Rayleigh distribution, $b(x)=\sqrt{\pi/2}|x|$, so the interworld
potential (3) becomes

$U_{b}({\bf x})=4\sum_{n=1}^{N}\left[{1\over
R_{2}(x_{n+1})-R_{2}(x_{n})}-{1\over
R_{2}(x_{n})-R_{2}(x_{n-1})}\right]^{2}x_{n}^{2}.$

From (4), the ground state is then a symmetric solution to the recursion
equation

$R_{2}(x_{n+1})=R_{2}(x_{n})-2\left(\sum_{i=1}^{n}{\rm
sign}(x_{i})\right)^{-1}.$ (5)

Further, it follows from analogous arguments in [16, Section 2] that the
ground state value of the Hamiltonian in this case is $4(N-1)$. Fig. 1
compares the empirical distribution of the solution of the above recursion
(for $N=22$ points) with the two-sided Rayleigh density, showing that the
agreement is remarkably accurate.

Figure 1: Histogram of the symmetric solution of the recursion (5) for $N=22$
compared with the two-sided Rayleigh density; the breaks in the histogram are
$x_{1},\ldots,x_{N}$

Returning now to the two-dimensional setting with Hamiltonian

$H_{2}({\bf x})=U_{2}({\bf x})+\sum_{j=1}^{M}\sum_{n=1}^{N_{j}}r_{nj}^{2},$
(6)

note that the $M$ components of the first term of the interworld potential
$U_{2}({\bf x})$ can be combined with the corresponding potential energy terms
and separately minimized (using the recursion (5)), giving a combined
contribution of

$\sum_{j=1}^{M}4(N_{j}-1)=4(N-M)$

to the Hamiltonian.

It remains to minimize the sum of the second and third terms in $U_{2}({\bf
x})$ to furnish the complete ground state solution. First consider $M$ as
fixed and note that by Cauchy–Schwarz

$\displaystyle M$ $\displaystyle=\sum_{j\sim
k}{|\theta_{j}-\theta_{k}|_{\pi}^{1/2}\over|\theta_{j}-\theta_{k}|_{\pi}^{1/2}}\leq\left(\sum_{j\sim
k}{|\theta_{j}-\theta_{k}|_{\pi}}\right)^{1/2}\left(\sum_{j\sim
k}{1\over|\theta_{j}-\theta_{k}|_{\pi}}\right)^{1/2}$
$\displaystyle=\left(\pi\sum_{j\sim
k}{1\over|\theta_{j}-\theta_{k}|_{\pi}}\right)^{1/2},$

so the minimum of the second term is $M^{2}$, attained by setting
$\theta_{j}=(j-1)\pi/M,\ j=1,\ldots,M$, because that produces equality
throughout the above display. With $M\leq N/2$ fixed, the third term in
$U_{2}({\bf x})$ is minimized by distributing the $N$ points as evenly as
possible among the $M$ directions (with at least two in each direction) and
then allotting the remaining ($N$ mod $M$) points to a sequence of neighboring
directions. This results in $N_{j}=N_{k}$ for all neighbors $j$ and $k$ except
possibly for two pairs of neighbors in which $|N_{j}-N_{k}|=1$.

The sum of the last two terms of $U_{2}({\bf x})$ therefore has the minimal
value $M^{2}+{N^{2}/M^{2}}$, and when combined with the contribution $4(N-M)$
arising from the first term in $U_{2}({\bf x})$ along with the potential
energy, as discussed earlier, the ground state minimizes
$4(N-M)+M^{2}+N^{2}/M^{2}$ as a function of $M$. Asymptotically, this minimum
is attained when $M\sim\sqrt{N}$. It follows also that the empirical
distribution of $\\{N_{j}/N,j=1,\ldots,M\\}$ in the ground state converges
weakly to uniform on $(0,1)$ as $N\to\infty$.

Figure 2: Ground states: $d=2$ (left), $N=22^{2}=484$ points, $M=22$
directions and $N_{j}=22$ points in each direction; $d=3$ (right), $N=2744$
points, $M=7$ shells (polar angles in $[0,\pi/2]$), $K_{j}=28$ azimuthal
angles in each shell, and $N_{jk}=14$ points in each radial direction

Combined with the fact that the empirical distribution (conditional on $M$) of
any set of minimizing angular coordinates $\\{\theta_{j},j=1,\ldots,M\\}$
converges weakly to uniform on $(0,\pi)$ as $M\to\infty$, we conclude that the
empirical distribution of the angular coordinates of the $N$ points in the
ground state has the same limit as $N\to\infty$. Thus, provided we can show
that the empirical distribution of the radial coordinates of the ground state
converges to the Rayleigh distribution, it will follow that the full empirical
distribution of the ground state (as illustrated by the left panel of Fig. 2)
converges to the standard 2D Gaussian distribution.

###  2.2 Interworld potentials for $d\geq 3$

In the three-dimensional case, a configuration ${\bf x}$ of points in
$\mathbb{R}^{3}$ is specified in terms of (signed) spherical coordinates as

$\displaystyle\\{(r_{njk},\theta_{j},\phi_{jk})\colon\ n=1,\ldots,N_{jk},\
j=1,\ldots,M,\ k=1,\ldots,K_{j}\\},$ (7)

where there are (distinct) directions corresponding to points with polar
angles $0\leq\theta_{1}<\theta_{2}<\ldots<\theta_{M}\leq\pi/2$ and azimuthal
angles $0\leq\phi_{j1}<\phi_{j2}<\ldots<\phi_{jK_{j}}<2\pi$, along with
$N_{jk}\geq 2$ points in each direction $(\theta_{j},\phi_{jk})$ determined by
their (signed) radial distances $r_{1jk}<\ldots<r_{N_{jk}jk}$. There are a
total of $N_{j\boldsymbol{\cdot}}=\sum_{k=1}^{K_{j}}N_{jk}$ points
corresponding to the $K_{j}$ directions in the $j$th “shell” with polar angle
$\theta_{j}$. Over all $M$ shells, there are a total of
$N=\sum_{j=1}^{M}N_{j\boldsymbol{\cdot}}$ points arising from
$K=\sum_{j=1}^{M}K_{j}$ different directions.

The spherical coordinates of a 3-dimensional standard Gaussian random vector
are independent. The signed radial component has the two-sided Maxwell
distribution with density $r^{2}\varphi(r),r\in\mathbb{R}$, and the azimuthal
angle is uniformly distributed on $[0,2\pi)$. The polar angle has cdf
$1-\cos(\theta)$, $\theta\in[0,\pi/2]$. Motivated by similar considerations as
in the two-dimensional case, the proposed interworld potential is taken to be

$\displaystyle U_{3}({\bf x})$
$\displaystyle=9\sum_{j=1}^{M}\sum_{k=1}^{K_{j}}\sum_{n=1}^{N_{jk}}\left({1\over
r_{n+1,jk}^{3}-r_{njk}^{3}}-{1\over
r_{njk}^{3}-r_{n-1,jk}^{3}}\right)^{2}r_{njk}^{4}$ $\displaystyle\ \ \
+\sum_{i\sim j}{1\over|\cos(\theta_{i})-\cos(\theta_{j})|}+{\pi\over
2}\sum_{j=1}^{M}\sum_{k\sim l}{1\over|\phi_{jk}-\phi_{jl}|_{2\pi}}$ (8)
$\displaystyle\ \ \ +\sum_{j=1}^{M}{N_{j\boldsymbol{\cdot}}^{2}\over
K_{j}^{2}}\,L\\!\\!\left(\sum_{k\sim l}|N_{jk}-N_{jl}|\right)+{N\over
4M}\,L\\!\\!\left(\sum_{i\sim
j}|N_{i\boldsymbol{\cdot}}-N_{j\boldsymbol{\cdot}}|\right).$

The last term above has a distinctive form that is not present in $U_{2}({\bf
x})$. The denominator in this term is designed to control the number of polar
angles ($M$) in the second term, which in turn needs to be balanced with the
$K_{j}$, since the polar angles range over an interval only one fourth the
length of that for the azimuthal angles. Numerical experiments show that the
resulting directional component in the ground state is close to uniformly
distributed over the sphere, as illustrated in the right panel of Fig. 2.

We now seek to minimize the Hamiltonian

$\displaystyle H_{3}({\bf x})=U_{3}({\bf
x})+\sum_{j=1}^{M}\sum_{k=1}^{K_{j}}\sum_{n=1}^{N_{jk}}r_{njk}^{2}.$ (9)

The components with distinct values of $(j,k)$ in the first term of the
interworld potential $U_{3}({\bf x})$ can be combined with the corresponding
potential energy terms and separately minimized (in terms of the recursion (4)
with $b(x)=x^{2}$ that was studied in [16]), giving a combined contribution of

$\sum_{j=1}^{M}\sum_{k=1}^{K_{j}}6(N_{jk}-1)=6(N-K)$

to the Hamiltonian.

The next step is to minimize each of the last three terms in (2.2) for fixed
$M$, $K_{1},\ldots,K_{M}$. Using a similar Cauchy–Schwarz argument to what we
used to minimize the second term in (2.1), the second term in (2.2) is
minimized for

$\displaystyle\theta_{j}=\cos^{-1}(1-(j-1)/(M-1)),\ j=1,\ldots,M,$ (10)

and the minimum is $(M-1)^{2}$. Similarly, the minimum of the third term in
(2.2) is attained by setting

$\phi_{jk}=2\pi(k-1)/K_{j},\ k=1,\ldots,K_{j},\ j=1,\ldots,M,$

and the minimum is $\sum_{j=1}^{M}K_{j}^{2}/4$. The minimum of the fourth term
given fixed values of the $K_{j}$ and $N_{j\boldsymbol{\cdot}}$ is found using
the same argument as in the two-dimensional case for the third term in (2.1),
resulting in the minimum
$\sum_{j=1}^{M}{N_{j\boldsymbol{\cdot}}^{2}/K_{j}^{2}}$ being attained when
$N_{jk}\sim N_{j\boldsymbol{\cdot}}/K_{j}$ for $k=1,\ldots,K_{j}$. This
reduces the minimal value of the Hamiltonian to

$6(N-K)+(M-1)^{2}+\sum_{j=1}^{M}[K_{j}^{2}/4+{N_{j\boldsymbol{\cdot}}^{2}/K_{j}^{2}}]+{N\over
4M}\,L\\!\\!\left(\sum_{i\sim
j}|N_{i\boldsymbol{\cdot}}-N_{j\boldsymbol{\cdot}}|\right).$

For fixed $M$, this expression is minimized when
$K_{j}\sim\sqrt{2N_{j\boldsymbol{\cdot}}}$ as
$N_{j\boldsymbol{\cdot}}\to\infty$, and it remains to minimize

$7N+(M-1)^{2}+{N\over 4M}\,L\\!\\!\left(\sum_{i\sim
j}|N_{i\boldsymbol{\cdot}}-N_{j\boldsymbol{\cdot}}|\right).$

Asymptotically, this expression is minimized by $N_{j\boldsymbol{\cdot}}\sim
N/M$ and setting the number of shells $M\sim N^{1/3}/2$. This implies

$N_{jk}\sim
N_{j\boldsymbol{\cdot}}/K_{j}\sim\sqrt{N_{j\boldsymbol{\cdot}}/2}\sim\sqrt{N/(2M)}\sim
N^{1/3},\ \ K_{j}\sim 2N^{1/3}.$

We conclude that the ground state asymptotically consists of $K_{j}=2N^{1/3}$
directions in each of $M=N^{1/3}/2$ shells, with $N_{jk}=N^{1/3}$ points
falling in each direction. An example with $N=2744$ points is provided in the
right panel of Fig. 2.

The extension to general $d\geq 3$ is straightforward. The radial component of
the $d$-dimensional standard Gaussian distribution has density
$p(r)=c_{d}|r|^{d-1}\varphi(r),\ r\in\mathbb{R}$, where $c_{d}$ is a
normalizing constant, and is independent of the (signed) directional
component, which is uniformly distributed (in the sense of Hausdorff measure)
over the unit-hemisphere $\\{x\in\mathbb{R}^{d}\colon{x}_{1}\geq
0,|{x}|=1\\}$. The configuration $\bf x$ is now specified using signed
hyperspherical coordinates, precisely as in (7), except each polar angle
$\theta_{j}$ is now required to be a vector $\theta_{j}=(\theta_{jl})$ of
$d-2$ polar angles, with each component satisfying
$0\leq\theta_{1l}<\theta_{2l}<\ldots<\theta_{Ml}\leq\pi/2$, $l=1,\ldots,d-2$.
The notation for the azimuthal angles, the $M$ shells and the $K$ directions
remain the same.

The $d-2$ polar angles are iid under the $d$-dimensional standard Gaussian,
and their cdf can be expressed using a simple formula for the area of a
hyperspherical cap [13]:

$F_{d}(\theta)=B\\!\left(\sin^{2}\theta,{d-1\over 2},{1\over 2}\right),\
\theta\in[0,\pi/2],$

where $B(x;a,b)=\int_{0}^{x}t^{a-1}(1-t)^{b-1}\,dt$ is the incomplete Beta
function. In particular, $F_{4}(\theta)=(2\theta-\sin(2\theta))/\pi$ and
$F_{5}(\theta)=1-(3/2)\cos\theta+(1/2)\cos^{3}\theta$. This leads to the
interworld potential

$\displaystyle U_{d}({\bf x})=$ $\displaystyle
d^{2}\sum_{j=1}^{M}\sum_{k=1}^{K_{j}}\sum_{n=1}^{N_{jk}}\left({1\over
R_{d}(r_{n+1,jk})-R_{d}(r_{njk})}-{1\over
R_{d}(r_{njk})-R_{d}(r_{n-1,jk})}\right)^{2}r_{njk}^{2(d-1)}$ $\displaystyle\
\ \ +{1\over(d-2)}\sum_{l=1}^{d-2}\sum_{i\sim
j}{1\over|F_{d}(\theta_{il})-F_{d}(\theta_{jl})|}+{\pi\over
2}\sum_{j=1}^{M}\sum_{k\sim l}{1\over|\phi_{jk}-\phi_{jl}|_{2\pi}}$ (11)
$\displaystyle\ \ \ +\sum_{j=1}^{M}{N_{j\boldsymbol{\cdot}}^{2}\over
K_{j}^{2}}\,L\\!\\!\left(\sum_{k\sim l}|N_{jk}-N_{jl}|\right)+{N\over
4(d-2)M^{d-2}}\,L\\!\\!\left(\sum_{i\sim
j}|N_{i\boldsymbol{\cdot}}-N_{j\boldsymbol{\cdot}}|\right).$

It follows by a routine extension of the argument we used in the $d=3$ case
that the ground state asymptotically consists of $K_{j}=2N^{1/2-1/(2d)}$
directions in each of $M=N^{1/d}/2$ shells, with $N_{jk}=N^{1/2-1/(2d)}$
points falling in each direction.

A possible alternative approach is to specify the directions in the unit-
hemisphere by a minimal Riesz energy point configuration, without reference to
polar and azimuthal components. As surveyed in [1], such configurations can
provide asymptotically uniformly distributed point sets with respect to
surface area (Hausdorff measure), with important applications in quasi-Monte
Carlo, approximation theory and material physics. However, there does not seem
to be a way of linking the numbers of radial points with the directions
obtained from minimizing Riesz energy that would lead to a Gaussian
approximation. In contrast, the Hamiltonian based on the interworld potential
(2.2) is readily minimized by following the same argument we used to minimize
(2.2); the only difference in the solution is that $\cos^{-1}$ in (10) is
replaced by $F_{d}^{-1}$ in the expression of each polar angle $\theta_{jl}$.
Moreover, as explained in the following sections, our proposed approach leads
to explicit rates of convergence (in terms of Wasserstein distance) of the
empirical measure of the ground state to standard Gaussian.

###  2.3 Interworld potentials for excited states

The eigenstates of the classical $d$-dimensional isotropic quantum harmonic
oscillator consist of products of $d$ one-dimensional eigenfunctions, with
separate euclidean coordinates in each component. In this section we show how
the approach in the previous sections extends naturally to the case when some
of these one-dimensional components are in excited (higher energy) states. The
various eigenstates of the full system are described by vectors of quantum
numbers indicating the energy level of each component (with 0 indicating the
ground state). When expressed in spherical or polar coordinates, there is a
separation of variables in the various eigenfunctions, which implies that the
corresponding densities have independent contributions from the radial and
angular components. This allows a separation of the radial and angular
components in the interworld potential, as we now make explicit in examples of
excited states for $d=2$ and $3$. The idea readily extends to any excited
state of the MIW quantum harmonic oscillator.

For $d=2$, the lowest three excited states have quantum numbers $(1,0)$,
$(0,1)$ and $(1,1)$. The density of the $(1,0)$-state in signed polar
coordinates is proportional to $\cos^{2}(\theta)|r|^{3}\varphi(r)$,
$\theta\in[0,\pi)$, $r\in\mathbb{R}$. The cdf of the angular component is
$G_{2}(\theta)=(\theta+\sin(\theta)\cos(\theta))/\pi$, and the density of the
signed radial component is proportional to $|r|^{3}\varphi(r)$. This leads to
the following interworld potential similar to (2.1):

$\displaystyle U_{2}({\bf x})$
$\displaystyle=16\sum_{j=1}^{M}\sum_{n=1}^{N_{j}}\left[{1\over
R_{4}(r_{n+1,j})-R_{4}(r_{nj})}-{1\over
R_{4}(r_{nj})-R_{4}(r_{n-1,j})}\right]^{2}r_{nj}^{6}$ $\displaystyle\ \ \
+\sum_{j\sim k}{1\over|G_{2}(\theta_{j})-G_{2}(\theta_{k})|_{1}}+{N^{2}\over
M^{2}}L\left(\sum_{j\sim k}|N_{j}-N_{k}|\right).$ (12)

Following similar arguments to those used earlier, the Hamiltonian (6) based
on (2.3) is minimized with the radial coordinates given by the symmetric
solution to the recursion (4) with $b(x)=|x|^{3}$, and setting
$\theta_{j}=G_{2}^{-1}((j-1)/M)$ with $M\sim\sqrt{N}$. Fig. 3 shows the
$(1,0)$ and $(1,1)$ excited states resulting from $N=484$ points. The $(0,1)$
state is the same as $(1,0)$ except rotated by 90 degrees.

Figure 3: Excited states $(1,0)$ (left), $(1,1)$ (right), for $d=2$,
$N=22^{2}=484$ points, $M=22$, $N_{j}=22$, cf. the ground state in Fig. 2

For $d=3$, the density of the $(1,0,0)$ excited state in signed spherical
coordinates is proportional to
$\sin^{3}(\theta)\cos^{2}(\phi)r^{4}\varphi(r)$, $\theta\in[0,\pi/2)$,
$\phi\in[0,2\pi)$, $r\in\mathbb{R}$. The cdf of the polar angle is
$G_{3}(\theta)=(\cos(3\theta)-9\cos(\theta))/8$, and the cdf of the azimuthal
angle is $A(\phi)=(\phi-\sin(\phi)\cos(\phi))/(2\pi)$. The interworld
potential is similar to (2.2):

$\displaystyle U_{3}({\bf x})$
$\displaystyle=25\sum_{j=1}^{M}\sum_{k=1}^{K_{j}}\sum_{n=1}^{N_{jk}}\left({1\over
r_{n+1,jk}^{5}-r_{njk}^{5}}-{1\over
r_{njk}^{5}-r_{n-1,jk}^{5}}\right)^{2}r_{njk}^{8}$ $\displaystyle\ \ \
+\sum_{i\sim j}{1\over|G_{3}(\theta_{i})-G_{3}(\theta_{j})|}+{\pi\over
2}\sum_{j=1}^{M}\sum_{k\sim l}{1\over|A(\phi_{jk})-A(\phi_{jl})|_{2\pi}}$
$\displaystyle\ \ \ +\sum_{j=1}^{M}{N_{j\boldsymbol{\cdot}}^{2}\over
K_{j}^{2}}\,L\\!\\!\left(\sum_{k\sim l}|N_{jk}-N_{jl}|\right)+{N\over
4M}\,L\\!\\!\left(\sum_{i\sim
j}|N_{i\boldsymbol{\cdot}}-N_{j\boldsymbol{\cdot}}|\right).$

The configurations of the $(1,0,0)$, $(1,1,0)$ and $(1,1,1)$ excited states
resulting from $N=2744$ points are displayed in Fig. 4.

Figure 4: Excited states $(1,0,0)$ (left), $(1,1,0)$ (middle) and $(1,1,1)$
(right) for $d=3$, $N=2744$ points, $M=7$, $K_{j}=28$, $N_{jk}=14$, cf. the
ground state in the right panel of Fig. 2

###  2.4 Wasserstein distance in spherical coordinates

The Wasserstein distance between probability measures $\mu$ and $\nu$ on a
Polish metric space $(S,\rho)$ can be defined equivalently as

$\mathop{d_{\mathrm{W}}}(\mu,\nu)=\sup_{h\in\mathcal{H}}\left|\int_{S}h\,d(\mu-\nu)\right|=\inf\mathbb{E}\rho(X,Y),$

where $\mathcal{H}$ is the family of Lipschitz functions $h\colon
S\rightarrow\mathbb{R}$ with $\mbox{Lip}(h)\leq 1$, and the (attained) infimum
is over all coupled $S$-valued random elements $X\sim\mu$ and $Y\sim\nu$. When
$S=\mathbb{R}$ and $\mu$ and $\nu$ have cdfs $F$ and $G$, their Wasserstein
distance coincides with the $L_{1}$-distance $\int|F(x)-G(x)|\,\mathrm{d}x$.
The following inequality bounds the Wasserstein distance between two product
measures $\mu=\prod_{i=1}^{d}\mu_{i}$ and $\nu=\prod_{i=1}^{d}\nu_{i}$ on
$\mathbb{R}^{d}$ endowed with the euclidean metric:

$\displaystyle\mathop{d_{\mathrm{W}}}(\mu,\nu)\leq\sum_{i=1}^{d}\mathop{d_{\mathrm{W}}}(\mu_{i},\nu_{i}),$

see [14, Lemma 3].

We now present a variation of this result to enable the study of the
convergence of the Wasserstein distance between the empirical distribution of
the energy minimizing configuration ${\bf x}$ and the distribution specified
by quantum theory. Concentrating on the case $d=3$ for simplicity, the
following result provides a bound on the Wasserstein distance between two
probability measures $\mu$ and $\nu$ on $\mathbb{R}^{3}$ (of the type that is
relevant here) in terms of the Wasserstein distance between the distributions
of their respective (signed) spherical coordinates (denoted $\mu_{r}$,
$\mu_{\theta}$, $\mu_{\phi}$, etc).

######  Lemma 2.1.

Let $X\sim\mu$ and $Y\sim\nu$ be random vectors in $\mathbb{R}^{3}$. If the
signed spherical coordinates of $X$ are independent, and the same holds for
$Y$, then

$\displaystyle\mathop{d_{\mathrm{W}}}(\mu,\nu)\leq\mathop{d_{\mathrm{W}}}(\mu_{r},\nu_{r})+\sqrt{m_{\mu}m_{\nu}}[\mathop{d_{\mathrm{W}}}(\mu_{\theta},\nu_{\theta})+\mathop{d_{\mathrm{W}}}(\mu_{\phi},\nu_{\phi})],$

where $m_{\mu}=\int|x|\,d\mu_{r}(x)$ is the mean absolute deviation of the
radial coordinate of $X$, similarly for $m_{\nu}$.

The result is the same in the case $d=2$, except the azimuthal contribution
$\mathop{d_{\mathrm{W}}}(\mu_{\phi},\nu_{\phi})$ is absent. The result
naturally extends to general $d\geq 3$, with additional terms of the form
$\mathop{d_{\mathrm{W}}}(\mu_{\theta},\nu_{\theta})$ representing each of the
($d-2$) polar angles. The directional contributions of the lowest energy
configuration will be shown to have a faster rate of convergence (as
$N\to\infty$) than the radial contribution, so it will suffice to restrict
attention to the convergence of the latter, as we do in the next section.

###### Proof.

For points $\underline{x}_{1},\underline{x}_{2}\in\mathbb{R}^{3}$, the
squared-euclidean distance between them in terms of their signed spherical
coordinates $(r_{i},\theta_{i},\phi_{i})$, $i=1,2$, is

$\displaystyle\|\underline{x}_{1}-\underline{x}_{2}\|^{2}$
$\displaystyle=r_{1}^{2}+r_{2}^{2}-2r_{1}r_{2}[\cos(\theta_{1}-\theta_{2})+\sin(\theta_{1})\sin(\theta_{2})(\cos(\phi_{1}-\phi_{2})-1)]$
$\displaystyle=(r_{1}-r_{2})^{2}-2r_{1}r_{2}[\cos(\theta_{1}-\theta_{2})-1+\sin(\theta_{1})\sin(\theta_{2})(\cos(\phi_{1}-\phi_{2})-1)]$
$\displaystyle\leq(r_{1}-r_{2})^{2}+r_{1}r_{2}[(\theta_{1}-\theta_{2})^{2}+(\phi_{1}-\phi_{2})^{2}],$

using the inequality $|\\!\cos(x)-1|\leq x^{2}/2$ for all $x\in\mathbb{R}$.
Then using the inequality
$(|a|+|b|+|c|)^{1/2}\leq|a|^{1/2}+|b|^{1/2}+|c|^{1/2}$ for any
$a,b,c\in\mathbb{R}$, we obtain

$\|\underline{x}_{1}-\underline{x}_{2}\|\leq|r_{1}-r_{2}|+\sqrt{|r_{1}||r_{2}|}[|\theta_{1}-\theta_{2}|+|\phi_{1}-\phi_{2}|].$
(13)

We can choose coupled random vectors $X^{\star}\sim\mu$ and $Y^{\star}\sim\nu$
such that

$\mathop{d_{\mathrm{W}}}(\mu_{r},\nu_{r})=\mathbb{E}|X_{r}^{\star}-Y_{r}^{\star}|,\
\
\mathop{d_{\mathrm{W}}}(\mu_{\theta},\nu_{\theta})=\mathbb{E}|X_{\theta}^{\star}-Y_{\theta}^{\star}|,\
\
\mathop{d_{\mathrm{W}}}(\mu_{\phi},\nu_{\phi})=\mathbb{E}|X_{\phi}^{\star}-Y_{\phi}^{\star}|,$

with the random vectors $(X_{r}^{\star},Y_{r}^{\star})$,
$(X_{\theta}^{\star},Y_{\theta}^{\star})$ and
$(X_{\phi}^{\star},Y_{\phi}^{\star})$ being independent. Substituting
$X^{\star}$ and $Y^{\star}$ for $\underline{x}_{1}$ and $\underline{x}_{2}$ in
(13) and taking expectations of both sides leads to

$\mathop{d_{\mathrm{W}}}(\mu,\nu)\leq\mathbb{E}\|X^{\star}-Y^{\star}\|\leq\mathbb{E}|X_{r}^{\star}-Y_{r}^{\star}|+[\mathbb{E}\sqrt{|X_{r}^{\star}||Y_{r}^{\star}}|][\mathbb{E}|X_{\theta}^{\star}-Y_{\theta}^{\star}|+\mathbb{E}|X_{\phi}^{\star}-Y_{\phi}^{\star}|].$

The proof is completed using Cauchy–Schwarz to obtain
$\mathbb{E}\sqrt{{|X_{r}^{\star}||Y_{r}^{\star}}|}\leq\sqrt{m_{\mu}m_{\nu}}$.
∎

##  3 Stein’s method for approximating radial distributions

We begin this section with some generalities concerning the version of Stein’s
method developed in [6, 5]. We only give a summary here and refer the reader
unfamiliar with Stein’s method to the Supplementary material for more
information and background.

Consider a standardized (i.e., mean 0, variance 1) random variable
$F_{\infty}$ with density $p_{\infty}$ (in the sequel $p_{\infty}$ will take
the form of a radial distribution), cdf
$P_{\infty}=\int_{-\infty}^{\cdot}p_{\infty}$ and survival function
$\bar{P}_{\infty}=1-P_{\infty}$. To $p_{\infty}$ we associate the “density
Stein operators”

$\displaystyle\mathcal{T}_{\infty}f(x)$
$\displaystyle=\frac{(f(x)p_{\infty}(x))^{\prime}}{p_{\infty}(x)}=f^{\prime}(x)+\frac{p_{\infty}^{\prime}(x)}{p_{\infty}(x)}f(x)$
$\displaystyle\mathcal{L}_{\infty}h(x)$
$\displaystyle=-\int_{-\infty}^{\infty}h^{\prime}(y)\frac{P_{\infty}(y\wedge
x)\bar{P}_{\infty}(y\lor x)}{p_{\infty}(x)}\mathrm{d}y$

acting on absolutely continuous functions $f$ and $h$ such that, moreover,
$f^{\prime}$ and $h$ are integrable with respect to $p_{\infty}$. One can show
that
$\mathcal{T}_{\infty}\mathcal{L}_{\infty}h(x)=h(x)-\mathbb{E}h(F_{\infty})$,
where $\mathcal{L}_{\infty}\mathcal{T}_{\infty}f(x)=f(x)$ for all appropriate
$h$ and $f$. The choice $h(x)=-x$ gives the _Stein kernel_

$\displaystyle\tau_{\infty}(x):=\mathcal{L}_{\infty}h(x)=\int_{-\infty}^{\infty}\frac{P_{\infty}(y\wedge
x)\bar{P}_{\infty}(y\lor x)}{p_{\infty}(x)}\mathrm{d}y.$

It is now well-documented that Stein kernels, if available, provide valuable
insight into the properties of the underlying densities. Following [12, 4], we
propose to apply Stein’s method based on the following Stein operator:

$\displaystyle\mathcal{A}_{\infty}g(x)=\mathcal{T}_{\infty}(\tau_{\infty}(x)g(x))=\tau_{\infty}(x)g^{\prime}(x)-xg(x).$

Note that, by design,
$\mathbb{E}\left[\mathcal{A}_{\infty}g(F_{\infty})\right]=0$ for all $g$ for
which the moments exist. The Stein equation for $p_{\infty}$ associated to
operator $\mathcal{A}_{\infty}$ is then

$\displaystyle\tau_{\infty}(x)g^{\prime}(x)-xg(x)=h(x)-\mathbb{E}h(F_{\infty})$
(14)

and one can easily check that the function

$\displaystyle g_{h}(x)=\frac{\mathcal{L}_{\infty}h(x)}{\tau_{\infty}(x)}$
(15)

solves (14) at all $x$ such that $\tau_{\infty}(x)\neq 0$, so that for all
real random variables $F$ such that ${P}(\tau_{\infty}(F)=0)=0$, it holds that

$\displaystyle\mathbb{E}h(F)-\mathbb{E}h(F_{\infty})$
$\displaystyle=\mathbb{E}\left[\tau_{\infty}(F)g_{h}^{\prime}(F)-Fg_{h}(F)\right]$
(16)

with $g_{h}$ as given in (15). By definition of the Wasserstein distance, we
conclude

$\mathop{d_{\mathrm{W}}}(\mathbb{F},\mathbb{F}_{\infty})\leq\sup_{h\in\mathcal{H}}\left|\mathbb{E}\left[\tau_{\infty}(F)g_{h}^{\prime}(F)-Fg_{h}(F)\right]\right|,$
(17)

where $\mathbb{F}$ is the probability distribution of $F$,
$\mathbb{F}_{\infty}$ is that of $F_{\infty}$ and $\mathcal{H}$ is the family
of Lipschitz functions with constant less than 1 (recall the first lines of
Section 2.4).

Specializing to densities of the form $p_{\infty}(x)=b(x)\varphi(x)$, with
$\varphi(x)$ the standard normal density, leads to the following non-uniform
bounds for $g$; they are based on results from [6] and proved in the
supplementary material (see Lemma 2.4 therein).

######  Proposition 3.1.

Let $p_{\infty}(x)=b(x)\varphi(x)$ be a standardized density on the real line,
with $b(x)$ a positive function such that $b(x)\neq 0$ for all $x\neq 0$.
Define

$R_{\infty}(x)=\frac{1}{{p_{\infty}(x)}}{\int_{-\infty}^{x}{P}_{\infty}(u)\mathrm{d}u\int_{x}^{\infty}\overline{P}_{\infty}(u)\mathrm{d}u}.$

Let $g(=g_{h})$ be as in (15), with $h\in\mathcal{H}$. Then

$\displaystyle|g(x)|\leq 1$ (18) $\displaystyle|g^{\prime}(x)|\leq
2\frac{R_{\infty}(x)}{\tau_{\infty}(x)^{2}}=:\Psi_{1}(x)$ (19)
$\displaystyle|g^{\prime\prime}(x)|\leq\frac{2}{\tau_{\infty}(x)}\left(1+\left|\frac{2x}{\tau_{\infty}(x)}-x+\frac{b^{\prime}(x)}{b(x)}\right|\frac{R_{\infty}(x)}{\tau_{\infty}(x)}\right)=:2\Psi_{2}(x)$
(20)

An important aspect of inequalities (18), (19) and (20), is that they hold
uniformly over all $h\in\mathcal{H}$; our approach to Stein’s method consists
in using this information in combination with (17) in order to bound the
Wasserstein distance.

We now further specialize to the choice $b(x)\propto|x|^{k}$, where $k=d-1$ in
the case of a $d$-dimensional radial distribution. We then have the following
explicit expression for the Stein kernel by straightforward integration of the
expressions involved.

######  Lemma 3.1.

Let $p_{\infty}(x)\propto|x|^{k}\varphi(x)$ for given $k\geq 0$ with
$\varphi(x)$ the standard Gaussian density. Then the Stein kernel

$\displaystyle\tau_{\infty}(x;k)=2^{k/2}e^{x^{2}/2}|x|^{-k}\Gamma(1+k/2,x^{2}/2),$
(21)

where $\Gamma(\alpha,x)=\int_{x}^{\infty}t^{\alpha-1}e^{-t}\mathrm{d}t$ is the
(upper) incomplete gamma function.

######  Example 3.2.

The following results are immediate from (21):

  * •

if $k=0$ (i.e., $p_{\infty}$ is standard Gaussian density) then
$\tau_{\infty}(x)=1$;

  * •

if $k=1$ (i.e., $p_{\infty}$ is two-sided Rayleigh density) then
$\tau_{\infty}(x)=1+\sqrt{\pi/2}|x|^{-1}$
$e^{x^{2}/2}\mathrm{Erfc}(x/\sqrt{2})$ (which behaves like $1/|x|+1$);

  * •

if $k=2$ (i.e., $p_{\infty}$ is two-sided Maxwell density) then
$\tau_{\infty}(x)=2/x^{2}+1$.

It is not hard to obtain expressions for the Stein kernel at any level
$k\in\mathbb{N}$; they are provided in the supplementary material (see Lemma
2.1 therein).

Figure 5: Functions $\Psi_{1}$ (left) and $\Psi_{2}$ (right) defined in (19)
and (20) for $b(x)\propto|x|^{k}$ with $k=0$ (blue), $k=1$ (orange), $k=2$
(green) and $k=21$ (red)

The following properties of the functions defined in (19) and (20) (with
$\tau_{\infty}$ as given in (21)) will be useful as well (see Fig. 5 for
illustration111The supplementary material contains more information on these
functions.): $\Psi_{1}$ and $\Psi_{2}$ are even and positive, and

  * •

if $k=0$ then $\Psi_{1}(x)$ is unimodal with maximum
$\Psi_{1}(0)=\sqrt{2/\pi}$ at $x=0$ and strictly decreasing towards 0 for
$x\geq 0$; $\Psi_{2}(x)$ is unimodal with minimum $\Psi_{2}(0)=1$ at $x=0$ and
strictly increasing towards 2 for $x\geq 0$.

  * •

If $k=1$ then $\Psi_{1}(x)$ is bimodal with minimum $\Psi_{1}(0)=0$ at $x=0$,
maximum value less than 1 after which it is strictly decreasing towards 0 for
$x\geq 0$; $\Psi_{2}(x)$ is unimodal with minimum $\Psi_{2}(0)=1/2$ at $x=0$
and strictly increasing towards 2 for $x\geq 0$.

  * •

If $k\geq 2$ then $\Psi_{1}(x)$ is bimodal with minimum $\Psi_{1}(0)=0$ at
$x=0$, maximum value less than 1 after which it is strictly decreasing towards
0 for $x\geq 0$; $\Psi_{2}(x)$ is unimodal with minimum $\Psi_{2}(0)=0$ at
$x=0$ and strictly increasing towards 2 for $x\geq 0$.

Having presented the key elements of Stein’s method for radial distributions,
we now give some properties of the corresponding approximating sequences. The
following lemma, ensuring a unique symmetric and monotone solution to (4), is
a routine extension of a result in [15, 16] to general $k$. In the sequel we
assume $N$ is even and define the median of an ordered set
$\\{x_{1},\ldots,x_{N}\\}$ to be $(x_{m}+x_{m+1})/2$, where $m=N/2$.

######  Lemma 3.3.

Every zero-median solution $\\{x_{1},\ldots,x_{N}\\}$ of (4) for a given
$k\geq 0$ satisfies:

  *     * (P1)

Zero-mean: $x_{1}+\ldots+x_{N}=0$.

    * (P2)

Variance-bound: $x_{1}^{2}+\ldots+x_{N}^{2}=(k+1)(N-1)$.

    * (P3)

Symmetry: $x_{n}=-x_{N+1-n}$ for $n=1,\ldots,N$.

Further, there exists a unique solution $\\{x_{1},\ldots,x_{N}\\}$ such that
(P1) and

  *     * (P4)

Strictly decreasing: $x_{1}>\ldots>x_{N}$

hold. This solution has the zero-median property, and thus also satisfies (P2)
and (P3).

Let ${\mathbb{F}}_{N}$ be the empirical distribution of the unique solution
$\left\\{x_{1},\ldots,x_{N}\right\\}$. For the case $k=0$, the following
result gives the optimal rate of convergence in Wasserstein distance of
${\mathbb{F}}_{N}$ to standard Gaussian. This result is not new, as the upper
bound follows by adapting the coupling argument in [16], and the lower bound,
which follows from a very delicate analysis of the $x_{i}$’s, was recently
proved in [3]. We provide a new proof for the upper bound in Example 3.7
below.

######  Proposition 3.2.

The empirical distribution $\mathbb{F}_{N}$ on the unique symmetric and
monotone solution to (1) (or (4) with $k=0$) satisfies

$\displaystyle\mathop{d_{\mathrm{W}}}({\mathbb{F}}_{N},{\mathbb{F}}_{\infty})\asymp\sqrt{\log
N}/N,$

where ${\mathbb{F}}_{\infty}$ is the standard Gaussian distribution.

Upper bounds on the rate of convergence for values of $k\geq 1$ can also be
obtained by adapting the arguments in [16] using the coupling approach (as
developed in detail in the Supplementary material). The following result for
the case $k=1$ is based on this approach (and is detailed in Section 2.3 of
the Supplementary material), but it is not expected to yield the optimal rate
which, as argued in the sequel, we expect to be the same as the Gaussian case
in the above proposition.

######  Proposition 3.3.

The empirical distribution $\mathbb{F}_{N}$ on the unique symmetric and
monotone solution to (5) (or (4) with $k=1$) satisfies

$\displaystyle\mathop{d_{\mathrm{W}}}({\mathbb{F}}_{N},{\mathbb{F}}_{\infty})=O({\log
N}/N),$

where ${\mathbb{F}}_{\infty}$ is the two-sided Rayleigh distribution.

The main result of this section is the following theorem based on a discrete
density approach to Stein’s method, cf. the approach of [8] in the case of
points on an integer grid. This result is entirely new and may be of
independent interest as the upper bound it provides does not rely on any
specific structure of $x_{1},\ldots,x_{N}$ apart from symmetry and
monotonicity, and therefore may be much more widely applicable. It leads to
the new proof of the tight upper bound in Proposition 3.2 for the case $k=0$,
as shown in Example 3.7 below. For $k\geq 1$, it becomes more difficult
analytically to determine the order of the bound, but it nevertheless provides
a tight bound that can be evaluated numerically and that agrees with the
asymptotic order of
$\mathop{d_{\mathrm{W}}}({\mathbb{F}}_{N},{\mathbb{F}}_{\infty})$ as
$N\to\infty$, as discussed in Example 3.8 below.

######  Theorem 3.4.

Let $N$ be an even integer and $x_{1},\ldots,x_{N}$ be any symmetric and
strictly monotone decreasing sequence. Set
$\tau_{N}(x_{i})=(x_{i}-x_{i+1})\sum_{j=1}^{i}x_{j}$ for $1\leq i\leq N-1$,
and $\tau_{N}(x_{N})=0$. Let ${\mathbb{F}}_{N}$ be the empirical distribution
of these points and let ${\mathbb{F}}_{\infty}$ be the probability measure
with density $p_{\infty}$ satisfying the assumptions of Proposition 3.1. Then

$\displaystyle\mathop{d_{\mathrm{W}}}({\mathbb{F}}_{N},{\mathbb{F}}_{\infty})$
$\displaystyle\leq\frac{1}{N}\sum_{i=1}^{N}\left|\tau_{\infty}(x_{i}){-\tau_{N}(x_{i})}\right|\Psi_{1}(x_{i})$
(22)
$\displaystyle\quad+\frac{1}{N}\sum_{i=1}^{N-1}{|x_{i}-x_{i+1}|\tau_{N}(x_{i})}\max\\{\Psi_{2}(x_{i}),\Psi_{2}(x_{i+1})\\},$
(23)

where $\Psi_{1}$ and $\Psi_{2}$ are given in (19) and (20), respectively.

######  Remark 3.5.

Equation (22) encourages us to consider the quantity
$\tau_{N}(x_{i})=(x_{i}-x_{i+1})\sum_{j=1}^{i}x_{j}$ as a Stein kernel for the
uniform distribution on the set $\left\\{x_{1},\ldots,x_{N}\right\\}$. This
function is non-negative and, under the conditions that we are considering, it
holds that $\tau_{N}(x)\approx\tau_{\infty}(x)$ throughout the support of
${\mathbb{F}}_{N}$ except close to the edges and the origin (see Figs. 6 and
7).

Figure 6:  Superimposition of $\tau_{N}(x_{i}),1\leq i\leq N$ (orange bars)
and $\tau_{\infty}(x_{i}),1\leq i\leq N$ (blue bars) when
$\tau_{\infty}(x)=1+2/x^{2}$ is the Stein kernel for the two-sided Maxwell
distribution and $(x_{i})_{1\leq i\leq N}$ is as in Lemma 3.3 with $k=2$ and
$N=20$ (left), $N=110$ (middle) and $N=300$ (right)

Figure 7:  Plot of $|\tau_{N}(x_{i})-\tau_{\infty}(x_{i})|,1\leq i\leq N$ when
$\tau_{\infty}(x)=1+2/x^{2}$ is the Stein kernel for the two-sided Maxwell
distribution and $(x_{i})_{1\leq i\leq N}$ is as in Lemma 3.3 with $k=2$ and
$N=20$ (left), $N=110$ (middle) and $N=300$ (right)

######  Remark 3.6.

Since both $\Psi_{1}$ and $\Psi_{2}$ are uniformly bounded by 2, we
immediately obtain the simpler upper bound

$\displaystyle\mathop{d_{\mathrm{W}}}({\mathbb{F}}_{N},{\mathbb{F}}_{\infty})$
$\displaystyle\leq\frac{1}{N}\sum_{i=1}^{N}\left|\tau_{\infty}(x_{i})-\tau_{N}(x_{i})\right|+\frac{1}{N}\sum_{i=1}^{N-1}|x_{i}-x_{i+1}|\tau_{N}(x_{i}).$

This simplicity comes at the cost of a loss of precision in the rate. For
instance, in the case that we consider here, the presence of the functions
$\Psi_{1}$ and $\Psi_{2}$ in the bounds controls the problems around the
origin.

######  Example 3.7 (Case $k=0$).

The sequence $\\{x_{n}\\}$ satisfies
$x_{n}-x_{n+1}=({\sum_{i=1}^{n}x_{i}})^{-1}$ so that $\tau_{N}(x_{i})=1$ for
all $i=1,\ldots,N-1$. Since $\tau_{\infty}(x_{i})=1$, (22) vanishes except for
the term corresponding to $i=N$ leading to $\Psi_{1}(x_{N})/N\leq 1/N$. Also,
(23) reads

$\displaystyle\frac{1}{N}\sum_{i=1}^{N-1}|x_{i+1}-x_{i}|\max\\{\Psi_{2}(x_{i}),\Psi_{2}(x_{i+1})\\}.$

Using $\Psi_{2}\leq 2$ along with the symmetry of the sequence, we then have

$\displaystyle\mathop{d_{\mathrm{W}}}({\mathbb{F}}_{N},{\mathbb{F}}_{\infty})\leq\frac{1}{N}+\frac{2}{N}\sum_{i=1}^{N-1}|x_{i+1}-x_{i}|=\frac{1+4x_{1}}{N}=O(\sqrt{\log
N}/N),$

where the last step follows by a version of the argument in [16, Lemma 4.8]
showing that $x_{1}=O(\sqrt{\log N})$. This proves the upper bound part of
Proposition 3.2.

######  Example 3.8 (Cases $k\geq 1$).

Evaluating the bound in Theorem 3.4 based on numerical solutions to (4) and
comparing the results with values of
$\mathop{d_{\mathrm{W}}}(\mathbb{F}_{N},{\mathbb{F}}_{\infty})$ calculated
using numerical integration, over a range of values of $N$, suggests that the
bound is of the same order as
$\mathop{d_{\mathrm{W}}}(\mathbb{F}_{N},{\mathbb{F}}_{\infty})$, and moreover
that

$\displaystyle\mathop{d_{\mathrm{W}}}({\mathbb{F}}_{N},{\mathbb{F}}_{\infty})\asymp\sqrt{\log
N}/N\mbox{ for }k=1,2,$
$\displaystyle\mathop{d_{\mathrm{W}}}({\mathbb{F}}_{N},{\mathbb{F}}_{\infty})\asymp({\log
N})^{6}/N^{2}\mbox{ for }k\geq 3.$

To obtain these rates we used Mathematica to compute the bounds and R to
compute the Wasserstein distance using its representation as the
$L_{1}$-distance between the cdfs (the programs are provided in the
Supplementary material).

######  Remark 3.9.

We emphasize that the bound in Theorem 3.4 holds irrespective of the sequence
that is chosen; it may be of interest to optimize the approximation of
$p_{\infty}$ by a sample $x_{1}<\ldots<x_{N}$ which minimizes the right hand
side. One simple way to do this is to require that the sequence satisfies the
recurrence

$x_{i+1}=x_{i}-\frac{\tau_{\infty}(x_{i})}{\sum_{j=1}^{i}x_{j}},$

thereby canceling out (22) and only leaving (23); note that in the Gaussian
case $\tau_{\infty}=1$ so we are back with the recurrence (1).

###### Proof of Theorem 3.4.

The idea is to apply a version of Stein’s density method to
$F_{N}\sim{\mathbb{F}}_{N}$. Note that a discrete “derivative” at $x_{i}$,
$i=1,\ldots,N-1$, is given by

$D_{N}f(x_{i})=f(x_{i})-f(x_{i+1}).$

Writing $\sigma(x_{i})=\sum_{j=1}^{i}x_{j}$, straightforward summation (along
with the identity $x_{N}=-\sum_{j=1}^{N-1}x_{j}$) shows

$\mathbb{E}\left[\sigma(F_{N})D_{N}g(F_{N})-F_{N}g(F_{N})\right]=0$

for all summable functions $g$. Subtracting this from (16) applied to
$F=F_{N}$ gives

$\displaystyle\mathbb{E}h(F_{N})-\mathbb{E}h(F_{\infty})$
$\displaystyle=\mathbb{E}\left[\tau_{\infty}(F_{N})g^{\prime}(F_{N})-F_{N}g(F_{N})\right]$
$\displaystyle=\mathbb{E}\left[\tau_{\infty}(F_{N})g^{\prime}(F_{N})-\sigma(F_{N})D_{N}g(F_{N})\right].$

From (17) it then follows that

$\mathop{d_{\mathrm{W}}}(\mathbb{F}_{N},\mathbb{F}_{\infty})\leq\sup_{h\in\mathcal{H}}\left|\mathbb{E}\left[\tau_{\infty}(F_{N})g^{\prime}(F_{N})-\sigma(F_{N})D_{N}g(F_{N})\right]\right|$

with $g(=g_{h})$ satisfying the bounds in Proposition 3.1. By Taylor’s
theorem,

$\displaystyle D_{N}f(x_{i})={(x_{i}-x_{i+1})f^{\prime}(x_{i})-}{\small{1\over
2}}{(x_{i}-x_{i+1})^{2}}f^{\prime\prime}(c_{i}){\mbox{ for all }1\leq i\leq
N-1,}$

where $c_{i}$ is between $x_{i+1}$ and $x_{i}$. Hence

$\displaystyle\mathbb{E}h(F_{N})-\mathbb{E}h(F_{\infty})$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{{N-1}}\left(\tau_{\infty}(x_{i})-{\sigma(x_{i})}{(x_{i}-x_{i+1})}\right)g^{\prime}(x_{i}){{-}\frac{1}{2N}\sum_{i=1}^{{N-1}}{\sigma(x_{i})}({x_{i}-x_{i+1}})^{2}g^{\prime\prime}(c_{i})}$
$\displaystyle\qquad+\tau_{\infty}(x_{N})g^{\prime}(x_{N})$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{{N}}\left(\tau_{\infty}(x_{i})-\tau_{N}(x_{i})\right)g^{\prime}(x_{i})-\frac{1}{2N}\sum_{i=1}^{{N-1}}(x_{i+1}-x_{i})\tau_{N}(x_{i})g^{\prime\prime}(c_{i})$

(recall that $\tau_{N}(x_{N})=0$). Taking absolute values, using (19) and
(20), along with the symmetry/unimodality of the function $\Psi_{2}$, we get
the result. ∎

##  4 Convergence of multidimensional ground states

In this section we apply the results of the previous section to the $N$-point
empirical distributions $\mathbb{P}_{N}$ of the $d$-dimensional ground states
discussed in Sections 2.1 and 2.2.

These states have (signed) radial components given by the unique zero-median
solution $x_{1}>\ldots>x_{R}$ to the recursion (4) with $b(x)=|x|^{k}$, where
$k=d-1\geq 1$ and $R=N^{1/2-1/(2d)}$ plays the role of $N$ in (4),
corresponding to $R$ points in each of $N^{1/2+1/(2d)}$ radial directions. By
translating the results in Example 3.8 into the $d$-dimensional setting, the
optimal convergence rates of the radial part of $\mathbb{P}_{N}$ to the radial
part of the $d$-dimensional normal distribution ${\cal N}_{d}$ are seen to be
given by

$d=2,3$: $\mathop{d_{\mathrm{W}}}(\mathbb{P}_{N}^{\rm radial},{\cal
N}_{d}^{\rm radial})\asymp\sqrt{\log R}/R$;

$d\geq 4$: $\mathop{d_{\mathrm{W}}}(\mathbb{P}_{N}^{\rm radial},{\cal
N}_{d}^{\rm radial})\asymp({\log R})^{6}/R^{2}$.

The number of distinct polar angles in the directional part of the ground
state (each with $d-2$ components taking values in $[0,\pi/2]$ when $d\geq 3$,
and a single component in $[0,\pi)$ when $d=2$) is of order $N^{1/d}$. For
$d\geq 3$, a larger number (namely $N^{1/2-1/(2d)}$) of azimuthal angles
corresponding to each polar angle are available, so the polar part has the
slower convergence rate. In general, the Wasserstein distance between the
directional part of $\mathbb{P}_{N}$ and that of ${\cal N}_{d}$ is thus of
order $N^{-1/d}$ (this follows easily using the representation of Wasserstein
distance as the $L_{1}$-distance between two cdfs). For $d\geq 4$, in view of
the upper bound in Lemma 2.1, namely a sum involving separate contributions
from the the radial part and the directional parts, we conclude that the
overall rate that $\mathbb{P}_{N}$ tends to ${\cal N}_{d}$ is of order
$N^{-1/d}$. For $d=3$, the radial component involving $N^{1/3}$ points
dominates, however, so the overall rate is of order $\sqrt{\log N}/N^{1/3}$.
Similarly, the radial component involving $N^{1/4}$ points dominates in the
case $d=2$, so the overall rate is of order $\sqrt{\log N}/N^{1/4}$.

In summary, we obtain

$d=2$: $\mathop{d_{\mathrm{W}}}(\mathbb{P}_{N},{\cal N}_{d})\asymp\sqrt{\log
N}/N^{1/4}$;

$d=3$: $\mathop{d_{\mathrm{W}}}(\mathbb{P}_{N},{\cal N}_{d})\asymp\sqrt{\log
N}/N^{1/3}$;

$d\geq 4$: $\mathop{d_{\mathrm{W}}}(\mathbb{P}_{N},{\cal N}_{d})\asymp
N^{-1/d}$.

Interestingly, the fastest rate of convergence to the ground state occurs in
three dimensions.

For the excited states discussed in Section 2.3, the rates of convergence are
of order $N^{-1/d}$; in this case the directional components dominate because
the radial component has the faster rate (as seen by applying Example 3.8 with
$k=d+1$ for $d=2,3$).

## References

  * [1] J. S. Brauchart and P. J. Grabner.  Distributing many points on spheres: Minimal energy and designs.  Journal of Complexity, 31(3):293 – 326, 2015.
  * [2] L. H. Y. Chen, L. Goldstein, and Q.-M. Shao.  Normal approximation by Stein’s method.  Probability and its Applications (New York). Springer, Heidelberg, 2011\.
  * [3] L. H. Chen and L. V. Thành.  Optimal bounds in normal approximation for many interacting worlds.  arXiv:2006.11027, 2020.
  * [4] C. Döbler.  Stein’s method of exchangeable pairs for the Beta distribution and generalizations.  Electronic Journal of Probability, 20(109):1–34, 2015.
  * [5] M. Ernst, G. Reinert, and Y. Swan.  First-order covariance inequalities via Stein’s method.  Bernoulli, 26(3):2051–2081, 2020.
  * [6] M. Ernst and Y. Swan.  Distances between distributions via Stein’s method.  Journal of Theoretical Probability, 2020.
  * [7] M. Ghadimi, M. J. W. Hall, and H. M. Wiseman.  Nonlocality in Bell’s theorem, in Bohm’s theory, and in many interacting worlds theorising.  Entropy, 20(8), 2018.
  * [8] L. Goldstein and G. Reinert.  Stein’s method for the Beta distribution and the Pólya-Eggenberger urn.  Journal of Applied Probability, 50(4):1187–1205, 2013.
  * [9] M. J. W. Hall, D.-A. Deckert, and H. M. Wiseman.  Quantum phenomena modeled by interactions between many classical worlds.  Physical Review X, 4:041013, Oct 2014.
  * [10] H. Herrmann, M. J. W. Hall, H. M. Wiseman, and D. A. Deckert.  Ground states in the many interacting worlds approach.  arXiv:1712.01918, 2018.
  * [11] G. Jameson.  The incomplete gamma functions.  The Mathematical Gazette, 100(548):298–306, 2016.
  * [12] C. Ley, G. Reinert, and Y. Swan.  Distances between nested densities and a measure of the impact of the prior in Bayesian statistics.  The Annals of Applied Probability, 27(1):216–241, 2017.
  * [13] S. Li.  Concise formulas for the area and volume of a hyperspherical cap.  Asian Journal of Mathematics and Statistics, 4:66–70, 2011.
  * [14] E. Mariucci and M. Reiß.  Wasserstein and total variation distance between marginals of Lévy processes.  Electronic Journal of Statistics, 12(2):2482–2514, 2018.
  * [15] I. W. McKeague and B. Levin.  Convergence of empirical distributions in an interpretation of quantum mechanics.  The Annals of Applied Probability, 26(4):2540–2555, 2016.
  * [16] I. W. McKeague, E. A. Peköz, and Y. Swan.  Stein’s method and approximating the quantum harmonic oscillator.  Bernoulli, 25(1):89–111, 02 2019.

##  Appendix A Supplementary material

All code for performing computations is available at the url

https://tinyurl.com/mcks2021code

The code is also available at the end of the supplementary material.

###  A.1 Preliminary remarks on Stein operators

Consider a probability distribution $\mathbb{P}$ with cdf $P$, pdf $p$ w.r.t.
the Lebesgue measure on $\mathbb{R}$. Suppose that $p$ itself is absolutely
continuous, with a.e. derivative $p^{\prime}(x)$. We denote $L^{1}(p)$ the
collection of functions such that
$\int_{-\infty}^{\infty}|h(x)|p(x)\mathrm{d}x<\infty$ and write
$P(h)=\mathbb{E}_{p}h=\int_{-\infty}^{\infty}h(x)p(x)\mathrm{d}x.$ We also
denote $\mathcal{F}^{(0)}(p)$ the collection of all mean 0 functions under
$p$. Following [5], to $p$ we associate the Stein operators

$\displaystyle\mathcal{T}_{p}f(x)=\frac{(f(x)p(x))^{\prime}}{p(x)}$ (24)
$\displaystyle\mathcal{L}_{p}h(x)=\frac{1}{p(x)}\int_{-\infty}^{x}(h(u)-P(h))p(u)\mathrm{d}u$
(25)

with the convention that $\mathcal{T}_{p}f(x)=\mathcal{L}_{p}h(x)=0$ for all
$x$ such that $p(x)=0$. In the sequel we denote
$\mathcal{S}(p)=\left\\{x\,|\,p(x)>0\right\\}$, $a=\inf\mathcal{S}(p)$ and
$b=\sup\mathcal{S}(p)$; we assume that $\mathcal{S}(p)$ is the union of a
finite number of intervals.

Of course (24) is only defined for functions $f$ such that $fp$ is absolutely
continuous. We denote $\mathcal{F}(p)$ the collection of functions $f$ such
that not only is $fp$ absolutely continuous, but also $(f(x)p(x))^{\prime}\in
L^{1}(p)$ and $\lim_{x\to a}f(x)p(x)=\lim_{x\to b}f(x)p(x)=0$. This class of
functions is important because $\mathcal{T}_{p}f\in\mathcal{F}^{(0)}(p)$ for
all $f\in\mathcal{F}(p)$; this is crucial for Stein’s method as it gives rise
to many “Stein identities” which can be used for a variety of purposes.
Similarly, (25) is only defined for functions $h\in L^{1}(p)$, in which case
$\mathcal{L}_{p}h\in\mathcal{F}(p)$ for all $h\in L^{1}(p)$.

As described in [12], it is interesting to “standardize” the operator (24) by
fixing some $c\in\mathcal{F}(p)$ and considering the family of “standardized
Stein operators”
$\mathcal{A}_{c}f(x)=\mathcal{T}_{p}(cf)(x)=c(x)f^{\prime}(x)+\Big{(}c^{\prime}(x)+{p^{\prime}(x)}/{p(x)}c(x)\Big{)}f(x)$
acting on some class $\mathcal{F}(\mathcal{A}_{c})$ made of all functions for
which $cf\in\mathcal{F}(p)$. Note that, by definition, we have
$\mathbb{E}[\mathcal{A}_{c}f(X)]=0$ for all
$f\in\mathcal{F}(\mathcal{A}_{c})$. It is important that $c$ be well-chosen to
ensure that $\mathcal{A}_{c}$ has a manageable expression; as is now well
known, there are many instances of densities $p$ (even intractable densities)
for which this turns out to be possible, leading to many powerful handles on
$p$ which can then serve for a variety of purposes including but not limited
to distributional approximation.

Given an operator $\mathcal{A}_{c}$, classical instantiations of Stein’s
method begin with a “Stein equation”, i.e. a differential equation of the form

$\mathcal{A}_{c}f(x)=h(x)-P(h)$ (26)

for $h$ some function belonging to a class $\mathcal{H}$ of test functions.
Typically, Stein’s method practitionners work with one of the following
classes: (i)
$h\in\mathcal{H}:=\mathrm{Kol}(\mathbb{R})=\left\\{\mathbb{I}(-\infty,z],z\in\mathbb{R}\right\\}$
the indicators of a lower half line; (ii)
$h\in\mathcal{H}:=\mathrm{TV}(\mathbb{R})$ the collection of functions such
that $\|h\|\leq 1$; (iii) $h\in\mathcal{H}:=\mathrm{Wass}(\mathbb{R})$ the
collection of Lipschitz functions such that $\|h^{\prime}\|\leq 1$. In the
sequel, we restrict our attention to $\mathcal{H}=\mathrm{Wass}(\mathbb{R})$,
and we assume that for each $h\in\mathcal{H}$ there exists a unique function
$f\in\mathcal{F}(\mathcal{A}_{c})$ for which (26) holds for all
$x\in\mathcal{S}(p)$. Under “reasonable assumptions on $p$” (to be verified on
a case-by-case basis) we can write
$\mathcal{T}_{p}\mathcal{L}_{p}h=h-\mathbb{E}_{p}h$ for all $h\in L^{1}(p)$
and in particular $\mathcal{T}_{p}\mathcal{L}_{p}=\mathrm{Id}$ over
$\mathcal{F}^{(0)}(p)$ ($\mathrm{Id}$ is the identity function). Similarly
$\mathcal{L}_{p}\mathcal{T}_{p}=\mathrm{Id}$ over $\mathcal{F}(p)$. In other
words, under “reasonable assumptions on $p$”, the solution to (26) is
$f_{h}(x)={\mathcal{L}_{p}h(x)}/{c(x)}$ at all $x\in\mathcal{S}(p)$ for which
$c(x)\neq 0$. Then, since the Wasserstein distance between two probability
measures $\mathbb{P}$ and $\mathbb{Q}$ can be written as
$\mathop{d_{\mathrm{W}}}(\mathbb{P},\mathbb{Q})=\sup_{h\in\mathrm{Wass}(\mathbb{R})}|\mathbb{E}h(X)-\mathbb{E}h(Y)|$
where $X\sim\mathbb{P}$ and $Y\sim\mathbb{Q}$, it holds that

$\mathop{d_{\mathrm{W}}}(\mathbb{P},\mathbb{Q})=\sup_{h\in\mathrm{Wass}(1)}\left|\mathbb{E}\left[\mathcal{A}_{c}f_{h}(Y)\right]\right|;$
(27)

Stein’s method in Wasserstein distance consists in exploiting this last
identity for the purpose of estimating the Wasserstein distance between the
laws $\mathbb{P}$ and $\mathbb{Q}$.

In order to be able to use (27) successfully, it is crucial to control
solutions $f_{h}$ and their derivatives. In [6] the following representations
for (25) are provided (recall that $h$ is Lipschitz with a.e. derivative
$h^{\prime}$):

$\displaystyle-\mathcal{L}_{p}h(x)$
$\displaystyle=-\mathbb{E}\left[({h}(X)-P(h))\frac{\mathbb{I}[X\leq
x]}{p(x)}\right]$
$\displaystyle=-\mathbb{E}\left[({h}(X)-P(h))\frac{\mathbb{I}[X\leq
x]-P(x)}{p(x)}\right]$
$\displaystyle=\mathbb{E}\left[(h(X_{2})-h(X_{1}))\frac{\mathbb{I}[X_{1}\leq
x\leq X_{2}]}{p(x)}\right]$ (28)
$\displaystyle=\mathbb{E}\left[h^{\prime}(X)\frac{P(x\wedge X)(1-P(x\lor
X))}{p(x)p(X)}\right]$ (29)

where, in (28), the random variables $X_{1},X_{2}$ are independent copies of
$X$. A simpler way to write (29) is

$\displaystyle-\mathcal{L}_{p}h(x)=\int_{-\infty}^{\infty}h^{\prime}(y)\frac{P_{\infty}(y\wedge
x)\bar{P}_{\infty}(y\lor x)}{p_{\infty}(x)}\mathrm{d}y$

It is also shown that

$\displaystyle\bar{h}(x):=h(x)-P(h)=\mathbb{E}\left[h^{\prime}(X)\frac{P(X)-\mathbb{I}[x\leq
X]}{p(X)}\right]$

for all $x\in\mathcal{S}(p)$. All these representations will be used in the
next section to control the solutions to the Stein equations; this in turn
will lead to the distributional approximation results.

###  A.2 Stein’s method for radial distributions

####  A.2.1 Notations and background

Before specializing to radial densities, it is enlightening to first widen the
scope somewhat and consider targets $F_{\infty}$ with density of the form
$p_{\infty}(x)=b(x)\gamma(x)$, for $\gamma$ some “basis density” and $b$ some
positive $\gamma$-integrable “tilting” function. This theory may also be of
independent interest.

First note that, in order for $p_{\infty}$ to be a density, it is necessary
that $b\geq 0\in L^{1}(\gamma)$ and $\mathbb{E}[b(Z)]=1$, where here and
throughout we denote $Z\sim\gamma$. We further impose the following
assumptions on $p_{\infty}$. First, we require that $\gamma$ is a
differentiable probability density function with support the full real line,
such that moreover $\gamma^{\prime}\in L^{1}(\mathrm{d}x)$ has exactly one
sign change (which, for simplicity, we fix at 0) and
$\int\gamma^{\prime}(x)\mathrm{d}x=0$. Second, we let $B(x)$ be an absolutely
continuous nondecreasing function with continuous derivative $b$, we denote
$\mathcal{S}_{B}=\left\\{x\in\mathbb{R}\right.$, such that
$\left.b(x)>0\right\\}$ and suppose that $\mathcal{S}_{B}$ is the union of a
finite number of intervals. Following [12], we also introduce
$\mathcal{F}(\gamma)$ the Stein class of $\gamma$; this is the class of
functions $f:\mathbb{R}\to\mathbb{R}$ such that $(f\gamma)^{\prime}\in
L^{1}(\mathrm{d}x)$ and
$\int_{\mathbb{R}}(f(x)\gamma(x))^{\prime}\mathrm{d}x=0$. We assume that
$B\in\mathcal{F}(\gamma)$; since $b\in L^{1}(\gamma)$, this assures us that
$\int b\gamma=-\int B\gamma^{\prime}$ so that integration by parts holds
without a remainder term. Finally, letting $F_{\infty}\sim p_{\infty}$, we
impose that $\mathbb{E}F_{\infty}(=\mathbb{E}[Zb(Z)])=0$.

With these assumptions we are now ready to provide a Stein’s method theory for
$p_{\infty}=b\gamma$; the backbone of our approach comes from [16].

######  Definition A.1 (Generalized $(b,\gamma)$-bias transformation).

Suppose that $F$ is such that $P(F\in\mathcal{S}_{{B}})=1$ and define

$\displaystyle\sigma^{2}_{B}(F)=\mathbb{E}\left[-\frac{\gamma^{\prime}(F)}{\gamma(F)}\frac{B(F)}{b(F)}\right].$

The random variable $F^{\star}$ satisfying

$\displaystyle\sigma^{2}_{B}(F)\mathbb{E}\left[\frac{f^{\prime}(F^{\star})}{b(F^{\star})}\right]=\mathbb{E}\left[-\frac{\gamma^{\prime}(F)}{\gamma(F)}\frac{f(F)}{b(F)}\right]$

for all $f$ such that both integrals exist is said to have the generalized
$(b,\gamma)$-bias distribution. The random variable $F^{\star}$ is the
generalized $(b,\gamma)$-bias transform of $F$.

By construction, we always have

$\displaystyle\sigma^{2}_{B}(F_{\infty})$
$\displaystyle=\mathbb{E}\left[-\frac{\gamma^{\prime}(F_{\infty})}{\gamma(F_{\infty})}\frac{B(F_{\infty})}{b(F_{\infty})}\right]=-\int_{\mathcal{S}_{B}}\gamma^{\prime}(x)B(x)\mathrm{d}x=\int_{\mathcal{S}_{B}}b(x)\gamma(x)\mathrm{d}x=\mathbb{E}b(Z)=1.$

Moreover, for any sufficiently regular function $f$:

$\displaystyle\mathbb{E}\left[-\frac{\gamma^{\prime}(F_{\infty})}{\gamma(F_{\infty})}\frac{f(F_{\infty})}{b(F_{\infty})}\right]=\mathbb{E}\left[\frac{f^{\prime}(F_{\infty})}{b(F_{\infty})}\right].$

Therefore $F_{\infty}=F^{\star}$, i.e. $p_{\infty}$ is a fixed point of the
generalized $(b,\gamma)$-bias transform. More generally, the following holds
true.

######  Lemma A.2.

If $\sigma^{2}_{B}(F)<\infty$ then its generalized $(b,\gamma)$-bias transform
$F^{\star}$ exits and is absolutely continuous with density

$\displaystyle
p^{\star}(x)=-\frac{b(x)}{\sigma^{2}_{B}(F)}\mathbb{E}\left[\frac{1}{b(F)}\frac{\gamma^{\prime}(F)}{\gamma(F)}\mathbb{I}[F\geq
x]\right].$

Moreover $F_{\infty}$ is the unique fixed point of this tranformation, in the
sense that if ${F}\stackrel{{\scriptstyle\mathcal{D}}}{{=}}F^{\star}$ then
$F\stackrel{{\scriptstyle\mathcal{D}}}{{=}}F_{\infty}$ (equality in
distribution).

###### Proof.

The proof of unicity is easy; the other points follow from arguments nearly
identical to those in [2, Proposition 2.1]. ∎

Now consider the function

$\displaystyle
f(x)=\frac{1}{\gamma(x)}\int_{-\infty}^{x}\big{(}h(u)-\mathbb{E}h(F_{\infty})\big{)}b(u)\gamma(u)\mathrm{d}u=:\frac{1}{\gamma(x)}\int_{-\infty}^{x}\bar{h}(u)b(u)\gamma(u)\mathrm{d}u$

which is solution to the differential equation

$\displaystyle(\bar{h}(x):=)h(x)-\mathbb{E}h(F_{\infty})=\frac{f^{\prime}(x)+\gamma^{\prime}(x)/\gamma(x)f(x)}{b(x)}$

for all $x\in\mathcal{S}_{B}$. Let $F$ be a random variable such that
$P(F\in\mathcal{S}_{B})=1$ and $\sigma^{2}_{B}(F)=1$. We have

$\displaystyle\mathbb{E}h(F)-\mathbb{E}h(F_{\infty})=\mathbb{E}\left[\frac{f^{\prime}(F)}{b(F)}+\frac{\gamma^{\prime}(F)}{\gamma(F)}\frac{f(F)}{b(F)}\right]=\mathbb{E}\left[\frac{f^{\prime}(F)}{b(F)}-\frac{f^{\prime}(F^{\star})}{b(F^{\star})}\right]$
(30)

and it remains to express the right hand side of (30) in terms of manageable
quantities, such as moments of $F,F^{\star}$ and $F-F^{\star}$. We cannot work
directly with the function $x\mapsto f^{\prime}(x)/b(x)$ because the latter is
unbounded at $x=0$. To bypass this difficulty, we introduce the notation

$\mathcal{L}_{\infty}h(x)=\frac{1}{p_{\infty}(x)}\int_{-\infty}^{x}\bar{h}(u)p_{\infty}(u)\mathrm{d}u$

(we stress that $\mathcal{L}_{\infty}\neq\mathcal{L}_{\gamma})$ and follow
[16] by introducing the function $g=g_{\eta,h}$ given by

$\displaystyle
g(x)=\frac{\mathcal{L}_{\infty}h(x)}{\mathcal{L}_{\infty}\eta(x)}=\frac{\int_{-\infty}^{x}(h(u)-\mathbb{E}h(F_{\infty}))b(u)\gamma(u)\mathrm{d}u}{\int_{-\infty}^{x}(\eta(u)-\mathbb{E}\eta(F_{\infty}))b(u)\gamma(u)\mathrm{d}u}$
(31)

at all $x\in\mathcal{S}_{B}$ where $h$ is fixed by the left hand side of (30)
but $\eta$ is kept unspecified, to be tuned to our needs at a later stage.
Obviously, the above relations are only defined at $x$ such that
$p_{\infty}(x)\neq 0$; we suppose this to be the case here and in the sequel.
The function $g$ from (31) is then solution to the Stein equation

$\displaystyle\left(\mathcal{L}_{\infty}\eta(x)\right)g^{\prime}(x)+\bar{\eta}(x)g(x)=h(x)-\mathbb{E}h(F_{\infty})$

at all $x$ inside the support of $p_{\infty}$. It will be useful to note that
the functions $g,h$ and $\eta$ satisfy the relations

$\displaystyle(\mathcal{L}_{\infty}\eta)g$
$\displaystyle=\mathcal{L}_{\infty}h$
$\displaystyle(\mathcal{L}_{\infty}\eta)g^{\prime}$
$\displaystyle=\overline{h}-\overline{\eta}g=\overline{h}-\overline{\eta}\frac{\mathcal{L}_{\infty}h}{\mathcal{L}_{\infty}\eta}$
$\displaystyle(\mathcal{L}_{\infty}\eta)g^{\prime\prime}$
$\displaystyle=h^{\prime}-(\overline{\eta}+(\mathcal{L}_{\infty}\eta)^{\prime})g^{\prime}-\eta^{\prime}g$
(32)
$\displaystyle=\left(\overline{\eta}\frac{\overline{\eta}+\left(\mathcal{L}_{\infty}\eta\right)^{\prime}}{\mathcal{L}_{\infty}\eta}-\eta^{\prime}\right)g-\left(\overline{h}\frac{\overline{\eta}+\left(\mathcal{L}_{\infty}\eta\right)^{\prime}}{\mathcal{L}_{\infty}\eta}-h^{\prime}\right).$

Straightforward manipulations of the definitions also lead to

$\displaystyle\mathcal{L}_{\infty}h(x)=\frac{\mathcal{L}_{\gamma}(hb)(x)}{b(x)}-\mathbb{E}[h(Z)b(Z)]\frac{\mathcal{L}_{\gamma}b(x)}{b(x)}.$

Finally we note that $g$ and $f$ are related through
$f(x)=\mathcal{L}_{\infty}\eta(x)g(x)b(x)$ so that

$\displaystyle\frac{f^{\prime}(x)}{b(x)}=\mathcal{L}_{\infty}\eta(x)g^{\prime}(x)+\left(\overline{\eta}(x)-\frac{\gamma^{\prime}(x)}{\gamma(x)}\mathcal{L}_{\infty}\eta(x)\right)g(x)=:\mathcal{L}_{\infty}\eta(x)g^{\prime}(x)+\Psi_{\infty}\eta(x)g(x).$

Identity (30) becomes

$\displaystyle\mathbb{E}h(F)-\mathbb{E}h(F_{\infty})$
$\displaystyle=\mathbb{E}\Big{[}\mathcal{L}_{\infty}\eta(F)g^{\prime}(F)-\mathcal{L}_{\infty}\eta(F^{\star})g^{\prime}(F^{\star})\Big{]}$
$\displaystyle\quad+\mathbb{E}\Big{[}\Psi_{\infty}\eta(F)g(F)-\Psi_{\infty}\eta(F^{\star})g(F^{\star})\Big{]}$
(33)

which is close to what is required. This is however not exactly what we need
because, although we shall see that for reasonable choices of $\eta$, the
function $g$ from (31) and its derivative $g^{\prime}$ are bounded, the second
derivative $g^{\prime\prime}$ is often not. In order to cater for this, we
introduce some further degrees of liberty in the expressions and rewrite (33)
as

$\displaystyle\mathbb{E}h(F)-\mathbb{E}h(F_{\infty})=\mathbb{E}\Big{[}\Big{(}r_{1}(F)\mathcal{L}_{\infty}\eta(F)-r_{1}(F^{\star})\mathcal{L}_{\infty}\eta(F^{\star})\Big{)}\frac{g^{\prime}(F^{\star})}{r_{1}(F^{\star})}\Big{]}$
$\displaystyle\quad+\mathbb{E}\Big{[}r_{1}(F)\mathcal{L}_{\infty}\eta(F)\Big{(}\frac{g^{\prime}(F)}{r_{1}(F)}-\frac{g^{\prime}(F^{\star})}{r_{1}(F^{\star})}\Big{)}\Big{]}$
$\displaystyle\quad+\mathbb{E}\Big{[}\Big{(}r_{2}(F)\Psi_{\infty}\eta(F)-r_{2}(F^{\star})\Psi_{\infty}\eta(F^{\star})\Big{)}\frac{g(F^{\star})}{r_{2}(F^{\star})}\Big{]}$
$\displaystyle\quad+\mathbb{E}\Big{[}r_{2}(F)\Psi_{\infty}\eta(F)\Big{(}\frac{g(F)}{r_{2}(F)}-\frac{g(F^{\star})}{r_{2}(F^{\star})}\Big{)}\Big{]},$

with $r_{1},r_{2}$ two functions left to be determined with the aim of
tempering the malevolent intentions of $g^{\prime}$ and $g^{\prime\prime}$.
These considerations lead to the main result of the Section.

######  Proposition A.1.

Let the previous notations and assumptions prevail. Then

$\displaystyle\left|\mathbb{E}h(F)-\mathbb{E}h(F_{\infty})\right|$
$\displaystyle\leq\kappa_{1}\mathbb{E}\Big{[}\Big{|}r_{1}(F)\mathcal{L}_{\infty}\eta(F)-r_{1}(F^{\star})\mathcal{L}_{\infty}\eta(F^{\star})\Big{|}\Big{]}+\kappa_{2}\mathbb{E}\Big{[}|r_{1}(F)\mathcal{L}_{\infty}\eta(F)||F-F^{\star}|\Big{]}$
$\displaystyle\quad+\kappa_{3}\mathbb{E}\Big{[}\Big{|}r_{2}(F)\Psi_{\infty}\eta(F)-r_{2}(F^{\star})\Psi_{\infty}\eta(F^{\star})\Big{|}\Big{]}+\kappa_{4}\mathbb{E}\Big{[}|r_{2}(F)\Psi_{\infty}\eta(F)||F-F^{\star}|\Big{]}$
(34)

where $\kappa_{j}=\sup_{x}\left|\kappa_{j}(x)\right|$ for $j=1,\ldots,4$ with

$\displaystyle\kappa_{1}(x)=\frac{g^{\prime}(x)}{r_{1}(x)},\kappa_{2}(x)=\left(\frac{g^{\prime}(x)}{r_{1}(x)}\right)^{\prime},\,\kappa_{3}(x)=\frac{g(x)}{r_{2}(x)}\mbox{
and }\kappa_{4}(x)=\left(\frac{g(x)}{r_{2}(x)}\right)^{\prime}.$ (35)

######  Remark A.3.

The functions $r_{1},r_{2}$ and $\eta$ can, for all intents and purposes, be
chosen freely. A good choice of function $\eta$ seems to be
$\eta(x)=\eta_{\gamma}(x)=\gamma^{\prime}(x)/\gamma(x)$, at least if
$\gamma^{\prime}(x)/\gamma(x)$ is non-increasing on $\mathbb{R}$ and
$\mathbb{E}[b^{\prime}(Z)]=0$ when $Z\sim\gamma$. Indeed in this case:

$\displaystyle\mathbb{E}[\eta_{\gamma}(F_{\infty})]=-\mathbb{E}\Big{[}\frac{b^{\prime}(F_{\infty})}{b(F_{\infty})}\Big{]}=-\mathbb{E}[b^{\prime}(Z)]=0\mbox{
and }\overline{\eta}_{\gamma}(x)=\eta_{\gamma}(x)$
$\displaystyle\mathcal{L}_{\infty}\eta_{\gamma}(x)=1-\frac{\mathcal{L}_{\gamma}b^{\prime}(x)}{b(x)}\mbox{
and
}\overline{\eta}_{\gamma}(x)-\frac{\gamma^{\prime}(x)}{\gamma(x)}\mathcal{L}_{\infty}\eta_{\gamma}(x)=\frac{\gamma^{\prime}(x)}{\gamma(x)}\frac{\mathcal{L}_{\gamma}b^{\prime}(x)}{b(x)}.$

Another natural choice (which turns out to be equivalent to the previous one
when $\gamma$ is the Gaussian density) is $\eta(x)=-\mathrm{Id}(x)=-x$ for
which

$\displaystyle-\mathbb{E}[\mathrm{Id}(F_{\infty})]=-\mathbb{E}\left[Zb(Z)\right]=\mathbb{E}[b^{\prime}(Z)]$
$\displaystyle-\mathcal{L}_{\infty}\mathrm{Id}(x)=\tau_{\infty}(x)\quad\mbox{
(the Stein kernel)}.$

Other choices are possible, depending on the properties of the density
$\gamma$; it may be worthwhile investigating this avenue, though we will not
do it here.

####  A.2.2 When the base distribution is standard Gaussian

We now specialize the previous construction to the case that $\gamma(x)$ is
the standard Gaussian density. As before, we suppose that $b$ is chosen in
such a way that $\mathbb{E}[F_{\infty}]=0$; note that if $Z\sim\gamma$ the
standard normal then we also have
$\mathbb{E}[F_{\infty}]=\mathbb{E}[Zb(Z)]=-\mathbb{E}[b^{\prime}(Z)]$. If
$\gamma$ is the Gaussian density then many of the previous expressions
simplify, because $\gamma^{\prime}(x)/\gamma(x)=-\mathrm{Id}(x):=-x$. For
instance $\sigma^{2}_{B}(F)=\mathbb{E}[F{B(F)}/{b(F)}]$ and taking
$\eta=-\mathrm{Id}$ we get $\mathcal{L}_{\gamma}\eta=1$. Also
$\mathcal{L}_{\infty}\eta=\tau_{\infty}$ is now the so-called Stein kernel of
$p_{\infty}$; this function is well known to have very good properties for the
analysis of $p_{\infty}$, see e.g. [5] for an overview. At this stage it
suffices to remark that $\tau_{\infty}(x)\geq 0$ for all $x$. We also have the
nice identity $\Psi_{\infty}\eta(x)=x(\tau_{\infty}(x)-1)$ so that (34)
becomes

$\displaystyle\left|\mathbb{E}h(F)-\mathbb{E}h(F_{\infty})\right|$
$\displaystyle\leq\kappa_{1}\mathbb{E}\Big{[}\Big{|}r_{1}(F)\tau_{\infty}(F)-r_{1}(F^{\star})\tau_{\infty}(F^{\star})\Big{|}\Big{]}+\kappa_{2}\mathbb{E}\Big{[}|r_{1}(F)|\tau_{\infty}(F)|F-F^{\star}|\Big{]}$
$\displaystyle\quad+\kappa_{3}\mathbb{E}\Big{[}\Big{|}Fr_{2}(F)(\tau_{\infty}(F)-1)-F^{\star}r_{2}(F^{\star})(\tau_{\infty}(F^{\star})-1)\Big{|}\Big{]}$
$\displaystyle\quad+\kappa_{4}\mathbb{E}\Big{[}|Fr_{2}(F)(\tau_{\infty}(F)-1)||F-F^{\star}|\Big{]}$
(36)

with $\kappa_{j},j=1,\ldots,4$ still as in (35). The following general result
permits to bound these constants.

######  Lemma A.4.

Let all above notations and assumptions prevail (in particular
$\|h^{\prime}\|\leq 1$). Then

$\displaystyle\kappa_{1}(x)\leq\frac{2}{|r_{1}(x)|}\frac{R_{\infty}(x)}{(\tau_{\infty}(x))^{2}}$
(37)
$\displaystyle\kappa_{2}(x)\leq\frac{2}{|r_{1}(x)|\tau_{\infty}(x)}\left(1+\left|\frac{2x}{\tau_{\infty}(x)}-x+\frac{b^{\prime}(x)}{b(x)}-\frac{r_{1}^{\prime}(x)}{r_{1}(x)}\right|\frac{R_{\infty}(x)}{\tau_{\infty}(x)}\right)$
$\displaystyle\kappa_{3}(x)\leq\frac{1}{|r_{2}(x)|}$ (38)
$\displaystyle\kappa_{4}(x)\leq\frac{1}{|r_{2}(x)|}\left(\frac{2R_{\infty}(x)}{(\tau_{\infty}(x))^{2}}+\left|\frac{r_{2}^{\prime}(x)}{r_{2}(x)}\right|\right)$

where
$R_{\infty}(x)={\int_{-\infty}^{x}{P}_{\infty}(u)\mathrm{d}u\int_{x}^{\infty}\overline{P}_{\infty}(u)\mathrm{d}u}/{p_{\infty}(x)}.$

###### Proof of Lemma A.4.

Let $g$ be defined in (31) with $\eta=-\mathrm{Id}$. Suppose that
$\mathbb{E}b^{\prime}(Z)=0$ and let $h$ be absolutely continuous. We start
with the fact that, from (32):

$\displaystyle\tau_{\infty}(x)g^{\prime\prime}(x)=h^{\prime}(x)-(-x+\tau_{\infty}^{\prime}(x))g^{\prime}(x)+g(x).$

Using

$\displaystyle\tau_{\infty}^{\prime}(x)=\left(x-\frac{b^{\prime}(x)}{b(x)}\right)\tau_{\infty}(x)-x$

we get

$\displaystyle
g^{\prime\prime}(x)=\frac{g(x)+h^{\prime}(x)}{\tau_{\infty}(x)}+\left(\frac{2x}{\tau_{\infty}(x)}-x+\frac{b^{\prime}(x)}{b(x)}\right)g^{\prime}(x).$
(39)

With all this we are ready to bound the different coefficients, whose
expressions we recall for ease of reference:

$\displaystyle\kappa_{1}(x)=\left|\frac{g^{\prime}(x)}{r_{1}(x)}\right|,\kappa_{2}(x)=\left|\left(\frac{g^{\prime}(x)}{r_{1}(x)}\right)^{\prime}\right|,\,\kappa_{3}(x)=\left|\frac{g(x)}{r_{2}(x)}\right|\mbox{
and }\kappa_{4}(x)=\left|\left(\frac{g(x)}{r_{2}(x)}\right)^{\prime}\right|.$

It follows immediately from (28) that

$g(x)=\frac{\mathbb{E}\left[(h(X_{1})-h(X_{2}))\mathbb{I}[X_{1}\leq x\leq
X_{2}]\right]}{\mathbb{E}\left[(X_{2}-X_{1})\mathbb{I}[X_{1}\leq x\leq
X_{2}]\right]}$

at all $x\in\mathcal{S}(p_{\infty})$, where $X_{1}$ and $X_{2}$ are
independent copies of $F_{\infty}$. Since, by assumption,
$|h(x)-h(y)|\leq\left|x-y\right|$, (38) follows. To pursue, we use [5, Lemma
2.25] to deduce

$\displaystyle g^{\prime}(x)$
$\displaystyle=\frac{\bar{h}(x)\,\mathcal{L}_{\infty}\eta(x)-\bar{\eta}(x)\,\mathcal{L}_{\infty}h(x)}{(\mathcal{L}_{\infty}\eta(x))^{2}}$
$\displaystyle=\frac{1}{p_{\infty}(x)\tau_{\infty}(x)^{2}}\left(\mathbb{E}\left[h^{\prime}(F_{\infty})F_{\infty}\frac{\bar{P}_{\infty}(F_{\infty})}{p_{\infty}(X)}\mathbb{I}[x\leq
F_{\infty}]\right]\mathbb{E}\left[\frac{P_{\infty}(F_{\infty})}{p_{\infty}(F_{\infty})}\mathbb{I}[F_{\infty}\leq
x]\right]\right.$
$\displaystyle\qquad\left.-\mathbb{E}\left[h^{\prime}(F_{\infty})\frac{P_{\infty}(F_{\infty})}{p_{\infty}(F_{\infty})}\mathbb{I}[F_{\infty}\leq
x]\right]\mathbb{E}\left[\frac{\bar{P}_{\infty}(F_{\infty})}{p_{\infty}(F_{\infty})}\mathbb{I}[x\leq
X]\right]\right)$

(where $\bar{P}$ is the survival function of cdf $P$). The bound on the
derivative then follows (see also [5, Equation (2.38)]):

$\displaystyle|g^{\prime}(x)|$
$\displaystyle\leq\|h^{\prime}\|\frac{2}{p_{\infty}(x)\tau_{\infty}(x)^{2}}\mathbb{E}\left[\frac{\overline{P}_{\infty}(F_{\infty})}{p_{\infty}(F_{\infty})}\mathbb{I}[x\leq
F_{\infty}]\right]\mathbb{E}\left[\frac{{P}_{\infty}(F_{\infty})}{p_{\infty}(F_{\infty})}\mathbb{I}[F_{\infty}\leq
x]\right]$ $\displaystyle\leq
2\frac{1}{\tau_{\infty}(x)^{2}}\frac{\int_{-\infty}^{x}{P}_{\infty}(u)\mathrm{d}u\int_{x}^{\infty}\overline{P}_{\infty}(u)\mathrm{d}u}{p_{\infty}(x)},$

which brings (37). Furthemore, simply by taking derivatives and using the
previous bounds, we get

$\displaystyle\left|\left(\frac{g(x)}{r_{2}(x)}\right)^{\prime}\right|\leq\frac{1}{|r_{2}(x)|}\left(\frac{2R_{\infty}(x)}{\tau_{\infty}(x)^{2}}+\left|\frac{r_{2}^{\prime}(x)}{r_{2}(x)}\right|\right)$

as well as (using (39) to express $g^{\prime\prime}$ in terms of the lower
order derivatives)

$\displaystyle\left(\frac{g^{\prime}(x)}{r_{1}(x)}\right)^{\prime}\leq\frac{2}{\tau_{\infty}(x)|r_{1}(x)|}\left(1+\left|\frac{2x}{\tau_{\infty}(x)}-x+\frac{b^{\prime}(x)}{b(x)}-\frac{r_{1}^{\prime}(x)}{r_{1}(x)}\right|\frac{R_{\infty}(x)}{\tau_{\infty}(x)}\right).$

All claims are therefore established. ∎

We now apply these results to the choice $B(x)\propto x|x|^{k}/(k+1)$ and
$b(x)\propto|x|^{k}$ with $k\in\mathbb{N}$. Then $\mathbb{E}F_{\infty}=0$ and

$\displaystyle\sigma^{2}_{B}(F)=\mathbb{E}\left[F\frac{B(F)}{b(F)}\right]=\frac{\mathbb{E}F^{2}}{k+1},$

so that our first assumptions become $P(F\neq 0)=1$ and
$\mathbb{E}[F^{2}]=k+1$. The following results then follow from direct
manipulations of the definitions. The first result is the same as Lemma 3.1 in
the main text.

######  Lemma A.5.

Let all above notations prevail, and set $\tau_{\infty}(\cdot;k)$ to be the
Stein kernel of $p_{\infty}(x;,k)=b(x)\gamma(x)$. If $x\neq 0$ then

$\displaystyle\tau_{\infty}(x;k)=2^{k/2}e^{x^{2}/2}|x|^{-k}\Gamma(1+k/2,x^{2}/2)$

where $\Gamma(\alpha,x)=\int_{x}^{\infty}t^{\alpha-1}e^{-t}\mathrm{d}t$ is the
(upper) incomplete gamma function.

Incomplete gamma functions are quite well understood. For instance, using
[11], we readily obtain the next result.

######  Lemma A.6.

The Stein kernel $\tau_{\infty}$ is strictly decreasing on $(0,\infty)$ and
satisfies

$\displaystyle\lim_{x\to
0}|x|^{k}\tau_{\infty}(x;k)=2^{k/2}\Gamma(1+k/2)\mbox{ and
}\lim_{|x|\to\infty}\tau_{\infty}(x;k)=1$

as well as the inequalities

$\displaystyle\frac{1}{|x|^{k}\tau_{\infty}(x;k)}\leq\frac{1}{2^{k/2}\Gamma(1+k/2)},\quad\left|\frac{2}{\tau_{\infty}(x;k)}-1\right|\leq
1$ (40) $\displaystyle\mbox{ and
}\frac{1}{|x|^{k-1}\tau_{\infty}(x;k)}\leq\frac{1}{2^{(k-1)/2}\Gamma(1+k/2)}$

for all $x\in\mathbb{R}$ and all $k\in\mathbb{N}_{0}$. Moreover, if we let
$P_{\infty}(\cdot;k)$ and $\overline{P}_{\infty}(\cdot;k)$ be the cdf and
survival function of $p_{\infty}(\cdot;k)$, and define
$R_{\infty}(x;k)={\int_{-\infty}^{x}{P}_{\infty}(u;k)\mathrm{d}u\int_{x}^{\infty}\overline{P}_{\infty}(u;k)\mathrm{d}u}/{p_{\infty}(x;k)}$
as in Lemma A.4, then

$\displaystyle\frac{R_{\infty}(x;k)}{\tau_{\infty}(x;k)}\leq\frac{\Gamma(k/2+1)}{\sqrt{2}\Gamma(k/2+1/2)}$
(41)

for all $x\in\mathbb{R}$ and all $k\in\mathbb{N}_{0}$.

######  Remark A.7.

The bounds (40) and (41) are sharp because they are attained at $x\to 0$.

Using Lemmas A.5 and A.6 we obtain the required bounds on the constants
$\kappa_{j}$.

######  Corollary A.1.

Set $r_{1}(x)=x^{k}$ and $r_{2}(x)=1$ (the sign of $x$) in Lemma A.4. Then

$\displaystyle\kappa_{1}(x)\leq\frac{2}{|x|^{k}}\frac{R_{\infty}(x;k)}{(\tau_{\infty}(x;k))^{2}}\leq\frac{2^{(1-k)/2}}{\Gamma((k+1)/2)}$
(42)
$\displaystyle\kappa_{2}(x)\leq\frac{2}{|x|^{k}\tau_{\infty}(x;k)}+\frac{2}{|x|^{k-1}\tau_{\infty}(x;k)}\left|\frac{2}{\tau_{\infty;k}(x)}-1\right|\frac{R_{\infty}(x;k)}{\tau_{\infty}(x;k)}\leq
3\frac{2^{-k/2}}{\Gamma((1+k)/2)}$ (43) $\displaystyle\kappa_{3}(x)\leq 1$
(44)
$\displaystyle\kappa_{4}(x)\leq\frac{2R_{\infty}(x;k)}{(\tau_{\infty}(x;k))^{2}}\leq
1$ (45)

######  Remark A.8.

The bounds in (42) and (44) are sharp; the other two are not.

####  A.2.3 When $k=0,1$ and a proof of Proposition 3.3

Upper bounding (42) by 1 and (43) by 2 (neither of these choices, nor the
bound 1 in (45), are optimal in $k$ because the true value goes to 0 as $k$
goes to $\infty$), inequality (36) leads to the following result.

######  Theorem A.9.

If $F_{\infty}\sim{\mathbb{F}}_{\infty}$ has density
$p_{\infty}(x)\propto|x|^{k}\varphi(x)$, for a given $k\in\mathbb{N}$ and
$F\sim{\mathbb{F}}$ is some random variable with mean 0 such that $P(F\neq
0)=1$ and $\mathbb{E}[F^{2}]=k+1$ then there exists a random variable
$F^{\star}$, called the generalized $k$-radial-bias transform of $F$, which
uniquely satisfies

$\mathbb{E}\left[\frac{f^{\prime}(F^{\star})}{|F^{\star}|^{k}}\right]=\mathbb{E}\left[\frac{f(F)}{|F|^{k-1}}\right]$

for all $f$ such that both integrals exist, and

$\displaystyle\mathop{d_{\mathrm{W}}}({\mathbb{F}},{\mathbb{F}}_{\infty})$
$\displaystyle\leq\mathbb{E}\Big{|}F^{k}\tau_{\infty}(F;k)-(F^{\star})^{k}\tau_{\infty}(F^{\star};k)\Big{|}+2\mathbb{E}\left[|F|^{k}\tau_{\infty}(F;k)|F-F^{\star}|\right]$
$\displaystyle\quad+\mathbb{E}\Big{|}F(\tau_{\infty}(F;k)-1)-F^{\star}(\tau_{\infty}(F^{\star};k)-1)\Big{|}$
$\displaystyle\quad+\mathbb{E}\left[\big{|}F(\tau_{\infty}(F;k)-1)\big{|}|F-F^{\star}|\right]$
(46)

where $\tau_{\infty}(x;k)$ is the Stein kernel given in (21).

The upper bounds in Propositions 3.2 and 3.3 from the main text are direct
corollaries of this result. We have already proved the upper bound from
Proposition 3.2 in Example 3.2 by other means. We therefore concentrate on
Proposition 3.3.

The following can be shown directly from the definitions:

$\displaystyle(x\tau_{\infty}(x,1))^{\prime}\leq
1,\quad|x(\tau_{\infty}(x,1)-1)|\leq 1\mbox{ and
}|\left(x(\tau_{\infty}(x,1)-1)\right)^{\prime}|\leq 1.$

Plugging these into (46) gives

$\displaystyle\mathop{d_{\mathrm{W}}}({\mathbb{F}},{\mathbb{F}}_{\infty})$
$\displaystyle\leq
3\mathbb{E}\left[|F-F^{\star}|\right]+2\mathbb{E}\left[(|F|(\tau_{\infty}(F;1)-1)|F-F^{\star}|\right]+2\mathbb{E}\left[|F||F-F^{\star}|\right]$
$\displaystyle=\mathbb{E}\left[\left(5+2|F|\right)|F-F^{\star}|\right]$
$\displaystyle\leq(5+2x_{1})(N-1)^{-1}\sum_{i=1}^{N-1}|x_{i+1}-x_{i}|$
$\displaystyle=(5+2x_{1})2{x_{1}}/(N-1)$ $\displaystyle=O(\log N/N),$

where we have used a $k=1$ version of the coupling argument given in [16,
Proof of Corollary 3.7] along with a version of [16, Lemma 4.8] showing that
$x_{1}=O(\sqrt{\log N})$, as required.

####  A.2.4 When $k\geq 2$

We now focus on the case $k\geq 2$. We once again call upon [11] to obtain the
next lemma.

######  Lemma A.10 (Stein kernel).

Let all previous notations prevail. Given $j,k$ two integers define

$a_{j}(k)={2^{j}}\frac{\Gamma(1+k/2)}{\Gamma(1+k/2-j)}$

with the convention that $a_{j}(k)=0$ for all $j\geq k$. Then, for all
$x\in\mathbb{R}$ and all $k\in\mathbb{N}$, we have

$\tau_{\infty}(x;k)=\sum_{j=0}^{\lfloor\frac{k}{2}\rfloor}\frac{a_{j}(k)}{x^{2j}}+\frac{a_{\lceil
k/2\rceil}(k)}{\sqrt{2}}\epsilon_{k}(x)$

where

$\epsilon_{k}(x)=\begin{cases}0&\mbox{ if }k=2\ell\\\
e^{x^{2}/2}|x|^{-(2\ell-1)}\Gamma(1/2,x^{2}/2)&\mbox{ if
}k=2\ell-1\end{cases}$

Moreover the remainder satisfies

$0\leq\epsilon_{k}(x)\leq
2^{(k+1)/2}x^{-(k+1)}\frac{\Gamma(1+k/2)}{\Gamma(1/2)}$

for all $k\in\mathbb{N}$ and all $x$.

###### Proof.

The claim follows from the following representation of the incomplete gamma
function (available e.g. from [11, Theorem 3 and Proposition 13]): for all
$a>0$ and all $x>0$,

$\displaystyle\Gamma(a,x)=e^{-x}x^{a-1}\sum_{j=0}^{\lfloor
a\rfloor-1}P_{j}(a)x^{-j}+r(a,x)$

where $P_{j}(a)={\Gamma(a)}/{\Gamma(a-j)}$ (and $P_{j}(a)=0$ for all $j\geq
a$), and, setting $[a]=a-\lfloor a\rfloor$, $r(a,x)=P_{\lfloor
a\rfloor}(a)\Gamma([a],x)$ which satisfies

$\displaystyle 0\leq r(a,x)\leq e^{-x}P_{\lfloor a\rfloor}(a)x^{[a]-1}.$

The claim follows.

∎

In [16] we considered $b(x)=x^{2}$. The argument from that paper is now
extended to arbitrary integer $k$ in the next result.

######  Theorem A.11.

Instate all previous notations and let

$a_{r}(k)=\frac{\Gamma(k/2+1)}{\Gamma(k/2-r+1)}2^{r}.$

It holds that if $k=2\ell$ is an even integer then

$\displaystyle d_{W}(F,F_{\infty})$
$\displaystyle\leq\sum_{j=0}^{\ell}a_{\ell-j}(2\ell)\Bigg{(}\mathbb{E}\left|F^{2j}-(F^{\star})^{2j}\right|+2\mathbb{E}\
\left[|F|^{2j}\left|F-F^{\star}\right|\right]\Bigg{)}$
$\displaystyle\quad\qquad+\sum_{j=1}^{\ell}a_{j}(2\ell)\Bigg{(}\mathbb{E}\left|\frac{1}{F^{2j-1}}-\frac{1}{(F^{\star})^{2j-1}}\right|+\mathbb{E}\left[\frac{1}{|F|^{2j-1}}\Big{|}F-F^{\star}\Big{|}\right]\Bigg{)}$

and, if $k=2\ell-1$ is an odd integer, then

$\displaystyle d_{W}(F,F_{\infty})$
$\displaystyle\leq\sum_{j=1}^{\ell}a_{\ell-j}(2\ell-1)\Bigg{(}\mathbb{E}\left|F^{2j-1}-(F^{\star})^{2j-1}\right|+2\mathbb{E}\
\left[|F|^{2j-1}\left|F-F^{\star}\right|\right]\Bigg{)}$
$\displaystyle\quad\qquad+\sum_{j=1}^{\ell-1}a_{j}(2\ell-1)\Bigg{(}\mathbb{E}\left|\frac{1}{F^{2j-1}}-\frac{1}{(F^{\star})^{2j-1}}\right|+\mathbb{E}\left[\frac{1}{|F|^{2j-1}}\Big{|}F-F^{\star}\Big{|}\right]\Bigg{)}$
$\displaystyle\quad\qquad+3a_{\ell}(2\ell-1)\mathbb{E}\left[\left(2+\frac{2}{|F|^{2(\ell-1)}}+\frac{1}{|F^{\star}|^{2(\ell-1)}}\right)\left|F-F^{\star}\right|\right]$
$\displaystyle\quad\qquad+3a_{\ell}(2\ell-1)\mathbb{E}\left[\left(|F|+|F|^{2\ell-1}\right)\left|\frac{1}{F^{2\ell-1}}-\frac{1}{(F^{\star})^{2\ell-1}}\right|\right]$

###### Proof.

The claim for even integers $k$ is immediate. If $k=2\ell-1$ is an odd integer
then

$\displaystyle d_{W}(F,F_{\infty})$
$\displaystyle\leq\sum_{j=1}^{\ell}a_{\ell-j}(2\ell-1)\Bigg{(}\mathbb{E}\left|F^{2j-1}-(F^{\star})^{2j-1}\right|+2\mathbb{E}\
\left[|F|^{2j-1}\left|F-F^{\star}\right|\right]\Bigg{)}$
$\displaystyle\quad\qquad+\sum_{j=1}^{\ell-1}a_{j}(2\ell-1)\Bigg{(}\mathbb{E}\left|\frac{1}{F^{2j-1}}-\frac{1}{(F^{\star})^{2j-1}}\right|+\mathbb{E}\left[\frac{1}{|F|^{2j-1}}\Big{|}F-F^{\star}\Big{|}\right]\Bigg{)}$
$\displaystyle\quad\qquad+\frac{a_{\ell}(2\ell-1)}{\sqrt{2}}\Psi_{\ell}(F,F^{\star})$

where

$\displaystyle\Psi_{\ell}(F,F^{\star})$
$\displaystyle=\mathbb{E}\left|F^{2\ell-1}\epsilon_{2\ell-1}(F)-(F^{\star})^{2\ell-1}\epsilon_{2\ell-1}(F^{\star})\right|+2\mathbb{E}\left[|F|^{2\ell-1}\epsilon_{2\ell-1}(F)|F-F^{\star}|\right]$
$\displaystyle\quad+\mathbb{E}\left|F\epsilon_{2\ell-1}(F)-(F^{\star})\epsilon_{2\ell-1}(F^{\star})\right|+\mathbb{E}\left[|F|\epsilon_{2\ell-1}(F)|F-F^{\star}|\right]$

Here we aim to bound

$\displaystyle\Psi_{\ell}(F,F^{\star})$
$\displaystyle=\mathbb{E}\left|F^{2\ell-1}\epsilon_{2\ell-1}(F)-(F^{\star})^{2\ell-1}\epsilon_{2\ell-1}(F^{\star})\right|+2\mathbb{E}\left[|F|^{2\ell-1}\epsilon_{2\ell-1}(F)|F-F^{\star}|\right]$
$\displaystyle\quad+\mathbb{E}\left|F\epsilon_{2\ell-1}(F)-(F^{\star})\epsilon_{2\ell-1}(F^{\star})\right|+\mathbb{E}\left[|F|\epsilon_{2\ell-1}(F)|F-F^{\star}|\right]$
$\displaystyle=:I+II+III+IV$

where $\epsilon_{2\ell-1}(x)=e^{x^{2}/2}|x|^{-(2\ell-1)}\Gamma(1/2,x^{2}/2)$.
We first note that, for $x>0$, the function
$x\mapsto\nu(x):=e^{x^{2}/2}\Gamma(1/2,x^{2/2})$ is strictly decreasing as
$|x|\to\infty$, with maximum value $\sqrt{\pi}$ at $x=0$. Hence

$\displaystyle II+IV$ $\displaystyle\leq
2\sqrt{\pi}\mathbb{E}\left[|F-F^{\star}|\right]+\sqrt{\pi}\mathbb{E}\left[\frac{1}{|F|^{2(\ell-1)}}|F-F^{\star}|\right].$

Also, $|\nu^{\prime}(x)|\leq\sqrt{2}$ so that
$|\nu(x)-\nu(y)|\leq\sqrt{2}|y-x|$, hence

$\displaystyle\left|\epsilon_{2\ell-1}(x)-\epsilon_{2\ell-1}(y)\right|$
$\displaystyle\leq|x|^{-(2\ell-1)}|\nu(x)-\nu(y)|+|x^{-(2\ell-1)}-y^{-(2\ell-1)}|\nu(y)$
$\displaystyle\leq\sqrt{2}|x|^{-(2\ell-1)}|x-y|+\sqrt{\pi}|x^{-(2\ell-1)}-y^{-(2\ell-1)}|.$

This gives

$\displaystyle I$
$\displaystyle\leq\mathbb{E}\left[\left|F\right|^{2\ell-1}\left|\epsilon_{2\ell-1}(F)-\epsilon_{2\ell-1}(F^{\star})\right|\right]+\mathbb{E}\left[\left|F^{2\ell-1}-(F^{\star})^{2\ell-1}\right|\epsilon_{2\ell-1}(F^{\star})\right]$
$\displaystyle\leq\mathbb{E}\left[\left|F\right|^{2\ell-1}\left|\epsilon_{2\ell-1}(F)-\epsilon_{2\ell-1}(F^{\star})\right|\right]+\sqrt{\pi}\mathbb{E}\left[\frac{1}{|F^{\star}|^{2\ell-1}}\left|F^{2\ell-1}-(F^{\star})^{2\ell-1}\right|\right]$
$\displaystyle\leq\sqrt{2}\mathbb{E}\left[\left|F-F^{\star}\right|\right]+\sqrt{\pi}\mathbb{E}\left[|F|^{2\ell-1}\left|F^{-(2\ell-1)}-(F^{\star})^{-(2\ell-1)}\right|\right]$
$\displaystyle\qquad+\sqrt{\pi}\mathbb{E}\left[\frac{1}{|F^{\star}|^{2\ell-1}}\left|F^{2\ell-1}-(F^{\star})^{2\ell-1}\right|\right].$

Similarly,

$\displaystyle III$
$\displaystyle\leq\mathbb{E}\left[\left|F\right|\left|\epsilon_{2\ell-1}(F)-\epsilon_{2\ell-1}(F^{\star})\right|\right]+\mathbb{E}\left[\left|F-F^{\star}\right|\epsilon_{2\ell-1}(F^{\star})\right]$
$\displaystyle\leq\mathbb{E}\left[\left|F\right|\left|\epsilon_{2\ell-1}(F)-\epsilon_{2\ell-1}(F^{\star})\right|\right]+\sqrt{\pi}\mathbb{E}\left[\frac{1}{|F^{\star}|^{2\ell-1}}\left|F-F^{\star}\right|\right]$
$\displaystyle\leq\sqrt{2}\mathbb{E}\left[\left|F\right|^{-2(\ell-1)}\left|F-F^{\star}\right|\right]+\sqrt{\pi}\mathbb{E}\left[|F|\left|F^{-(2\ell-1)}-(F^{\star})^{-(2\ell-1)}\right|\right]$
$\displaystyle\quad+\sqrt{\pi}\mathbb{E}\left[\frac{1}{|F^{\star}|^{2\ell-1}}\left|F-F^{\star}\right|\right].$

Combining these bounds leads to

$\displaystyle\Psi_{\ell}(F,F^{\star})$ $\displaystyle\leq
2\sqrt{\pi}\mathbb{E}\left[|F-F^{\star}|\right]+\sqrt{\pi}\mathbb{E}\left[\frac{1}{|F|^{2(\ell-1)}}|F-F^{\star}|\right]$
$\displaystyle\qquad+\sqrt{2}\mathbb{E}\left[\left|F-F^{\star}\right|\right]+\sqrt{\pi}\mathbb{E}\left[|F|^{2\ell-1}\left|F^{-(2\ell-1)}-(F^{\star})^{-(2\ell-1)}\right|\right]$
$\displaystyle\qquad+\sqrt{\pi}\mathbb{E}\left[\frac{1}{|F^{\star}|^{2\ell-1}}\left|F^{2\ell-1}-(F^{\star})^{2\ell-1}\right|\right]$
$\displaystyle\qquad+\sqrt{2}\mathbb{E}\left[\left|F\right|^{-2(\ell-1)}\left|F-F^{\star}\right|\right]+\sqrt{\pi}\mathbb{E}\left[|F|\left|F^{-(2\ell-1)}-(F^{\star})^{-(2\ell-1)}\right|\right]$
$\displaystyle\qquad+\sqrt{\pi}\mathbb{E}\left[\frac{1}{|F^{\star}|^{2\ell-1}}\left|F-F^{\star}\right|\right]$

which, after bounding all constants by $3\sqrt{2}$ for simplicity, gives

$\displaystyle\Psi_{\ell}(F,F^{\star})$ $\displaystyle\leq
3\sqrt{2}\mathbb{E}\left[\left(2+\frac{2}{|F|^{2(\ell-1)}}+\frac{1}{|F^{\star}|^{2(\ell-1)}}\right)\left|F-F^{\star}\right|\right]$
$\displaystyle+3\sqrt{2}\mathbb{E}\left[\left(|F|+|F|^{2\ell-1}\right)\left|F^{-(2\ell-1)}-(F^{\star})^{-(2\ell-1)}\right|\right],$

which is the claim. ∎

######  Corollary A.2.

Instate all notations from Theorem A.11. The following bounds hold:

  * •

$k=2$ (Maxwell case): $p_{\infty}(x)=x^{2}\varphi(x)$, $\mathbb{E}[F^{2}]=3$,
$\mathbb{E}[f^{\prime}(F^{\star})/(F^{\star})^{2}]=\mathbb{E}[f(F)/|F|]$ for
all $f$, and

$\displaystyle d_{W}(F,F_{\infty})$
$\displaystyle\leq\mathbb{E}\left[\left(3+2|F|+\frac{2}{|F|}+\frac{2}{|F||F^{\star}|}\right)|F-F^{\star}|\right].$

  * •

$k=3$: $p_{\infty}(x)\propto|x|^{3}\varphi(x)$, $\mathbb{E}[F^{2}]=4$,
$\mathbb{E}[f^{\prime}(F^{\star})/|F^{\star}|^{3}]=\mathbb{E}[f(F)/F^{2}]$ for
all $f$, and

$\displaystyle d_{W}(F,F_{\infty})$
$\displaystyle\leq\mathbb{E}\left[\left(21+6|F|+2F^{2}+\frac{3}{|F|}+\frac{18}{F^{2}}+\frac{9}{(F^{\star})^{2}}+\frac{3}{|F||F^{\star}|}\right)|F-F^{\star}|\right]$
$\displaystyle\quad+\mathbb{E}|F^{2}-(F^{\star})^{2}|+9\mathbb{E}\left[(|F|+|F|^{3})\left|\frac{1}{F^{3}}-\frac{1}{(F^{\star})^{3}}\right|\right].$

######  Remark A.12.

We note that the Maxwell bound ($k=2$) is the same as in [16, Equation (24)]
but with improved constants.

###  A.3 Upper bounds on the rate of convergence via coupling

In this section we apply Theorem A.11 to give explicit upper bounds on the
accuracy of $F=F_{N}\sim{\mathbb{F}}_{N}$ (defined after Lemma 3.3 of the main
text) in approximating the radial distributions with density
$p_{\infty}(x)\propto|x|^{k}\varphi(x)$ for $k=2,\ldots,14$.

Two crucial steps in controlling the various “error” terms in Theorem A.11
when $k\geq 1$ (in which case there is a singularity in $p_{\infty}$ at the
origin) are to obtain good approximations to $x_{1}$ and $x_{m}$ (where
$m=N/2$), which (when necessary) we write as $x_{1,k}$ and $x_{m,k}$, to
reflect their dependence on $k$. It is easily shown by extending [16, Lemma
4.8] that

$\displaystyle x_{1,k}=O(\sqrt{\log N})$

for all $k\geq 0$. Concerning $x_{m,k}$, we have $x_{m,k}\geq 1/\sqrt{N}$ for
all $k\geq 1$, cf. [16, proof of Corollary 3.7]. However, for $k\geq 3$ a
lower bound of order $1/\sqrt{N}$ is too small to control all error terms in
Theorem A.11, and it seems very challenging to improve this lower bound
analytically.

A way around this difficulty is to examine the numerical behavior of
$x_{m,k}$. We find that the following scaling relation provides a remarkably
accurate approximation:

$x_{m,k}\asymp N^{-1/r_{k}},\ \ \ {\rm where}\ r_{k}=3/2+4k/5$ (47)

for $k=2,3,\ldots,14$, see Figure 8. We expect this scaling relation hold for
all $k\geq 2$.

Figure 8: The graph of $r_{k}=3/2+4k/5$ (blue solid line) obtained by fitting
the scaling relation (47) for $k=2,\ldots,14$. The fitted line approximates
data points (circles) derived from estimates of the slope when regressing
$\log x_{m,k}$ against $\log N$, $N=14,24,\ldots,114$, separately for each
value of $k$. It is remarkable that there is virtually no scatter around any
of the linear fits. For $k$ odd, the fitted line needs to fall above the
diagonal (red dashed) line to ensure adequate control of (48) when $l=k$. For
$k$ even, the fitted line needs to fall above the dot-dashed (cyan) line to
obtain control when $l=k-1$.

The next step is to use the above scaling relation to obtain upper bounds on
the error terms in Theorem A.11 of the form

$\mathbb{E}\left|\frac{1}{F^{l}}-\frac{1}{(F^{\star})^{l}}\right|,\ \ 1\leq
l\leq k$ (48)

with $F\sim\mathbb{F}_{N}$, the empirical distribution of
$x_{1}>\ldots>x_{N}$, and $F^{\star}\sim\mathbb{F}^{\star}_{N}$ is the
corresponding $k$-radial-bias distribution, as given in the following lemma,
which is a consequence of [16, Proposition 3.5].

######  Lemma A.13.

The $k$-radial-bias distribution ${\mathbb{F}}^{\star}_{N}$ of
$\mathbb{F}_{N}$ is defined, and has density

$p^{\star}(x)\propto|x|^{k}\left[\sum_{i=1}^{n}\frac{x_{i}}{|x_{i}|^{k}}\right]$

for $x_{n+1}<x\leq x_{n}$ ($n=1,\ldots,N-1$), and $p^{\star}(x)=0$ if
$x>x_{1}$ or $x\leq x_{N}$.

From Lemma A.13 and the recursion satisfied by $x_{1}>\ldots>x_{N}$, it
follows that $p^{\star}(x)$ puts mass $1/(N-1)$ on each interval between
successive $x_{n}$, so there exists a coupling of $F\sim\mathbb{P}_{N}$ with
$F^{\star}\sim p^{\star}(x)$ such that

$|F-F^{\star}|\leq|x_{n}-x_{n+1}|$

when $F\in[x_{n+1},x_{n}]$. For a detailed proof of such a coupling, see the
construction given in [15]. Now decompose (48) as

$\mathbbm{E}\left|\frac{1}{F^{l}}-\frac{1}{(F^{\star})^{l}}\right|=\mathbbm{E}\left|\frac{1}{F^{l}}-\frac{1}{(F^{\star})^{l}}\right|1_{F^{\star}\in(x_{m+1},x_{m}]}+2\sum_{n=1}^{m-1}\mathbbm{E}\left|\frac{1}{F^{l}}-\frac{1}{(F^{\star})^{l}}\right|1_{F^{\star}\in(x_{n+1},x_{n}]}.$

From Proposition A.13 note that $p^{\star}(x)\propto|x|^{k}$ for
$x\in(x_{m+1},x_{m}]$. Using the fact that $p^{\star}(x)$ puts mass $1/(N-1)$
on this interval, the first term above can be written

$\frac{2(k+1)}{x_{m}^{k+1}(N-1)}\int_{0}^{x_{m}}\left({1\over x^{l}}-{1\over
x_{m}^{l}}\right)x^{k}\,dx\asymp\frac{1}{x_{m}^{l}N}\asymp N^{l/r_{k}-1}\to
0,$

provided $l<3/2+4k/5$ by (47). The second term is bounded above by the
telescoping sum

${2\over N-1}\sum_{n=1}^{m-1}\left({1\over x_{n+1}^{l}}-{1\over
x_{n}^{l}}\right)={2\over N-1}\left({1\over x_{m}^{l}}-{1\over
x_{1}^{l}}\right)=O\left(N^{l/r_{k}-1}\right),$

so we conclude

$\mathbbm{E}\left|\frac{1}{F^{l}}-\frac{1}{(F^{\star})^{l}}\right|=O\left(N^{l/r_{k}-1}\right).$

This bound gives the desired control of (48) for any $l\leq k\leq 7$, see
Figure 8. For even $k$, we only need to consider $l\leq k-1$, so we have
control for $k=8,10$ and $12$, as well.

##  Appendix B Codes

###  B.1 Mathematica code

See McKS2021-Mathematica.pdf

###  B.2 R Code

See McKS2021-Rfile.pdf