# A note on concatenation of quasi-Monte Carlo and plain Monte Carlo rules in
high dimensions††thanks: Keywords: quasi-Monte Carlo, plain Monte Carlo, high-
dimensional integration, weighted function space, rank-1 lattice rule

Takashi Goda School of Engineering, University of Tokyo, 7-3-1 Hongo, Bunkyo-
ku, Tokyo 113-8656, Japan<EMAIL_ADDRESS>
###### Abstract

In this note, we study a concatenation of quasi-Monte Carlo and plain Monte
Carlo rules for high-dimensional numerical integration in weighted function
spaces. In particular, we consider approximating the integral of periodic
functions defined over the $s$-dimensional unit cube by using rank-1 lattice
point sets only for the first $d\,(<s)$ coordinates and random points for the
remaining $s-d$ coordinates. We prove that, by exploiting a decay of the
weights of function spaces, almost the optimal order of the mean squared
worst-case error is achieved by such a concatenated quadrature rule as long as
$d$ scales at most linearly with the number of points. This result might be
useful for numerical integration in extremely high dimensions, such as partial
differential equations with random coefficients for which even the standard
fast component-by-component algorithm is considered computationally expensive.

##  1 Introduction

We study numerical integration of functions defined over the $s$-dimensional
unit cube with large $s$. For an integrable function
$f:[0,1)^{s}\to\mathbb{R}$, we write

$I(f):=\int_{[0,1)^{s}}f(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}.$

For a point set $P\subset[0,1)^{s}$, the integral $I(f)$ is approximated by an
equally-weighted quadrature rule

$I(f;P):=\frac{1}{|P|}\sum_{\boldsymbol{x}\in P}f(\boldsymbol{x}).$

We call $I(f;P)$ a quasi-Monte Carlo (QMC) rule when the point set $P$ is
chosen in a deterministic manner, while we call $I(f;P)$ a plain Monte Carlo
(MC) rule when each point $\boldsymbol{x}\in P$ is generated independently
from the uniform distribution over the domain $[0,1)^{s}$. In this paper we
consider an intermediate rule between QMC and MC rules, in the sense that the
projection of $P$ onto the first $d$ coordinates with some $d<s$ is given by a
deterministic QMC point set, whereas the projection onto the remaining $s-d$
coordinates is simply a set of random points. That is, we _concatenate_ QMC
and MC rules.

Our motivation behind introducing such a concatenated rule mainly comes from
applications to extremely high-dimensional problems such as partial
differential equations with random coefficients [12, 11, 8, 6]. For such
applications, an infinite Karhunen-Loève expansion needs to be truncated up to
some finite terms to perform computations, whose number is the problem
dimension $s$. Therefore, $s$ can be very large in order to keep the
truncation error small. Even so, the notion of weighted function spaces from
[21] can be successfully applied, often resulting in a dimension-independent
error bound supported by the (fast) component-by-component (CBC) construction
of QMC point sets [20, 10, 15, 2].

Typically, for the case of product weights, the fast CBC construction of
rank-1 lattice point sets, a special class of QMC point sets, requires
$O(sN\log N)$ arithmetic operations with $O(N)$ memory [15], where $N$ denotes
the number of points. When $s$ is as large as or even much larger than $N$, we
may need some techniques to reduce the construction cost from the standard
fast CBC algorithm. Some relevant ideas in this direction have been introduced
in [4, 3]. For instance, in [4], the authors assume a fast decay of the weight
parameters and introduce a fast CBC algorithm which reduces the size of the
search space for coordinates with less relative importance. The search space
finally contains only one possible choice after some coordinates, so that we
do not need to search any more. This idea has been further explored in [9] as
a truncation algorithm, in which the latter coordinates are fixed at an anchor
$c\in[0,1)$. In this paper, instead of reducing the size of the search space
or fixing at an anchor, we apply the standard fast CBC algorithm up to the
first $d$ coordinates, and then the remaining $s-d$ coordinates with less
relative importance are filled by random points. Hence the necessary
construction cost is of $O(dN\log N)$, which may lead to a substantial cost
saving if $d$ can be much smaller than $s$.

Here we point out that the idea of concatenating quasi-Monte Carlo and plain
Monte Carlo rules goes back to the work of Spanier [22] and then has been
investigated in terms of discrepancy of concatenated, or mixed/hybrid, point
sets and sequences, see [16, 7, 14, 1] among others. In [17, 18], Owen studied
the variance of concatenated rules between quasi-Monte Carlo and Latin
hypercube samplings. We also refer to [23, Section 4.5] for some results on
concatenating two randomized sequences.

The rest of this note is organized as follows. In the next section, we
introduce the function space of our interest, namely, the weighted Korobov
space consisting of smooth periodic functions over $[0,1)^{s}$. In Section 3,
we describe our concatenated rule made up of a rank-1 lattice point set and
random points. As the main result of this note, we show in Section 4 that
almost the optimal order of the mean squared worst-case error is achieved if
$d$ scales at most linearly with the number of points. We conclude this paper
with comparing the result for our concatenated rule with the results for the
reduced CBC algorithm and for a truncation algorithm, respectively.

##  2 Weighted Korobov spaces

Let $f:[0,1)^{s}\to\mathbb{R}$ be periodic and given by its absolutely
convergent Fourier series

$f(\boldsymbol{x})=\sum_{\boldsymbol{h}\in\mathbb{Z}^{s}}\hat{f}(\boldsymbol{h})\exp\left(2\pi
i\boldsymbol{h}\cdot\boldsymbol{x}\right),$

where the dot product $\cdot$ denotes the usual inner product of two vectors
on the Euclidean space $\mathbb{R}^{s}$ and $\hat{f}(\boldsymbol{h})$ denotes
the $\boldsymbol{h}$-th Fourier coefficient of $f$:

$\hat{f}(\boldsymbol{h}):=\int_{[0,1)^{s}}f(\boldsymbol{x})\exp\left(-2\pi
i\boldsymbol{h}\cdot\boldsymbol{x}\right)\,\mathrm{d}\boldsymbol{x}.$

We measure the smoothness of periodic functions by a parameter $\alpha>1/2$. A
non-increasing sequence of _weights_ $\gamma_{1}\geq\gamma_{2}\geq\cdots>0$
plays a role in moderating the relative importance of different variables and
we write
$\boldsymbol{\gamma}=(\gamma_{1},\gamma_{2},\ldots)\in\mathbb{R}_{>0}^{\mathbb{N}}$.
For a vector $\boldsymbol{h}\in\mathbb{Z}^{s}$, we define

$r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{h}):=\prod_{\begin{subarray}{c}j=1\\\
h_{j}\neq 0\end{subarray}}^{s}\frac{\gamma_{j}}{|h_{j}|^{2\alpha}},$

where the empty product is set to 1. Then the weighted Korobov space, denoted
by $H_{s,\alpha,\boldsymbol{\gamma}}$, is a reproducing kernel Hilbert space
with the reproducing kernel

$K_{s,\alpha,\boldsymbol{\gamma}}(\boldsymbol{x},\boldsymbol{y})=1+\sum_{\boldsymbol{h}\in\mathbb{Z}^{s}\setminus\\{\boldsymbol{0}\\}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{h})\exp\left(2\pi
i\boldsymbol{h}\cdot(\boldsymbol{x}-\boldsymbol{y})\right).$

and the inner product

$\langle
f,g\rangle_{s,\alpha,\boldsymbol{\gamma}}=\hat{f}(\boldsymbol{0})\overline{\hat{g}(\boldsymbol{0})}+\sum_{\boldsymbol{h}\in\mathbb{Z}^{s}\setminus\\{\boldsymbol{0}\\}}\frac{\hat{f}(\boldsymbol{h})\overline{\hat{g}(\boldsymbol{h})}}{r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{h})}.$

The parameter $\alpha$ not only moderates the decay of the Fourier
coefficients, but also coincides precisely with the number of available
square-integrable partial mixed derivatives in each variable when it is an
integer [5, Section 5.8].

The worst-case error of an equally-weighted quadrature rule with a fixed point
set $P$ for the weighted Korobov space is defined by

$e^{\mathrm{wor}}(H_{s,\alpha,\boldsymbol{\gamma}},P):=\sup_{\begin{subarray}{c}f\in
H_{s,\alpha,\boldsymbol{\gamma}}\\\ \|f\|_{s,\alpha,\boldsymbol{\gamma}}\leq
1\end{subarray}}\left|I(f;P)-I(f)\right|.$

It is well-known (see, for instance, [5, Theorem 3.5]) that we have

$\displaystyle(e^{\mathrm{wor}}(H_{s,\alpha,\boldsymbol{\gamma}},P))^{2}$
$\displaystyle=\int_{[0,1)^{s}}\int_{[0,1)^{s}}K_{s,\alpha,\boldsymbol{\gamma}}(\boldsymbol{x},\boldsymbol{y})\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}\boldsymbol{y}$
$\displaystyle\quad-\frac{2}{|P|}\sum_{\boldsymbol{x}\in
P}\int_{[0,1)^{s}}K_{s,\alpha,\boldsymbol{\gamma}}(\boldsymbol{x},\boldsymbol{y})\,\mathrm{d}\boldsymbol{y}+\frac{1}{|P|^{2}}\sum_{\boldsymbol{x},\boldsymbol{y}\in
P}K_{s,\alpha,\boldsymbol{\gamma}}(\boldsymbol{x},\boldsymbol{y})$
$\displaystyle=-1+\frac{1}{|P|^{2}}\sum_{\boldsymbol{x},\boldsymbol{y}\in
P}K_{s,\alpha,\boldsymbol{\gamma}}(\boldsymbol{x},\boldsymbol{y}),$ (1)

where the second equality comes from the equality

$\int_{[0,1)^{s}}K_{s,\alpha,\boldsymbol{\gamma}}(\boldsymbol{x},\boldsymbol{y})\,\mathrm{d}\boldsymbol{y}=1,$

which holds for any $\boldsymbol{x}\in[0,1)^{s}$. As a reference value, the
initial error is defined and given by

$e^{\mathrm{wor}}(H_{s,\alpha,\boldsymbol{\gamma}},0):=\sup_{\begin{subarray}{c}f\in
H_{s,\alpha,\boldsymbol{\gamma}}\\\ \|f\|_{s,\alpha,\boldsymbol{\gamma}}\leq
1\end{subarray}}\left|I(f)\right|=\left(\int_{[0,1)^{s}}\int_{[0,1)^{s}}K_{s,\alpha,\boldsymbol{\gamma}}(\boldsymbol{x},\boldsymbol{y})\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}\boldsymbol{y}\right)^{1/2}=1.$

##  3 Concatenation of QMC and MC rules

As stated in the introduction, we consider using a rank-1 lattice point set
for the first $d$ coordinates with some $0<d<s$ and random points for the
remaining $s-d$ coordinates. First, rank-1 lattice point sets are defined as
follows.

######  Definition 3.1 (Rank-1 lattice point sets).

Let $N,d$ be positive integers and $\boldsymbol{z}\in\\{1,\ldots,N-1\\}^{d}$.
The rank-1 lattice point set determined by $N,d$ and $\boldsymbol{z}$ is given
by

$P_{N,d,\boldsymbol{z}}=\left\\{\left(\left\\{\frac{nz_{1}}{N}\right\\},\ldots,\left\\{\frac{nz_{d}}{N}\right\\}\right)\mid
0\leq n\leq N-1\right\\}\subset[0,1)^{d},$

where $\\{x\\}:=x-\lfloor x\rfloor$ denotes the fractional part of a real
number $x$. The vector $\boldsymbol{z}$ is called a generating vector.

Evidently we need to construct a good generating vector $\boldsymbol{z}$. The
standard (fast) CBC algorithm is a greedy algorithm which successively
searches for one component $z_{j}$ from the set $\\{1,\ldots,N-1\\}$ at a time
while keeping previous components $z_{1},\ldots,z_{j-1}$ unchanged. To be
precise, the CBC algorithm using the squared worst-case error as a criterion
proceeds as follows:

  1. 1.

Set $z_{j}=1$.

  2. 2.

For $j=2,3,\ldots,d$, choose

$z_{j}=\arg\min_{z\in\\{1,\ldots,N-1\\}}(e^{\mathrm{wor}}(H_{j,\alpha,\boldsymbol{\gamma}},P_{N,j,(z_{1},\ldots,z_{j-1},z)}))^{2}.$

As already pointed out, error bounds resulting from the CBC algorithm have
been well-established in the literature, such as [20, 10, 15, 2] among many
others.

The concept of the dual lattice and its associated character property,
described below, play a central role in our error analysis.

######  Definition 3.2 (Dual lattice).

Let $N,d$ be positive integers and $\boldsymbol{z}\in\\{1,\ldots,N-1\\}^{d}$.
The dual lattice of the rank-1 lattice point set $P_{N,d,\boldsymbol{z}}$ is
defined by

$P^{\perp}_{N,d,\boldsymbol{z}}=\left\\{\boldsymbol{h}\in\mathbb{Z}^{d}\mid\boldsymbol{h}\cdot\boldsymbol{z}\equiv
0\pmod{N}\right\\}\subset\mathbb{Z}^{d}.$

######  Lemma 3.3 (Character property).

Let $N,d$ be positive integers and $\boldsymbol{z}\in\\{1,\ldots,N-1\\}^{d}$.
With the notation above, it holds for $\boldsymbol{h}\in\mathbb{Z}^{d}$ that

$\frac{1}{N}\sum_{\boldsymbol{x}\in P_{N,d,\boldsymbol{z}}}\exp(2\pi
i\boldsymbol{h}\cdot\boldsymbol{x})=\begin{cases}1&\text{if $\boldsymbol{h}\in
P^{\perp}_{N,d,\boldsymbol{z}}$,}\\\ 0&\text{otherwise.}\end{cases}$

Now our concatenated point set is given as follows:

$P_{N,s,d,\boldsymbol{z}}=\left\\{\left(\left\\{\frac{nz_{1}}{N}\right\\},\ldots,\left\\{\frac{nz_{d}}{N}\right\\},x_{n,d+1},\ldots,x_{n,s}\right)\mid
0\leq n\leq N-1\right\\}\subset[0,1)^{s},$

where $x_{n,j},\,n=0,\ldots,N-1,\,j=d+1,\ldots,s$ are sampled independently
from the uniform distribution on $[0,1)$.

##  4 An error bound

In what follows, $\mathbb{E}$ denotes the expected value with respect to the
probability measure, under which the random variables
$x_{n,j},\,n=0,\ldots,N-1,\,j=d+1,\ldots,s$ are independently uniformly
distributed on $[0,1)$. Moreover, for two integers $a,b$ with $a\leq b$, we
write $a:b=\\{a,a+1,\ldots,b\\}$. First we prove the following lemma.

######  Lemma 4.1.

Let $N,s,d$ be positive integers with $d<s$ and $\alpha>1/2$ be a real number.
The mean squared worst-case error of an equally-weighted quadrature rule using
a point set $P_{N,s,d,\boldsymbol{z}}$ for the weighted Korobov space
$H_{s,\alpha,\boldsymbol{\gamma}}$ is given by

$\displaystyle\mathbb{E}\left[(e^{\mathrm{wor}}(H_{s,\alpha,\boldsymbol{\gamma}},P_{N,s,d,\boldsymbol{z}}))^{2}\right]=\frac{1}{N}\sum_{\begin{subarray}{c}\boldsymbol{h}\in\mathbb{Z}^{s}\setminus\\{\boldsymbol{0}\\}\\\
\boldsymbol{h}_{d+1:s}\neq\boldsymbol{0}\end{subarray}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{h})+\sum_{\boldsymbol{h}_{1:d}\in
P^{\perp}_{N,d,\boldsymbol{z}}\setminus\\{\boldsymbol{0}\\}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{h}_{1:d},\boldsymbol{0}).$
(2)

###### Proof.

Using (1), the linearity of expectation and Lemma 3.3, we have

$\displaystyle\mathbb{E}\left[(e^{\mathrm{wor}}(H_{s,\alpha,\boldsymbol{\gamma}},P_{N,s,d,\boldsymbol{z}}))^{2}\right]$
$\displaystyle=\mathbb{E}\left[-1+\frac{1}{N^{2}}\sum_{\boldsymbol{x},\boldsymbol{y}\in
P_{N,s,d,\boldsymbol{z}}}K_{s,\alpha,\boldsymbol{\gamma}}(\boldsymbol{x},\boldsymbol{y})\right]$
$\displaystyle=\mathbb{E}\left[\frac{1}{N^{2}}\sum_{\boldsymbol{x},\boldsymbol{y}\in
P_{N,s,d,\boldsymbol{z}}}\sum_{\boldsymbol{h}\in\mathbb{Z}^{s}\setminus\\{\boldsymbol{0}\\}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{h})\exp\left(2\pi
i\boldsymbol{h}\cdot(\boldsymbol{x}-\boldsymbol{y})\right)\right]$
$\displaystyle=\sum_{\boldsymbol{h}\in\mathbb{Z}^{s}\setminus\\{\boldsymbol{0}\\}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{h})\cdot\mathbb{E}\left[\frac{1}{N^{2}}\sum_{\boldsymbol{x},\boldsymbol{y}\in
P_{N,s,d,\boldsymbol{z}}}\exp\left(2\pi
i\boldsymbol{h}\cdot(\boldsymbol{x}-\boldsymbol{y})\right)\right]$
$\displaystyle=\sum_{\boldsymbol{h}\in\mathbb{Z}^{s}\setminus\\{\boldsymbol{0}\\}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{h})\cdot\mathbb{E}\left[\frac{1}{N^{2}}\sum_{\boldsymbol{x}\in
P_{N,s,d,\boldsymbol{z}}}\exp\left(2\pi
i\boldsymbol{h}\cdot\boldsymbol{0}\right)+\frac{1}{N^{2}}\sum_{\begin{subarray}{c}\boldsymbol{x},\boldsymbol{y}\in
P_{N,s,d,\boldsymbol{z}}\\\
\boldsymbol{x}\neq\boldsymbol{y}\end{subarray}}\exp\left(2\pi
i\boldsymbol{h}\cdot(\boldsymbol{x}-\boldsymbol{y})\right)\right]$
$\displaystyle=\sum_{\boldsymbol{h}\in\mathbb{Z}^{s}\setminus\\{\boldsymbol{0}\\}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{h})\cdot\left[\frac{1}{N}+\frac{1}{N^{2}}\sum_{\begin{subarray}{c}\boldsymbol{x},\boldsymbol{y}\in
P_{N,d,\boldsymbol{z}}\\\
\boldsymbol{x}\neq\boldsymbol{y}\end{subarray}}\exp\left(2\pi
i\boldsymbol{h}_{1:d}\cdot(\boldsymbol{x}-\boldsymbol{y})\right)\prod_{j=d+1}^{s}\int_{0}^{1}\int_{0}^{1}\exp\left(2\pi
ih_{j}\cdot(x_{j}-y_{j})\right)\,\mathrm{d}x_{j}\,\mathrm{d}y_{j}\right]$
$\displaystyle=\sum_{\boldsymbol{h}\in\mathbb{Z}^{s}\setminus\\{\boldsymbol{0}\\}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{h})\cdot\left[\frac{1}{N}+\frac{1}{N^{2}}\sum_{\begin{subarray}{c}\boldsymbol{x},\boldsymbol{y}\in
P_{N,d,\boldsymbol{z}}\\\
\boldsymbol{x}\neq\boldsymbol{y}\end{subarray}}\exp\left(2\pi
i\boldsymbol{h}_{1:d}\cdot(\boldsymbol{x}-\boldsymbol{y})\right)\prod_{j=d+1}^{s}1_{h_{j}=0}\right]$
$\displaystyle=\frac{1}{N}\sum_{\boldsymbol{h}\in\mathbb{Z}^{s}\setminus\\{\boldsymbol{0}\\}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{h})+\frac{1}{N^{2}}\sum_{\begin{subarray}{c}\boldsymbol{h}\in\mathbb{Z}^{s}\setminus\\{\boldsymbol{0}\\}\\\
\boldsymbol{h}_{d+1:s}=\boldsymbol{0}\end{subarray}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{h})\sum_{\begin{subarray}{c}\boldsymbol{x},\boldsymbol{y}\in
P_{N,d,\boldsymbol{z}}\\\
\boldsymbol{x}\neq\boldsymbol{y}\end{subarray}}\exp\left(2\pi
i\boldsymbol{h}_{1:d}\cdot(\boldsymbol{x}-\boldsymbol{y})\right)$
$\displaystyle=\frac{1}{N}\sum_{\begin{subarray}{c}\boldsymbol{h}\in\mathbb{Z}^{s}\setminus\\{\boldsymbol{0}\\}\\\
\boldsymbol{h}_{d+1:s}\neq\boldsymbol{0}\end{subarray}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{h})+\frac{1}{N^{2}}\sum_{\begin{subarray}{c}\boldsymbol{h}\in\mathbb{Z}^{s}\setminus\\{\boldsymbol{0}\\}\\\
\boldsymbol{h}_{d+1:s}=\boldsymbol{0}\end{subarray}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{h})\sum_{\boldsymbol{x},\boldsymbol{y}\in
P_{N,d,\boldsymbol{z}}}\exp\left(2\pi
i\boldsymbol{h}_{1:d}\cdot(\boldsymbol{x}-\boldsymbol{y})\right)$
$\displaystyle=\frac{1}{N}\sum_{\begin{subarray}{c}\boldsymbol{h}\in\mathbb{Z}^{s}\setminus\\{\boldsymbol{0}\\}\\\
\boldsymbol{h}_{d+1:s}\neq\boldsymbol{0}\end{subarray}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{h})+\sum_{\begin{subarray}{c}\boldsymbol{h}\in\mathbb{Z}^{s}\setminus\\{\boldsymbol{0}\\}\\\
\boldsymbol{h}_{1:d}\in P^{\perp}_{N,d,\boldsymbol{z}}\\\
\boldsymbol{h}_{d+1:s}=\boldsymbol{0}\end{subarray}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{h})$
$\displaystyle=\frac{1}{N}\sum_{\begin{subarray}{c}\boldsymbol{h}\in\mathbb{Z}^{s}\setminus\\{\boldsymbol{0}\\}\\\
\boldsymbol{h}_{d+1:s}\neq\boldsymbol{0}\end{subarray}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{h})+\sum_{\boldsymbol{h}_{1:d}\in
P^{\perp}_{N,d,\boldsymbol{z}}\setminus\\{\boldsymbol{0}\\}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{h}_{1:d},\boldsymbol{0}).$

Thus we are done. ∎

Now we show the main result of this paper.

######  Theorem 4.2.

Let $N,s,d$ be positive integers with $d<s$ and $\alpha>1/2$ be a real number.
Assume that there exists $0<\lambda^{*}<1$ such that

$\displaystyle\sum_{j=1}^{\infty}\gamma_{j}^{\lambda^{*}}<\infty.$ (3)

Let $\boldsymbol{z}\in\\{1,\ldots,N-1\\}^{d}$ be constructed by the (fast)
component-by-component algorithm with the quality criterion
$(e^{\mathrm{wor}}(H_{d,\alpha,\boldsymbol{\gamma}},P_{N,d,\boldsymbol{\gamma}}))^{2}$.
Then the mean squared worst-case error of an equally-weighted quadrature rule
using a point set $P_{N,s,d,\boldsymbol{z}}$ for the weighted Korobov space
$H_{s,\alpha,\boldsymbol{\gamma}}$ is bounded above by

$\displaystyle\mathbb{E}\left[(e^{\mathrm{wor}}(H_{s,\alpha,\boldsymbol{\gamma}},P_{N,s,d,\boldsymbol{z}}))^{2}\right]$
$\displaystyle\leq\left[\frac{1}{\varphi(N)}\left(-1+\exp\left(2C^{\lambda}\zeta(2\alpha\lambda)\zeta(\lambda/\lambda^{*})\right)\right)\right]^{1/\lambda}$
$\displaystyle\quad+\frac{\exp\left(2C\zeta(2\alpha)\zeta(1/\lambda^{*})\right)}{Nd^{1/\lambda^{*}-1}}\left[-1+\exp\left(\frac{2C\zeta(2\alpha)}{1/\lambda^{*}-1}\right)\right],$
(4)

for any $\lambda\in(\max(\lambda^{*},1/(2\alpha)),1]$, where $C>0$ is a
constant depending only on $\boldsymbol{\gamma}$, and $\varphi(\cdot)$ and
$\zeta(\cdot)$ denote the Euler totient function and the Riemann zeta
function, respectively,

######  Remark 4.3.

It is shown in [19, Theorem 15] that

$\frac{1}{\varphi(N)}\leq\frac{1}{N}\left(e^{c}\log\log
N+\frac{2.50637}{\log\log N}\right)$

for any $N\geq 3$ with $c=0.577\ldots$ being the Euler’s constant.

###### Proof of Theorem 4.2.

It suffices from Lemma 4.1 to give upper bounds on the two terms on the right-
hand side of (2), respectively. The second term on the right-hand side of (2)
is exactly the squared worst-case error for a rank-1 lattice rule in
$H_{d,\alpha,\boldsymbol{\gamma}}$ (not in $H_{s,\alpha,\boldsymbol{\gamma}}$,
of course). By choosing a generating vector $\boldsymbol{z}$ by the standard
(fast) component-by-component algorithm, it can be bounded above by

$\sum_{\boldsymbol{h}_{1:d}\in
P^{\perp}_{N,d,\boldsymbol{z}}\setminus\\{\boldsymbol{0}\\}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{h}_{1:d},\boldsymbol{0})\leq\left[\frac{1}{\varphi(N)}\left(-1+\prod_{j=1}^{d}\left[1+2\gamma_{j}^{\lambda}\zeta(2\alpha\lambda)\right]\right)\right]^{1/\lambda},$

for any $\lambda\in(1/(2\alpha),1]$, see [5, Theorem 5.12]. Given the
summability condition (3), we can assume the existence of $C>0$, which depends
only on $\boldsymbol{\gamma}$, such that

$\displaystyle\gamma_{j}\leq Cj^{-1/\lambda^{*}}$ (5)

for all $j\in\mathbb{N}$. By using the inequality $\log(1+x)\leq x$ for any
$x>0$, it holds for any $\lambda\in(\max(\lambda^{*},1/(2\alpha)),1]$ that

$\displaystyle\sum_{\boldsymbol{h}_{1:d}\in
P^{\perp}_{N,d,\boldsymbol{z}}\setminus\\{\boldsymbol{0}\\}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{h}_{1:d},\boldsymbol{0})$
$\displaystyle\leq\left[\frac{1}{\varphi(N)}\left(-1+\exp\left(\sum_{j=1}^{d}\log(1+2\gamma_{j}^{\lambda}\zeta(2\alpha\lambda))\right)\right)\right]^{1/\lambda}$
$\displaystyle\leq\left[\frac{1}{\varphi(N)}\left(-1+\exp\left(2\zeta(2\alpha\lambda)\sum_{j=1}^{d}\gamma_{j}^{\lambda}\right)\right)\right]^{1/\lambda}$
$\displaystyle\leq\left[\frac{1}{\varphi(N)}\left(-1+\exp\left(2C^{\lambda}\zeta(2\alpha\lambda)\sum_{j=1}^{\infty}j^{-\lambda/\lambda^{*}}\right)\right)\right]^{1/\lambda}$
$\displaystyle=\left[\frac{1}{\varphi(N)}\left(-1+\exp\left(2C^{\lambda}\zeta(2\alpha\lambda)\zeta(\lambda/\lambda^{*})\right)\right)\right]^{1/\lambda}.$

Let us look at the first term on the right-hand side of (2). We have

$\displaystyle\frac{1}{N}\sum_{\begin{subarray}{c}\boldsymbol{h}\in\mathbb{Z}^{s}\setminus\\{\boldsymbol{0}\\}\\\
\boldsymbol{h}_{d+1:s}\neq\boldsymbol{0}\end{subarray}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{h})$
$\displaystyle=\frac{1}{N}\sum_{\begin{subarray}{c}\boldsymbol{h}_{1:d}\in\mathbb{Z}^{d}\\\
\boldsymbol{h}_{d+1:s}\in\mathbb{Z}^{s-d}\setminus\\{\boldsymbol{0}\\}\end{subarray}}\prod_{\begin{subarray}{c}j=1\\\
h_{j}\neq 0\end{subarray}}^{s}\frac{\gamma_{j}}{|h_{j}|^{2\alpha}}$
$\displaystyle=\frac{1}{N}\left(\sum_{\boldsymbol{h}_{1:d}\in\mathbb{Z}^{d}}\prod_{\begin{subarray}{c}j=1\\\
h_{j}\neq
0\end{subarray}}^{d}\frac{\gamma_{j}}{|h_{j}|^{2\alpha}}\right)\left(\sum_{\boldsymbol{h}_{d+1:s}\in\mathbb{Z}^{s-d}\setminus\\{\boldsymbol{0}\\}}\prod_{\begin{subarray}{c}j=d+1\\\
h_{j}\neq 0\end{subarray}}^{s}\frac{\gamma_{j}}{|h_{j}|^{2\alpha}}\right)$
$\displaystyle=\frac{1}{N}\prod_{j=1}^{d}\left[1+2\gamma_{j}\zeta(2\alpha)\right]\left(-1+\prod_{j=d+1}^{s}\left[1+2\gamma_{j}\zeta(2\alpha)\right]\right).$

Then, by taking account of the decay (5) and using the inequality
$\log(1+x)\leq x$ again, we have

$\displaystyle\prod_{j=1}^{d}\left[1+2\gamma_{j}\zeta(2\alpha)\right]$
$\displaystyle=\exp\left(\sum_{j=1}^{d}\log(1+2\gamma_{j}\zeta(2\alpha))\right)$
$\displaystyle\leq\exp\left(2\zeta(2\alpha)\sum_{j=1}^{d}\gamma_{j}\right)$
$\displaystyle\leq\exp\left(2C\zeta(2\alpha)\sum_{j=1}^{\infty}j^{-1/\lambda^{*}}\right)$
$\displaystyle=\exp\left(2C\zeta(2\alpha)\zeta(1/\lambda^{*})\right),$

and

$\displaystyle-1+\prod_{j=d+1}^{s}\left[1+2\gamma_{j}\zeta(2\alpha)\right]$
$\displaystyle=-1+\exp\left(\sum_{j=d+1}^{s}\log(1+2\gamma_{j}\zeta(2\alpha))\right)$
$\displaystyle\leq-1+\exp\left(2\zeta(2\alpha)\sum_{j=d+1}^{s}\gamma_{j}\right)$
$\displaystyle\leq-1+\exp\left(2C\zeta(2\alpha)\sum_{j=d+1}^{\infty}j^{-1/\lambda^{*}}\right)$
$\displaystyle\leq-1+\exp\left(2C\zeta(2\alpha)\int_{d}^{\infty}x^{-1/\lambda^{*}}\,\mathrm{d}x\right)$
$\displaystyle=-1+\exp\left(\frac{2C\zeta(2\alpha)}{1/\lambda^{*}-1}d^{-1/\lambda^{*}+1}\right)$
$\displaystyle\leq\left[-1+\exp\left(\frac{2C\zeta(2\alpha)}{1/\lambda^{*}-1}\right)\right]d^{-1/\lambda^{*}+1}.$

Here the last inequality is obtained as follows: given that the function
$f(x)=-1+\exp(x)$ is convex, Jensen’s inequality leads to

$\displaystyle-1+\exp\left(\frac{2C\zeta(2\alpha)}{1/\lambda^{*}-1}d^{-1/\lambda^{*}+1}\right)$
$\displaystyle=f\left(\frac{2C\zeta(2\alpha)}{1/\lambda^{*}-1}d^{-1/\lambda^{*}+1}\right)$
$\displaystyle\leq(1-d^{-1/\lambda^{*}+1})f(0)+d^{-1/\lambda^{*}+1}f\left(\frac{2C\zeta(2\alpha)}{1/\lambda^{*}-1}\right)$
$\displaystyle=\left[-1+\exp\left(\frac{2C\zeta(2\alpha)}{1/\lambda^{*}-1}\right)\right]d^{-1/\lambda^{*}+1}.$

Altogether the first term on the right-hand side of (2) is bounded by

$\frac{1}{N}\sum_{\begin{subarray}{c}\boldsymbol{h}\in\mathbb{Z}^{s}\setminus\\{\boldsymbol{0}\\}\\\
\boldsymbol{h}_{d+1:s}\neq\boldsymbol{0}\end{subarray}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{h})\leq\frac{\exp\left(2C\zeta(2\alpha)\zeta(1/\lambda^{*})\right)}{Nd^{1/\lambda^{*}-1}}\left[-1+\exp\left(\frac{2C\zeta(2\alpha)}{1/\lambda^{*}-1}\right)\right].$

Since both the terms on the right-hand side of (2) are now bounded above as
desired, we are done. ∎

It is important that our obtained upper bound on the mean squared worst-case
error is independent of the problem dimension $s$. Assuming $s\gg N,d$ or
$s\to\infty$, let us discuss a choice of $d$.

  1. 1.

If $\lambda^{*}\leq 1/(2\alpha)$, the bound (4) holds for any
$\lambda\in(1/(2\alpha),1]$. Since it is known that we cannot achieve a
convergence rate better than $O(N^{-\alpha})$ for the worst-case error, see
[13, Theorem 4.22] for a lower bound on the worst-case error, the order of the
first term of (4) is almost optimal, as evident from Remark 4.3. Here we
recall that, although we consider the average of the squared worst-case error
in this note, the square of a lower bound on the deterministic worst-case
error still applies; otherwise a contradiction occurs. Thus, in order for the
second term not to be dominant, it suffices to choose

$d\propto N^{(2\alpha-1)/(1/\lambda^{*}-1)}<N.$

For such a choice, the required construction cost for the first coordinates is
of order

$dN\log N\propto N^{1+(2\alpha-1)/(1/\lambda^{*}-1)}\log N.$

A substantial cost saving can be expected for smaller $\lambda^{*}$, as
compared to the standard algorithm of order $sN\log N$.

  2. 2.

If $1/(2\alpha)<\lambda^{*}\leq 1$, on the other hand, the bound (4) holds
only for the range $\lambda\in(\lambda^{*},1]$, which does not lead to an
optimal order error bound. Nevertheless, in order for the second term not to
be dominant, it suffices to choose

$d\propto N$

regardless of the value of $\lambda^{*}$, resulting in the necessary
construction cost of order

$dN\log N\propto N^{2}\log N.$

This might be useful when $s$ is extremely large.

###  4.1 Comparison with the reduced CBC algorithm

Let us consider rank-1 lattice rules constructed the reduced CBC algorithm
introduced in [4], in which $\alpha$ should be replaced by $2\alpha$ to be
consistent with this note. The worst-case error is bounded independently of
the dimension $s$ if there exist $0<\lambda^{*}<1$ and $0\leq w_{1}\leq
w_{2}\leq\cdots$ with $w_{j}\in\mathbb{N}\cup\\{0\\}$ such that

$\displaystyle\sum_{j=1}^{\infty}\gamma_{j}^{\lambda^{*}}b^{w_{j}}<\infty$ (6)

holds for a fixed prime $b$. The squared worst-case error decays with order
$N^{-1/\lambda}$ for any $\lambda\in(\max(\lambda^{*},1/(2\alpha)),1]$. Note
that such results for the reduced CBC algorithm are deterministic on the
contrary to our concatenated rule. Also, because of the shift-invariance of
the kernel $K_{s,\alpha,\boldsymbol{\gamma}}$, the same result holds for the
mean squared worst-case error of a randomly shifted rank-1 lattice rule. When
the number of points is given by $N=b^{m}$ for some $m\in\mathbb{N}$, the
required construction cost for the reduced CBC algorithm, improved in [6], is
of order

$\displaystyle\sum_{j=1}^{\min(s,s^{*})}(m-w_{j})b^{m-w_{j}},$ (7)

where $s^{*}$ is defined to be the largest integer such that $w_{s^{*}}<m$. In
the following, assuming that the weights are given by
$\gamma_{j}=j^{-2c\alpha}$ for some $c>1/(2\alpha)$, we show that this
construction cost is in general not comparable to that for our concatenated
rule.

#### Slowly decaying weights

First let us consider the case $1/(2\alpha)<c<1$. The condition (3) holds for
$\lambda^{*}=1/(2c\alpha)+\epsilon$ with arbitrarily small $\epsilon>0$. As
discussed above, the necessary construction cost for our concatenated rule is
of $O(N^{2}\log N)$. Regarding the reduced CBC algorithm, the condition (6)
holds trivially for $\lambda^{*}=1/(2c\alpha)+\epsilon$ if
$w_{1}=w_{2}=\cdots=0$. Then the reduced CBC algorithm coincides with the
standard CBC algorithm, requiring the necessary construction cost of $O(sN\log
N)$. By setting $w_{j}=\lfloor\beta\log_{b}j\rfloor$ for some
$0<\beta<2c\alpha-1$, we can set $s^{*}\leq N^{1/\beta}$ and it follows from
(7) that the order of construction cost becomes at most

$\displaystyle\sum_{j=1}^{\min(s,s^{*})}(m-w_{j})b^{m-w_{j}}\leq
mb^{m}\sum_{j=1}^{N^{1/\beta}}\frac{b}{j^{\beta}}\leq\begin{cases}b\zeta(\beta)N\log_{b}N&\text{for
$\beta>1$,}\\\ (b\log b)N(\log_{b}Ne)^{2}&\text{for $\beta=1$,}\\\
bN^{1/\beta}\log_{b}N&\text{otherwise.}\end{cases}$ (8)

Thus the dividing case in terms of construction cost is $\beta=1/2$. However,
for any $\beta>0$, the condition (6) holds for
$\lambda^{*}=(\beta+1)/(2c\alpha)+\epsilon$, which deteriorates the decay of
the worst-case error.

#### Fast decaying weights

Let us consider the case $c>1$ next. The condition (3) holds for
$\lambda^{*}=1/(2c\alpha-\epsilon)$ with arbitrarily small $\epsilon>0$. The
necessary construction cost for our concatenated rule is of order $N^{a}\log
N$ with

$a=1+\frac{2\alpha-1}{2c\alpha-\epsilon-1}.$

Regarding the reduced CBC algorithm, by setting
$w_{j}=\lfloor\beta\log_{b}j\rfloor$ for some $0<\beta<2c\alpha-1$, the
condition (6) holds for $\lambda^{*}=(\beta+1)/(2c\alpha)+\epsilon$. Looking
at (8), the dividing case in terms of construction cost is

$\beta=\frac{1}{a}=\left(1+\frac{2\alpha-1}{2c\alpha-\epsilon-1}\right)^{-1}<1.$

However, in this case, the reduced CBC algorithm can achieve almost the
optimal order of the worst-case error as long as
$\beta<c(1-2\alpha\epsilon)-1$. Thus, if $c$ is large enough, we can set
$\beta\geq 1$, providing a smaller construction cost than our concatenated
rule.

###  4.2 Comparison with truncation algorithm

Finally we compare the result for our concatenated rule with that for a
truncation algorithm inspired by [9]. Let us consider the point set

$P_{N,s,d,\boldsymbol{z},\boldsymbol{c}}^{\mathrm{tr}}=\left\\{\left(\left\\{\frac{nz_{1}}{N}\right\\},\ldots,\left\\{\frac{nz_{d}}{N}\right\\},c_{d+1},\ldots,c_{s}\right)\mid
0\leq n\leq N-1\right\\},$

for an anchor vector $\boldsymbol{c}=(c_{d+1},\ldots,c_{s})\in[0,1)^{s-d}$,
which can be fixed or chosen independently uniformly from $[0,1)^{s-d}$. We
obtain the following worst-case error bound for the truncation algorithm.
Although a similar result can be found, for instance, in [9, Theorem 13], we
give a proof for the sake of completeness.

######  Theorem 4.4.

Let $N,s,d$ be positive integers with $d<s$ and $\alpha>1/2$ be a real number.
Assume that there exists $0<\lambda^{*}<1$ such that (3) holds. Let
$\boldsymbol{z}\in\\{1,\ldots,N-1\\}^{d}$ be constructed by the (fast)
component-by-component algorithm with the quality criterion
$(e^{\mathrm{wor}}(H_{d,\alpha,\boldsymbol{\gamma}},P_{N,d,\boldsymbol{\gamma}}))^{2}$.
Then, for any $\boldsymbol{c}\in[0,1)^{s-d}$, the squared worst-case error of
an equally-weighted quadrature rule using a point set
$P_{N,s,d,\boldsymbol{z},\boldsymbol{c}}^{\mathrm{tr}}$ for the weighted
Korobov space $H_{s,\alpha,\boldsymbol{\gamma}}$ is bounded above by

$\displaystyle(e^{\mathrm{wor}}(H_{s,\alpha,\boldsymbol{\gamma}},P_{N,s,d,\boldsymbol{z},\boldsymbol{c}}^{\mathrm{tr}}))^{2}$
$\displaystyle\leq\left[\frac{1}{\varphi(N)}\left(-1+\exp\left(2C^{\lambda}\zeta(2\alpha\lambda)\zeta(\lambda/\lambda^{*})\right)\right)\right]^{1/\lambda}\exp\left(2C\zeta(2\alpha)\zeta(1/\lambda^{*})\right)$
$\displaystyle\quad+\frac{1}{d^{1/\lambda^{*}-1}}\left[-1+\exp\left(\frac{2C\zeta(2\alpha)}{1/\lambda^{*}-1}\right)\right],$
(9)

for any $\lambda\in(\max(\lambda^{*},1/(2\alpha)),1]$, where $C>0$ is a
constant depending only on $\boldsymbol{\gamma}$.

###### Proof.

Because of the shift-invariance of the kernel
$K_{s,\alpha,\boldsymbol{\gamma}}$, it suffices to prove the statement for the
case $\boldsymbol{c}=\boldsymbol{0}$. Following an argument similar to the
proof of Lemma 4.1, we have

$\displaystyle(e^{\mathrm{wor}}(H_{s,\alpha,\boldsymbol{\gamma}},P_{N,s,d,\boldsymbol{z},\boldsymbol{0}}^{\mathrm{tr}}))^{2}$
$\displaystyle=\frac{1}{N^{2}}\sum_{\boldsymbol{x},\boldsymbol{y}\in
P_{N,s,d,\boldsymbol{z},\boldsymbol{0}}^{\mathrm{tr}}}\sum_{\boldsymbol{h}\in\mathbb{Z}^{s}\setminus\\{\boldsymbol{0}\\}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{h})\exp\left(2\pi
i\boldsymbol{h}\cdot(\boldsymbol{x}-\boldsymbol{y})\right)$
$\displaystyle=\sum_{\boldsymbol{h}\in\mathbb{Z}^{s}\setminus\\{\boldsymbol{0}\\}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{h})\frac{1}{N^{2}}\sum_{\boldsymbol{x},\boldsymbol{y}\in
P_{N,d,\boldsymbol{z}}}\exp\left(2\pi
i\boldsymbol{h}_{1:d}\cdot(\boldsymbol{x}-\boldsymbol{y})\right)$
$\displaystyle=\sum_{\begin{subarray}{c}\boldsymbol{h}\in\mathbb{Z}^{s}\setminus\\{\boldsymbol{0}\\}\\\
\boldsymbol{h}_{1:d}\in
P^{\perp}_{N,d,\boldsymbol{z}}\end{subarray}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{h})$
$\displaystyle=\sum_{\begin{subarray}{c}\boldsymbol{h}_{1:d}\in
P^{\perp}_{N,d,\boldsymbol{z}}\setminus\\{\boldsymbol{0}\\}\\\
\boldsymbol{h}_{d+1:s}\in\mathbb{Z}^{s-d}\end{subarray}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{h})+\sum_{\boldsymbol{h}_{d+1:s}\in\mathbb{Z}^{s-d}\setminus\\{\boldsymbol{0}\\}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{0},\boldsymbol{h}_{d+1:s}).$

As in the proof of Theorem 4.2, with $C$ being a positive constant depending
only on $\boldsymbol{\gamma}$, the first and second terms are bounded above by

$\displaystyle\sum_{\begin{subarray}{c}\boldsymbol{h}_{1:d}\in
P^{\perp}_{N,d,\boldsymbol{z}}\setminus\\{\boldsymbol{0}\\}\\\
\boldsymbol{h}_{d+1:s}\in\mathbb{Z}^{s-d}\end{subarray}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{h})$
$\displaystyle=\sum_{\boldsymbol{h}_{1:d}\in
P^{\perp}_{N,d,\boldsymbol{z}}\setminus\\{\boldsymbol{0}\\}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{h}_{1:d},\boldsymbol{0})\sum_{\boldsymbol{h}_{d+1:s}\in\mathbb{Z}^{s-d}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{0},\boldsymbol{h}_{d+1:s})$
$\displaystyle\leq\left[\frac{1}{\varphi(N)}\left(-1+\exp\left(2C^{\lambda}\zeta(2\alpha\lambda)\zeta(\lambda/\lambda^{*})\right)\right)\right]^{1/\lambda}\prod_{j=d+1}^{s}\left[1+2\gamma_{j}\zeta(2\alpha)\right]$
$\displaystyle\leq\left[\frac{1}{\varphi(N)}\left(-1+\exp\left(2C^{\lambda}\zeta(2\alpha\lambda)\zeta(\lambda/\lambda^{*})\right)\right)\right]^{1/\lambda}\exp\left(2C\zeta(2\alpha)\zeta(1/\lambda^{*})\right),$

for any $\lambda\in(\max(\lambda^{*},1/(2\alpha)),1]$ and

$\displaystyle\sum_{\boldsymbol{h}_{d+1:s}\in\mathbb{Z}^{s-d}\setminus\\{\boldsymbol{0}\\}}r_{\alpha,\boldsymbol{\gamma}}(\boldsymbol{0},\boldsymbol{h}_{d+1:s})$
$\displaystyle=-1+\prod_{j=d+1}^{s}\left[1+2\gamma_{j}\zeta(2\alpha)\right]$
$\displaystyle\leq\left[-1+\exp\left(\frac{2C\zeta(2\alpha)}{1/\lambda^{*}-1}\right)\right]d^{-1/\lambda^{*}+1},$

respectively, which completes the proof. ∎

In contrast to the bound (4) for our concatenated rule, the one (9) we show
here is a deterministic worst-case error bound, which is a clear advantage of
the truncation algorithm. On the other hand, in order for the second term of
(9) not to be dominant, the truncation dimension $d$ should be larger than
that for our concatenated rule. This is simply because the $N$ factor is
missing in the denominator of the second term as compared to (4). In fact, a
sufficient choice for the truncation algorithm is given by

$d\propto\begin{cases}N^{2\alpha/(1/\lambda^{*}-1)}&\text{if
$0<\lambda^{*}\leq 1/(2\alpha)$,}\\\ N^{1/(1-\lambda^{*})}&\text{if
$1/(2\alpha)<\lambda^{*}<1$,}\end{cases}$

which can make a significant difference of the necessary construction cost for
the first $d$ coordinates. This is a point where randomization helps for our
high-dimensional integration problems.

## Acknowledgements

The author would like to thank Kei Ishikawa (ETH Zürich) for bringing the
problem to his attention. He is very grateful to Michael Gnewuch (Osnabrück)
for pointing out some relevant literature and also an error contained in the
first version of this manuscript. He also appreciates the comments and
suggestions from the anonymous reviewers. The work of T.G. is supported by
JSPS KAKENHI Grant Number 20K03744.

## References

  * [1] C. Aistleitner and M. Hofer.  Probabilistic error bounds for the discrepancy of mixed sequences.  Monte Carlo Methods and Applications, 18(2):181–200, 2012.
  * [2] J. Dick and T. Goda.  Stability of lattice rules and polynomial lattice rules constructed by the component-by-component algorithm.  Journal of Computational and Applied Mathematics, 382:113062, 2021\.
  * [3] J. Dick and P. Kritzer.  On a projection-corrected component-by-component construction.  Journal of Complexity, 32(1):74–80, 2016.
  * [4] J. Dick, P. Kritzer, G. Leobacher, and F. Pillichshammer.  A reduced fast component-by-component construction of lattice points for integration in weighted spaces with fast decreasing weights.  Journal of Computational and Applied Mathematics, 276:1–15, 2015\.
  * [5] J. Dick, F. Y. Kuo, and I. H. Sloan.  High-dimensional integration: the quasi-Monte Carlo way.  Acta Numerica, 22:133–288, 2013.
  * [6] A. Ebert, P. Kritzer, and D. Nuyens.  Constructing QMC finite element methods for elliptic PDEs with random coefficients by a reduced CBC construction.  In B. Tuffin and P. L’Ecuyer, editors, Monte Carlo and Quasi-Monte Carlo Methods, pages 183–205. Springer Proceedings in Mathematics & Statistics 324, Springer, 2020.
  * [7] M. Gnewuch.  On probabilistic results for the discrepancy of a hybrid-Monte Carlo sequence.  Journal of Complexity, 25(4):312–317, 2009.
  * [8] Y. Kazashi.  Quasi–Monte Carlo integration with product weights for elliptic PDEs with log-normal coefficients.  IMA Journal of Numerical Analysis, 39(3):1563–1593, 2019.
  * [9] P. Kritzer, F. Pillichshammer, and G. W. Wasilkowski.  Very low truncation dimension for high dimensional integration under modest error demand.  Journal of Complexity, 36:63–85, 2016.
  * [10] F. Y. Kuo.  Component-by-component constructions achieve the optimal rate of convergence for multivariate integration in weighted Korobov and Sobolev spaces.  Journal of Complexity, 19(3):301–320, 2003.
  * [11] F. Y. Kuo and D. Nuyens.  Application of quasi-Monte Carlo methods to elliptic PDEs with random diffusion coefficients: a survey of analysis and implementation.  Foundations of Computational Mathematics, 16(6):1631–1696, 2016\.
  * [12] F. Y. Kuo, C. Schwab, and I. H. Sloan.  Quasi-Monte Carlo finite element methods for a class of elliptic partial differential equations with random coefficients.  SIAM Journal on Numerical Analysis, 50(6):3351–3374, 2012.
  * [13] G. Leobacher and F. Pillichshammer.  Introduction to Quasi-Monte Carlo Integration and Applications.  Springer International Publishing, 2014.
  * [14] H. Niederreiter.  On the discrepancy of some hybrid sequences.  Acta Arithmetica, 138(4):373–398, 2009.
  * [15] D. Nuyens and R. Cools.  Fast algorithms for component-by-component construction of rank-1 lattice rules in shift-invariant reproducing kernel Hilbert spaces.  Mathematics of Computation, 75(254):903–920, 2006.
  * [16] G. Ökten.  A probabilistic result on the discrepancy of a hybrid-Monte Carlo sequence and applications.  Monte Carlo Methods and Applications, 2(4):250–270, 1996.
  * [17] A. B. Owen.  Lattice sampling revisited: Monte Carlo variance of means over randomized orthogonal arrays.  Annals of Statistics, 22(2):930–945, 1994.
  * [18] A. B. Owen.  Latin supercube sampling for very high-dimensional simulations.  ACM Transactions on Modeling and Computer Simulation, 8(1):71–102, 1998.
  * [19] J. B. Rosser and L. Schoenfeld.  Approximate formulas for some functions of prime numbers.  Illinois Journal of Mathematics, 6(1):64–94, 1962.
  * [20] I. H. Sloan and A. V. Reztsov.  Component-by-component construction of good lattice rules.  Mathematics of Computation, 71(237):263–273, 2002.
  * [21] I. H. Sloan and H. Woźniakowski.  When are quasi-Monte Carlo algorithms efficient for high-dimensional integrals?  Journal of Complexity, 14(1):1–33, 1998.
  * [22] J. Spanier.  Quasi-Monte Carlo methods for particle transport problems.  In H. Niederreiter and P. J.-S. Shiue, editors, Monte Carlo and Quasi-Monte Carlo Methods in Scientific Computing, pages 121–148. Lecture Notes in Statistics 106, Springer, 1995.
  * [23] M. Wnuk, M. Gnewuch, and N. Hebbinghaus.  On negatively dependent sampling schemes, variance reduction, and probabilistic upper discrepancy bounds.  In D. Bylik, J. Dick, and F. Pillichshammer, editors, Discrepancy Theory, pages 43–68. Radon Series on Computational and Applied Mathematics 26, De Gruyter, 2020.


# Generating ultra compact boson stars with modified scalar potentials

Sarah Louisa Pitz<EMAIL_ADDRESS>Institut für Theoretische Physik,
Goethe Universität, Max-von-Laue Straße 1, D-60438 Frankfurt am Main, Germany
Jürgen Schaffner-Bielich<EMAIL_ADDRESS>Institut für
Theoretische Physik, Goethe Universität, Max-von-Laue Straße 1, D-60438
Frankfurt am Main, Germany

###### Abstract

The properties of self-interacting boson stars with different scalar
potentials going beyond the commonly used $\phi^{4}$ ansatz are studied. The
scalar potential is extended to different values of the exponent $n$ of the
form $V\propto\phi^{n}$. Two stability mechanism for boson stars are
introduced, the first being a mass term and the second one a vacuum term. We
present analytic scale-invariant expressions for these two classes of EOS. The
resulting properties of the boson star configurations differ considerably from
previous calculations. We find three different categories of mass-radius
relation: the first category resembles the mass-radius curve of self-bound
stars, the second one those of neutron stars and the third one is the well
known constant radius case from the standard $\phi^{4}$ potential. We
demonstrate that the maximal compactness can reach extremely high values going
to the limit of causality $C_{\text{max}}=0.354$ asymptotically for
$n\to\infty$. The maximal compactnesses exceed previously calculated values of
$C_{\text{max}}=0.16$ for the standard $\phi^{4}$-theory and
$C_{\text{max}}=0.21$ for vector-like interactions and is in line with
previous results for solitonic boson stars. Hence, boson stars even described
by a simple modified scalar potential in the form of $V\propto\phi^{n}$ can be
ultra compact black hole mimickers where the photon ring is located outside
the radius of the star.

##  I Introduction

Dark matter plays a crucial role in explaining certain phenomena in cosmology
and astrophysics on large scales, where self-interacting dark matter provides
an explanation for small-scale structure observations such as the core-cusp
and the missing satellites problem [1]. By including self-interactions dark
matter can form compact objects, such as boson stars. Boson stars are self-
gravitating spheres, described by complex scalar fields, see [2] for a review.
These kind of compact objects were discussed for the first time by Wheeler in
1955 [3] with self-gravitating, non-interacting spheres made of bosons, which
he named geons. However, these configurations turned out to be unstable. In
order to obtain stable stars one needs time dependent solutions of the Klein-
Gordon equation [4]. These stable boson stars were usually of microscopic
sizes for the non-interacting case [5, 6]. Boson stars of astrophysical sizes
were found for solitonic boson stars [7, 8, 9, 10], by introducing a self-
interaction potential of the form $V=\lambda\phi^{4}$ [11], and for repulsive
vector-like self-interactions proportional to the density squared [12].
Generic self-interactions have been considered in [13]. The properties of
boson stars with self-interactions have been investigated in detail using a
$\phi^{4}$-potential or vector-like self-interactions [14, 15, 16, 17]. Boson
stars can be built via standard structure formation from the early universe
[18, 19]. So far the LIGO-Virgo collaboration has only detected one event of a
neutron star-neutron star merger, GW170817 [20], confirmed by a gamma-ray
burst and the optical afterglow. A significant amount of the gravitational
wave sources measured by the LIGO-Virgo collaboration is located in the mass
range between the lightest black hole of $5M_{\odot}$ and the heaviest neutron
star [21]. These compact objects could possibly be exotic stars since the most
massive neutron stars measured so far have masses of $2.01\pm 0.04M_{\odot}$
respectively $2.08\pm 0.07M_{\odot}$, constrained by observations of the radio
pulsars PSR J0348+0432 and PSR J0740+6620 [22, 23, 24] as well as $2.35\pm
0.17M_{\odot}$, constrained from optical observations of the black widow
pulsar PSR J0952-0607 [25]. For example, in the gravitational event GW190814
one compact object has a mass of $2.6M_{\odot}$, exceeding all known neutron
star masses considerably [26]. Boson stars can be observed by the emission of
gravitational waves from a merger of two boson stars or one boson star with a
neutron star. Moreover gravitational waves from boson star-boson star merger
can be distinguished from other mergers by their echo in gravitational waves
[27, 28, 29, 30, 31, 32, 33, 34]. Also collisions of neutron stars with boson
stars can send out gravitational waves that are detectable by future
telescopes [35]. Here an unusual deviation by e.g. an extremely small value of
the tidal deformability $\Lambda$ can indicate such an exotic encounter. We
will find below that the maximum compactness $C=GM/R$ for boson stars is
$C=0.354$. Since the tidal deformability goes like $\Lambda\sim k_{2}C^{-5}$
with a Love number of $k_{2}\approx 0.03$ or smaller the value for boson stars
can be as low as of the order of $1$, which is much smaller compared to the
standard values for ordinary neutron stars. For comparison we note that the
limit $\Lambda<720$ has been extracted from the neutron star merger event
GW170817 [36]. Another possibility of gravitational wave sources are ordinary
neutron stars containing a bosonic dark matter fluid [37, 38, 39, 40, 41, 42,
43]. These two-fluid compact objects could contain dark matter cores which are
influencing the macroscopic properties of the star. With the advent of
observing supermassive black holes by radio observations of the Event Horizon
Telescope (EHT) collaboration one can study the accretion disk of these
objects in detail [44, 45]. The precision of these measurements allows to
distinguish between black holes and other compact objects, as e.g. boson
stars. In comparison to black holes boson stars don’t have a shadow since
light does not interact with dark matter and thus passes through the star.
This feature has been used in the analysis of the first pictures from black
hole M87* and SGR A* by radio observations to constrain models of boson stars
as black hole mimickers [46]. In our work we are extending the standard
$\phi^{4}$ self-interaction potential of bosons to a general one of the form
$\phi^{n}$ with arbitrary values of the power $n$. Examples like the three
gluon interaction for $n=3$ and the Wess-Zumino-Witten term in the chiral
Lagrangian to describe the coupling between two kaons and three pions for
$n=5$, motivate the study of other values of $n$. By choosing this generalized
self-interaction potential we find that compact stars with massive bosons can
have different forms of their mass-radius relations, with some of them being
similar to those of neutron stars or self-bound stars. We find curves that are
not constrained to the constant radius case for the $\phi^{4}$-potential and
the one of vector-like self-interactions. We will also show that their
compactness can exceed those of neutron stars and goes asymptotically to the
limit of causality with a compactness of $C=0.354$ for large values of $n$.
These extreme values of compactness are similar to the ones found in a recent
work on solitonic boson stars [47] and lead to new opportunities in the search
of self-interacting boson stars via merging boson stars and thus in the search
of self-interacting dark matter. For a first investigation of ultra compact
solitonic boson star mergers we refer to [34]. The outline of the paper is as
follows: first the theoretical basis from classical field theory for complex
scalar fields and general relativity is summarized. Then equations of state
(EOS) with different stability mechanisms for boson stars will be introduced.
We discuss two stabilizing mechanisms: a standard mass term in the Lagrangian
and a vacuum energy in the potential without a mass term. Finally the mass-
radius curves as well as the compactness are presented.

##  II Theoretical Framework

###  II.1 Scale-invariant equation of state from a classical scalar potential

Assuming a complex scalar field for the description of a bosonic matter, a
suitable Lagrangian reads as follows:

$\mathcal{L}=\partial_{\mu}\phi^{*}\partial^{\mu}\phi+m^{2}\phi^{*}\phi-V$ (1)

where $V$ represents the potential, $m$ the mass of the boson and $\phi$ the
complex scalar field. The equation of motion is then given by

$(\partial_{\mu}\partial^{\mu}+m^{2})\phi=-\frac{\partial
V}{\partial\phi^{*}}.$ (2)

For the investigation of the potential we are going to present analytic
general EOS which only depend on the exponent $n$. In this work we assume an
ideal fluid for the bosons and calculate the EOS adopting a flat space-time.
Assuming flat space-time is justifiable in a local density approach. Without
any interaction potential the radius of curvature of the boson star is of the
same order as the Compton wavelength of the massive boson ($r\propto
M_{P}/m_{b}^{2}$ with the Planck mass $M_{P}$ and the mass of the boson
$m_{b}$) [11], which means space-time is strongly curved. The opposite is
given with a strong interaction potential with an interaction strength
$\lambda$. The radius of curvature increases with $M_{P}/m_{b}\sqrt{\lambda}$
[11], which allows to consider flat space-time. Moreover the scalar field only
varies on a large scale, so the gradient of the field can be neglected. This
is given, when $M_{P}/m_{b}\sqrt{\lambda}\gg 1$ [11]. Starting off with the
energy-momentum tensor to calculate the equation of state

$T_{\mu\nu}=-\eta_{\mu\nu}{\cal
L}+\sum_{\phi,\phi^{*}}\frac{\partial\mathcal{L}}{\partial(\partial^{\mu}\phi_{i})}\partial_{\nu}\phi_{i}.$
(3)

Together with equation (1) we get

$T_{\mu\nu}=-\eta_{\mu\nu}{\cal
L}+\partial_{\mu}\phi^{*}\partial_{\nu}\phi+\partial_{\mu}\phi\partial_{\nu}\phi^{*}+\frac{\partial
V}{\partial(\partial^{\mu}\phi)}\partial_{\nu}\phi+\frac{\partial
V}{\partial(\partial^{\mu}\phi^{*})}\partial_{\nu}\phi^{*}.$ (4)

The last two terms in equation (4) vanish since the potentials used in this
work do not depend on the derivatives of the scalar field. The energy-momentum
tensor reduces to the form for an ideal fluid:

$\displaystyle T^{\mu}_{\nu}=\begin{pmatrix}\varepsilon&0&0&0\\\ 0&p&0&0\\\
0&0&p&0\\\ 0&0&0&p\\\ \end{pmatrix}$ (5)

with the energy density $\varepsilon$ and the pressure $p$. Calculating
$T_{00}$ and $T_{ii}$ and making use of the ansatz $\phi=\phi_{0}e^{-i\omega
t}$ leads to the following equations:

$\displaystyle T_{00}$ $\displaystyle=$
$\displaystyle\phi_{0}^{2}\omega^{2}+m^{2}\phi_{0}^{2}+V$ (6) $\displaystyle
T_{ii}$ $\displaystyle=$
$\displaystyle\phi_{0}^{2}\omega^{2}-m^{2}\phi_{0}^{2}-V$ (7)

where $\omega$ and $t$ denote energy and time respectively. Using the ansatz
together with equation (2) gives a relation between $\omega$ and $m$:

$\omega^{2}\phi_{0}^{2}=\phi_{0}^{2}m^{2}+\frac{\partial
V}{\partial\phi^{*}}\phi^{*}.$ (8)

Making use of this relation leads to the desired expressions for $\varepsilon$
and $p$:

$\displaystyle T_{00}$ $\displaystyle=$
$\displaystyle\varepsilon=2m^{2}\phi_{0}^{2}+\frac{\partial
V}{\partial\phi^{*}}\phi^{*}+V$ (9) $\displaystyle T_{ii}$ $\displaystyle=$
$\displaystyle p=\frac{\partial V}{\partial\phi^{*}}\phi^{*}-V.$ (10)

To obtain dimensionless quantities we divide equations (9) and (10) by a
factor $1/m^{4}$

$\displaystyle\varepsilon^{\prime}$ $\displaystyle=$ $\displaystyle
2\phi_{0}^{{}^{\prime}2}+\frac{\partial
V^{\prime}}{\partial\phi^{*}}\phi^{*}+V^{\prime}$ (11) $\displaystyle p$
$\displaystyle=$ $\displaystyle\frac{\partial
V^{\prime}}{\partial\phi^{*}}\phi^{*}-V^{\prime}$ (12)

with $\varepsilon^{\prime}=\varepsilon/m^{4}$, $V^{\prime}=V/m^{4}$ and
$p^{\prime}=p/m^{4}$.

###  II.2 Tolman-Oppenheimer-Volkoff (TOV) equations

By solving the TOV equations one obtains the mass and the radius of an compact
object, which is necessary in order to calculate the mass-radius curves:

$\displaystyle\frac{dp}{dr}$ $\displaystyle=$
$\displaystyle-G\frac{m_{r}(r)\varepsilon(r)}{r^{2}}\biggl{(}1+\frac{p(r)}{\varepsilon(r)}\biggr{)}\biggl{(}1+\frac{4\pi
r^{3}p(r)}{m_{r}(r)}\biggr{)}\biggl{(}1-\frac{2Gm_{r}(r)}{r}\biggr{)}^{-1}$
(13) $\displaystyle\frac{dm_{r}(r)}{dr}$ $\displaystyle=$ $\displaystyle 4\pi
r^{2}\varepsilon(r)$ (14)

with the pressure $p(r)$, radius $r$, the gravitational constant $G$, the mass
$m_{r}$ inside a sphere of radius $r$ and the energy density $\varepsilon(r)$.
Since all the calculations are dimensionless, the TOV equations need to be
rescaled. Applying the following scaling relations for the mass and the radius

$\displaystyle m_{r}$ $\displaystyle=$
$\displaystyle(G^{3}\cdot\varepsilon_{0})^{-1/2}m_{r}^{\prime}$ (15)
$\displaystyle r$ $\displaystyle=$
$\displaystyle(G\cdot\varepsilon_{0})^{-1/2}r^{\prime}$ (16)

leads to the following form:

$\displaystyle\frac{dp^{\prime}}{dr^{\prime}}$ $\displaystyle=$
$\displaystyle-\frac{m_{r}^{\prime}\varepsilon^{\prime}}{r^{\prime
2}}\biggl{(}1+\frac{p^{\prime}}{\varepsilon^{\prime}}\biggr{)}\biggl{(}1+\frac{4\pi
r^{\prime
3}p^{\prime}}{m_{r}^{\prime}}\biggr{)}\biggl{(}1-\frac{2m_{r}^{\prime}}{r^{\prime}}\biggr{)}^{-1}$
(17) $\displaystyle\frac{dm_{r}^{\prime}}{dr^{\prime}}$ $\displaystyle=$
$\displaystyle 4\pi r^{\prime 2}\varepsilon^{\prime}.$ (18)

with $\varepsilon_{0}$ being a constant with dimension of an energy density.
Also the pressure and the energy density need to be rescaled:

$\displaystyle p^{\prime}$ $\displaystyle=$
$\displaystyle\frac{p}{\varepsilon_{0}}$ (19)
$\displaystyle\varepsilon^{\prime}$ $\displaystyle=$
$\displaystyle\frac{\varepsilon}{\varepsilon_{0}}.$ (20)

Please note that in natural units the energy density and the pressure have the
same dimension. For e.g. the case of non-interacting massive bosons with a
mass $m$ one would choose $\varepsilon_{0}=m^{4}$.

###  II.3 Dimensionless Equation of State

Since all calculations are independent on units, we need to derive an equation
of state that satisfies this requirement. We use in the following two
different kind of generalized scalar potentials: one with a mass term and one
with a vacuum term but without a mass term. It turns out that these two kinds
of potentials lead to stable compact star configurations. We note in passing
that without a mass and a vacuum term the compact star configurations are
unstable. For the two different stability mechanisms the EOS differ
considerably as discussed below. The starting point for the following
calculations is equation (9), respectively equation (10). For the scalar
potential of the form

$V=\frac{\lambda}{2^{n/2}}\left(\phi^{*}\phi\right)^{n/2}$ (21)

we get

$\displaystyle p^{\prime}$ $\displaystyle=$
$\displaystyle\frac{\lambda^{\prime}}{2^{n/2}}(\phi_{0}^{\prime})^{n}\left(\frac{n}{2}-1\right)$
(22) $\displaystyle\varepsilon^{\prime}$ $\displaystyle=$ $\displaystyle
2\phi_{0}^{\prime
2}+\frac{\lambda^{\prime}}{2^{n/2}}(\phi_{0}^{\prime})^{n}\left(\frac{n}{2}+1\right)$
(23)

with the dimensionless coupling strength $\lambda^{\prime}=\lambda/m^{4-n}$
and the dimensionless scalar field $\phi_{0}^{\prime}=\phi_{0}/m$. In order to
obtain an analytic expression for the EOS, we can express the pressure in
terms of the scalar field and insert it into the expression for the energy
density. The EOS can be further simplified by rescaling the pressure and the
energy density with a factor $\lambda^{\prime}(n/2-1)$ to the form

$\varepsilon^{\prime}=p^{\prime 2/n}+\frac{n+2}{n-2}\,p^{\prime}.$ (24)

It is evident that the EOS (24) is restricted to $n>2$.

For the second form of the EOS studied we introduce a vacuum term $V_{0}$. We
start with equation (9) and (10) but setting the mass term to zero. The
potential in this case is given by

$V=\frac{\lambda}{2^{n/2}}\left(\phi^{*}\phi\right)^{n/2}+V_{0}.$ (25)

We obtain for the pressure and energy density

$\displaystyle\varepsilon$ $\displaystyle=$
$\displaystyle\frac{\lambda}{2^{n/2}}\phi_{0}^{n}\left(\frac{n}{2}+1\right)+V_{0}$
(26) $\displaystyle p$ $\displaystyle=$
$\displaystyle\frac{\lambda}{2^{n/2}}\phi_{0}^{n}\left(\frac{n}{2}-1\right)-V_{0}.$
(27)

Combining these two equations gives an EOS which is independent on the
interaction strength $\lambda$

$\varepsilon=\frac{n+2}{n-2}\,p+\frac{2n}{n-2}\,V_{0}.$ (28)

Interestingly, this EOS is the one for self-bound stars of the form

$p=c_{s}^{2}\left(\varepsilon-\varepsilon_{\text{vac}}\right)$ (29)

where the pressure vanishes at a nonvanishing vacuum energy density
$\varepsilon_{\text{vac}}$. The prefactor $c_{s}^{2}$ stands for the speed of
sound. As one can see from eq. (28) different values of $c_{s}^{2}$, i.e.
different stiffnesses of the EOS, emerge from the chosen value of the power
$n$ in the scalar potential. Note that the self-bound EOS for has been derived
from interacting bosonic matter. One arrives at the MIT bag EOS by setting
$n=4$ in eq. (28) so that $c_{s}^{2}=1/3$. These similar EOS are based on a
completely different descriptions of the matter, demonstrating that the EOS is
composition blind so that the results from general relativity do not depend on
the underlying microphysics, as dictated from the strong equivalence
principle. Rescaling the EOS with $\varepsilon^{\prime}=\varepsilon/V_{0}$ and
$p^{\prime}=p/V_{0}$ results in a dimensionless EOS. A further rescaling with
the factor $2n/(n-2)$ gives the final dimensionless EOS of the form

$\displaystyle\varepsilon^{\prime}=\frac{n+2}{n-2}\,p^{\prime}+1.$ (30)

This EOS is quite similar to the one with the mass term (see equation (24)).
The part linear in $p$ stays the same but the second part differs and is now a
constant.

##  III Results: Mass-Radius Curves and Compactness

By solving the TOV equations together with the derived EOS, we obtain the
corresponding mass-radius curves.

Figure 1: Mass-radius curve for $V\propto\phi^{3}$ with a mass term in the
Lagrangian in dimensionless units. The mass decreases with increasing radius
in the stable branch as $M\propto R^{-1}$.

Figures 1 and 2 show the mass-radius curves for the case with a mass term for
different values of the power $n$. The case $n=3$ is plotted in a separate
figure, see Fig. 1, due to a different magnitude of the maximum mass and the
corresponding radius. Additionally this curve differs in its shape compared to
curves with larger values of $n$, since the mass decreases with increasing
radius, in a fashion known for e.g. neutron stars. Nevertheless its
compactness is lower than the ones for larger values of the power $n$. Note
that the solutions to the left of the maximum are unstable as depicted in Fig.
1. The mass-radius curves for $n=4$ starts for vanishing mass at a
nonvanishing radius, as it is typical for a mass-radius relation with a radius
independent on the mass. On the other hand, the mass-radius curves for the
cases with $n=5$ and larger start at the origin, i.e. at vanishing mass and
radius. The shape of these mass-radius curves look like the ones of self-bound
stars where the mass increases with $R^{3}$. However, the underlying EOS does
not exhibit a nonvanishing energy density at vanishing pressure.

Figure 2: Mass-radius curves for $V\propto\phi^{n}$ for $n=4$ up to $n=100$
with a mass term in the scalar potential, shown in dimensionless units. The
maximum mass and radius decrease with $n$. The mass-radius curve for $n=4$,
i.e. for the standard $\phi^{4}$-potential, has a different shape compared to
the curves with $n>4$ since the radius goes to a constant value when the mass
goes to zero. The other mass-radius curves resemble those of self-bound stars
as they start at the origin.

In summary three types of solutions are identifiable, for $n=3$, $n=4$ and
$n>4$. The curves differ in their behaviour in the limit of small masses. The
curve with $n=3$ goes to infinite radius, the curve with $n=4$ goes to a
constant value and the curves with $n>4$ are going to zero for small masses.
In order to understand their behaviour it is useful to consider the limit of
small pressure $p$ for the EOS. Equation (24) then simplifies to:

$\varepsilon\approx p^{2/n}.$ (31)

These different shapes can be understood by having a look at the mass-radius
relation of a sphere in hydrostatic equilibrium for a polytropic EOS of the
form $p\propto\rho^{\Gamma}$, with $\rho$ being the mass density and $\Gamma$
a constant. In the nonrelativistic limit at low density, the energy density is
simply given by the mass density $\varepsilon=\rho$. The mass $M$ and the
radius $R$ for a polytrope are related by (see e.g. [48]):

$M^{2-\Gamma}\cdot R^{3\Gamma-4}\propto\text{const.}$ (32)

Equation (31) describes a polytropic EOS and thus gives the possibility to
make use of equation (32) by setting $\Gamma=n/2$.

$\boldsymbol{n}$ | $\boldsymbol{3}$ | $\boldsymbol{4}$ | $\boldsymbol{5}$ | $\boldsymbol{6}$ | $\boldsymbol{7}$ | $\boldsymbol{8}$ | $\boldsymbol{\infty}$
---|---|---|---|---|---|---|---
$\boldsymbol{\Gamma}$ | $\frac{3}{2}$ | $2$ | $\frac{5}{2}$ | $3$ | $\frac{7}{2}$ | $4$ | $\infty$
$\boldsymbol{M}$-$\boldsymbol{R}$ relation | $M\propto R^{-1}$ | $R\propto\text{const.}$ | $M\propto R^{7}$ | $M\propto R^{5}$ | $M\propto R^{13/3}$ | $M\propto R^{4}$ | $M\propto R^{3}$
Table 1: Mass-radius relations for different values of $n$ derived by using
equation (32).

In table 1 we calculated the mass-radius relations for different values of
$n$. We can reproduce that the mass decreases in the limit of large radii in
the case of $n=3$ as $M\propto{}R^{-1}$. Furthermore table 1 confirms a
constant value of the radius $R$ for small pressures, respectively small
masses, for the case $n=4$. The remaining cases (with $n>4$) are described by
slightly different mass-radius relations. However, the mass vanishes in all
these cases in the limit of $R\rightarrow 0$, making them look like the mass-
radius curves of self-bound stars. In the limit of $n\to\infty$ one recovers
the mass-radius relation of an incompressible fluid, i.e.
$M/R^{3}=\text{const.}$ as for self-bound stars. Please note that these
relations are valid for small pressure $p$ which corresponds to small masses.
For higher pressures $p$ the first term in equation (24) dominates so that the
mass-radius relations shown in the table 1 do not hold anymore. As seen in
Figs. 1 and 2, this is the case for the configurations close to the maximum
mass.

Another way to understand the shape of the curves is by having a look at their
slope. By rearranging equation (32) we obtain the relation

$\displaystyle\frac{\text{d}\log M}{\text{d}\log
R}=\frac{3\Gamma-4}{\Gamma-2}.$ (33)

We can now identify the right-hand side of the equation with the slope of the
curve $m$. For $n=3$ one finds $m=-1$, for $n=4$ the slope goes to infinity
($m\rightarrow\infty$) and the case $n>4$ gives constant positive values for
$m$. Again we can confirm the shape of the curves with our numerical results.
As already mentioned before the curves for $n>4$ look similar to those of
self-bound stars, where gravity is not needed to ensure stability. self-bound
stars are characterized by a nonvanishing value of the energy density at zero
pressure. The shape of the mass-radius relation is determined by a constant
energy density, ensuring that $M\propto R^{3}$ holds. This explains why the
mass-radius curves of those stars are located at the origin. However, we
stress that the stars here at not purely self-bound, since they need gravity
to remain stable.

Figure 3: Mass-radius curves for $V\propto\phi^{n}$ stabilized with a vacuum
constant in dimensionless units. All curves have a similar shape and the same
behavior in the limit of small masses, i.e. the radius vanishes for small
masses. The maximum mass and radius increase with $n$.

The mass-radius curves for the scalar potential without a mass term but with a
vacuum term are depicted in Fig. 3. Contrary to the mass term case, these
mass-radius curves only lead to one type of solution. The mass vanishes for
small radii and increases with increasing radius. The mass-radius curves for
the different values of $n$ are lying on top of each other for small radii.
This feature originates from the dominant behaviour of the vacuum term in the
EOS at low densities, which is independent on the value of $n$. The energy
density stays nearly constant as it is just given by the vacuum energy density
resulting in a mass-radius curve of the form $M/R^{3}=\text{const.}$, i.e. the
familiar one for self-bound stars. The maximum masses and radii increase with
higher values of $n$. For both types of EOS we also investigated their
dependence on the speed of sound squared $c_{s}^{2}$, which is defined as the
derivative of the pressure with respect to the energy density:

$\displaystyle c_{s}^{2}=\frac{\partial p}{\partial\varepsilon}.$ (34)

This gives the following simple relation between the equation of state and
$c_{s}^{2}$ for the vacuum term case:

$\displaystyle c_{s}^{2}$ $\displaystyle=$ $\displaystyle\frac{n-2}{n+2}$ (35)

while for the case with a mass term one arrives at the relation:

$\displaystyle
c_{s}^{2}=\left(\frac{2}{n}\,p^{\frac{2-n}{n}}+\frac{n+2}{n-2}\right)^{-1}.$
(36)

One realizes that the EOS with a vacuum term results in a constant
$c_{s}^{2}$. The EOS with a mass term has an increasing $c_{s}^{2}$ with
increasing pressure (strictly speaking for $n>2$ which is the case in our
studies). In the limit of $p\to\infty$ one recovers the speed of sound squared
of the case with a vacuum term. Hence, for a given power $n$, $c_{s}^{2}$ will
always be larger for the EOS with a vacuum term compared to the one with a
mass term. Thereby, the EOS with a vacuum term will be stiffer, for a fixed
value of $n$, compared to the one with a mass term. We also see, that
$c_{s}^{2}$ increases monotonically with the power $n$ reaching $c_{s}^{2}=1$
in the limit $n\to\infty$.

Figure 4: The speed of sound squared $c_{s}^{2}$ plotted against the exponent
of the potential $n$. $c_{s}^{2}$ increases with increasing $n$ and goes
asymptotically to $c_{s}^{2}=1$. This is the stiffest possible equation of
state, where the speed of sound is equal to the speed of light.

In figure 4 we plot the results for the maximal value of $c_{s}^{2}$ in the
center of the maximum mass configuration, compared for both cases of the EOS.
We find that the values of $c_{s}^{2}$ are in good agreement but not identical
as expected. This feature indicates that the maximum mass configuration for
the EOS with a mass term has a $c_{s}^{2}$ close to their asymptotic limit for
large pressures. For large values of $n$ the speed of sound squared reaches
the causal limit $c_{s}^{2}=1$ where the speed of sound squared equals the
speed of light.

In addition to $c_{s}^{2}$ we also studied the compactness, respectively the
maximum compactness $C_{\text{max}}$, of the the resulting compact star
configurations, defined as:

$\displaystyle C=\frac{M}{R},$ (37)

with $M$ and $R$ representing the dimensionless quantities calculated by
solving the TOV equations. In the following we demonstrate that
$C_{\text{max}}$ can exceed the maximum value for quark stars
($C_{\text{max}}=0.271$), as well as the compactness needed to place the
photon orbit outside the star ($C_{\text{max}}=1/3$) and that the maximum
compactness goes asymptotically to the limit of causality
$C_{\text{max}}=0.354$ for large values of the power $n$. The existence of a
photon sphere for boson stars increases the probability of misinterpretating
them as black holes, since black holes always have a photon sphere around them
as a tell-tale signature [49].

Figure 5: The maximum compactness for the $V\propto\phi^{n}$ potential with a
mass and a vacuum term plotted against $n$. Furthermore different limits for
compact objects are included in the plot. The red line represents the limit
for the photon sphere being outside the radius of the star ($C=1/3$). The
green line displays the maximum compactness for causal EOS ($C=0.354$) (causal
in the sense that the speed of sound cannot exceed the speed of light). The
blue line represents the maximum compactness for self-bound stars with a
conformal EOS where $c_{s}^{2}=1/3$ ($C=0.271$).

The results for both cases of the EOS are depicted in figure 5. We find that
for higher values of $n$ the maximum compactness increases monotonically. The
steepest slope of the curves is in the range of $n=3$ up to $n=10$. From that
point on the maximum compactness goes asymptotically to an upper limit. This
boundary is given by the limit of causality $c_{s}^{2}=1$ with $C=0.354$,
where the speed of sound is equal to the speed of light.

This feature implies that the compactness can be greater than $C=1/3$. The
photon sphere for the Schwarzschild metric lies at $R=3M$ with a corresponding
minimal compactness of $C=M/R=1/3$. According to Fig. 5 one realizes that the
photon orbit can lie outside of the boson star for sufficiently large values
of $n$, producing a light ring. A similar feature has been seen for solitonic
boson stars where the maximum compactness also reaches asymptotically the one
for the causal limit [47]. We conclude that solitonic boson stars as well as
self-interacting boson stars studied here can constitute black hole mimickers
with a maximal compactness in excess of $C=1/3$.

Figure 5 shows that boson stars can reach higher $C_{\text{max}}$ than
expected. These values exceed those of standard neutron stars (without a phase
transition typically $C_{\text{max}}=0.2$ to 0.3) or other exotic compact
objects like quark stars where $C=0.271$, (see e.g. the discussions in [50,
48]). We find that the maximal compactness for $n=4$ in the vacuum term case
is $C=0.271$ which is in agreement with the values for quark stars quoted
above. This was expected since equation (28) for this value of $n$ is
$\varepsilon=3p+\text{const.}$ which is the equation of state for an
ultrarelativistic ideal gas of quarks with nonvanishing vacuum energy, given
by the MIT bag constant, or more general a conformal EOS with $c_{s}^{2}=1/3$.

##  IV Summary and Conclusions

We investigated the properties of boson stars with a generalized scalar
potential by extending the power-law potential $V\propto\lambda\phi^{n}$ to an
arbitrary value of the exponent $n$. We derive analytic expressions for the
equation of state as input to the TOV equations for static spheres of fluids.
We introduced two ways to stabilize the boson star configurations: by
including a mass term for the bosons and by including a vacuum term in the
scalar potential without a mass term. For both cases we calculated the mass-
radius curves, the speed of sound, as well as their maximal compactness. We
found three different categories of the mass-radius curves for the EOS with a
mass term. They differ in their behaviour in the limit of small masses. The
radius $R$ for the case $n=3$ goes to infinity when the mass $M$ goes to zero.
For the classical $\phi^{4}$-potential, i.e. for $n=4$, the radius goes to a
constant value for small masses. The mass-radius curves for larger values of
$n$ resemble the mass-radius curves of self-bound stars as their mass vanishes
when the radius goes to zero. Nevertheless, these compact star configurations
are not self-bound because gravity is needed to stabilize them. However, we
find that the mass-radius curves for the EOS with a vacuum term turn out to be
like those of self-bound stars. The mass vanishes for small radii but in this
case the EOS has a nonvanishing energy density for a vanishing pressure so
that these compact star configurations are bound without gravity, i.e. they
constitute self-bound star configurations. For both types of EOSs studied we
find that the speed of sound squared $c_{s}^{2}$ goes asymptotically to the
limit of causality $c_{s}^{2}=1$. The case for the EOS with a vacuum constant
and the one with a mass term case differ slightly in their absolute values for
the same value of the power $n$ while the former one turns out to be always
stiffer compared to the latter one. Both cases show an increase of the speed
of sound squared with the power $n$, so that higher powers in $n$ lead to a
stiffer EOS. We calculated also the compactness for several values of the
exponent $n$ in the range between $n=3$ and $n=100$. We demonstrate that the
maximal compactness $C_{\text{max}}$ increases continuously with the power
$n$. The highest increase occurs between $n=3$ and $n=10$. For higher values
of $n$ the maximal compactness goes asymptotically to the limit of causality
for both types of EOS, in line with the increase of the speed of sound squared
$c_{s}^{2}$ with $n$ seen before. The highest values of $C_{\text{max}}$ for a
given power $n$ are reached with the EOS with a vacuum term. Even for low
values of $n$ the compactness can be already close to the conformal limit,
where $c_{s}^{2}=1/3$. Specifically, for the EOS with a vacuum term and $n=4$
we reproduce the result from the equation of state of an ultrarelativistic
ideal gas of quarks with a vacuum term, the MIT bag model, with a compactness
of $C=0.271$. Furthermore, for the EOS with a mass term and the one with a
vacuum term the maximal compactness can reach and even surpass a compactness
of $C=1/3$ marking the compactness needed to place the photon ring outside of
the compact object. This makes boson stars described by a scalar potential
with a power law with the power of $n\gtrsim 20$ black hole mimickers. In
summary, we establish that by extending the scalar potential
$V=\lambda\phi^{n}$ to an arbitrary value of the exponent $n$ we change the
properties, in particular the compactness, of the resulting compact
configurations drastically. We point out that we have investigated properties
like mass and radius of the boson star independently of the mass of the scalar
field boson $m$ as well as the coupling strength $\lambda$, which enables a
generic approach. For a given dark matter model of self-interacting bosonic
dark matter our results can be rescaled to physical units by the simple
rescaling laws given in the derivation of the EOS. By choosing a suitable mass
scale and interaction strength one can obtain boson stars with masses and
radii comparable to (supermassive) black holes and neutron stars, for example.
Our results could be used to study now boson star mergers with the generalized
scalar potential. By choosing different values of the power of the scalar
potential one is now able to study within the same setting boson star
configurations with entirely different properties concerning the mass-radius
relation, the speed of sound squared, and the maximal compactness up to black
hole mimickers and to delineate the impact on the pattern of emitted
gravitational waves and possible signals for their detection in present and
future gravitational wave observatories.

###### Acknowledgements.

The authors acknowledge support by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) through the CRC-TR 211 ’Strong-interaction matter
under extreme conditions’– project number 315477589 – TRR 211.

## References

  * Tulin and Yu [2018] S. Tulin and H.-B. Yu, Phys. Rept. 730, 1 (2018), arXiv:1705.02358 [hep-ph] .
  * Liebling and Palenzuela [2023] S. L. Liebling and C. Palenzuela, Living Rev. Rel. 26, 1 (2023), arXiv:1202.5809 [gr-qc] .
  * Wheeler [1955] J. A. Wheeler, Phys. Rev. 97, 511 (1955).
  * Derrick [1964] G. H. Derrick, J. Math. Phys. 5, 1252 (1964).
  * Kaup [1968] D. J. Kaup, Phys. Rev. 172, 1331 (1968).
  * Ruffini and Bonazzola [1969] R. Ruffini and S. Bonazzola, Phys. Rev. 187, 1767 (1969).
  * Lee and Pang [1987] T. D. Lee and Y. Pang, Phys. Rev. D 35, 3678 (1987).
  * Lee [1987] T. D. Lee, Phys. Rev. D 35, 3637 (1987).
  * Friedberg  _et al._ [1987a] R. Friedberg, T. D. Lee, and Y. Pang, Phys. Rev. D 35, 3640 (1987a).
  * Friedberg  _et al._ [1987b] R. Friedberg, T. D. Lee, and Y. Pang, Phys. Rev. D 35, 3658 (1987b).
  * Colpi  _et al._ [1986] M. Colpi, S. L. Shapiro, and I. Wasserman, Phys. Rev. Lett. 57, 2485 (1986).
  * Agnihotri  _et al._ [2009] P. Agnihotri, J. Schaffner-Bielich, and I. N. Mishustin, Phys. Rev. D 79, 084033 (2009), arXiv:0812.2770 [astro-ph] .
  * Schunck and Torres [2000] F. E. Schunck and D. F. Torres, Int. J. Mod. Phys. D 9, 601 (2000), arXiv:gr-qc/9911038 .
  * Eby  _et al._ [2016] J. Eby, C. Kouvaris, N. G. Nielsen, and L. C. R. Wijewardhana, JHEP 02, 028, arXiv:1511.04474 [hep-ph] .
  * Kouvaris and Nielsen [2015] C. Kouvaris and N. G. Nielsen, Phys. Rev. D 92, 063526 (2015), arXiv:1507.00959 [hep-ph] .
  * Maselli  _et al._ [2017a] A. Maselli, P. Pnigouras, N. G. Nielsen, C. Kouvaris, and K. D. Kokkotas, Phys. Rev. D 96, 023005 (2017a), arXiv:1704.07286 [astro-ph.HE] .
  * Cassing  _et al._ [2023] M. Cassing, A. Brisebois, M. Azeem, and J. Schaffner-Bielich, Astrophys. J. 944, 130 (2023), arXiv:2210.13697 [gr-qc] .
  * Khlopov  _et al._ [1985] M. Khlopov, B. A. Malomed, and I. B. Zeldovich, Mon. Not. Roy. Astron. Soc. 215, 575 (1985).
  * Chang  _et al._ [2019] J. H. Chang, D. Egana-Ugrinovic, R. Essig, and C. Kouvaris, JCAP 03, 036, arXiv:1812.07000 [hep-ph] .
  * Abbott  _et al._ [2017] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 119, 161101 (2017), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ [2021] R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run (2021), arXiv:2111.03606 [gr-qc] .
  * Antoniadis  _et al._ [2013] J. Antoniadis _et al._ , Science 340, 6131 (2013), arXiv:1304.6875 [astro-ph.HE] .
  * Fonseca  _et al._ [2021] E. Fonseca _et al._ , Astrophys. J. Lett. 915, L12 (2021), arXiv:2104.00880 [astro-ph.HE] .
  * Riley  _et al._ [2021] T. E. Riley _et al._ , Astrophys. J. Lett. 918, L27 (2021), arXiv:2105.06980 [astro-ph.HE] .
  * Romani  _et al._ [2022] R. W. Romani, D. Kandel, A. V. Filippenko, T. G. Brink, and W. Zheng, Astrophys. J. Lett. 934, L17 (2022), arXiv:2207.05124 [astro-ph.HE] .
  * Abbott  _et al._ [2020] R. Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 896, L44 (2020), arXiv:2006.12611 [astro-ph.HE] .
  * Cardoso  _et al._ [2016] V. Cardoso, S. Hopper, C. F. B. Macedo, C. Palenzuela, and P. Pani, Phys. Rev. D 94, 084031 (2016), arXiv:1608.08637 [gr-qc] .
  * Maselli  _et al._ [2017b] A. Maselli, S. H. Völkel, and K. D. Kokkotas, Phys. Rev. D 96, 064045 (2017b), arXiv:1708.02217 [gr-qc] .
  * Mark  _et al._ [2017] Z. Mark, A. Zimmerman, S. M. Du, and Y. Chen, Phys. Rev. D 96, 084002 (2017), arXiv:1706.06155 [gr-qc] .
  * Sennett  _et al._ [2017] N. Sennett, T. Hinderer, J. Steinhoff, A. Buonanno, and S. Ossokine, Phys. Rev. D 96, 024002 (2017), arXiv:1704.08651 [gr-qc] .
  * Palenzuela  _et al._ [2017] C. Palenzuela, P. Pani, M. Bezares, V. Cardoso, L. Lehner, and S. Liebling, Phys. Rev. D 96, 104058 (2017), arXiv:1710.09432 [gr-qc] .
  * Bezares  _et al._ [2017] M. Bezares, C. Palenzuela, and C. Bona, Phys. Rev. D 95, 124005 (2017), arXiv:1705.01071 [gr-qc] .
  * Toubiana  _et al._ [2021] A. Toubiana, S. Babak, E. Barausse, and L. Lehner, Phys. Rev. D 103, 064042 (2021), arXiv:2011.12122 [gr-qc] .
  * Bezares  _et al._ [2022] M. Bezares, M. Bošković, S. Liebling, C. Palenzuela, P. Pani, and E. Barausse, Phys. Rev. D 105, 064067 (2022), arXiv:2201.06113 [gr-qc] .
  * Dietrich  _et al._ [2019] T. Dietrich, S. Ossokine, and K. Clough, Class. Quant. Grav. 36, 025002 (2019), arXiv:1807.06959 [gr-qc] .
  * Abbott  _et al._ [2019] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 9, 011001 (2019), arXiv:1805.11579 [gr-qc] .
  * Bauswein  _et al._ [2023] A. Bauswein, G. Guo, J.-H. Lien, Y.-H. Lin, and M.-R. Wu, Phys. Rev. D 107, 083002 (2023), arXiv:2012.11908 [astro-ph.HE] .
  * Lee  _et al._ [2021] B. K. K. Lee, M.-c. Chu, and L.-M. Lin, Astrophys. J. 922, 242 (2021), arXiv:2110.05538 [astro-ph.HE] .
  * Karkevandi  _et al._ [2022] D. R. Karkevandi, S. Shakeri, V. Sagun, and O. Ivanytskyi, Phys. Rev. D 105, 023001 (2022), arXiv:2109.03801 [astro-ph.HE] .
  * Giangrandi  _et al._ [2022] E. Giangrandi, V. Sagun, O. Ivanytskyi, C. Providência, and T. Dietrich, The effects of self-interacting bosonic dark matter on neutron star properties (2022), arXiv:2209.10905 [astro-ph.HE] .
  * Rutherford  _et al._ [2023] N. Rutherford, G. Raaijmakers, C. Prescod-Weinstein, and A. Watts, Phys. Rev. D 107, 103051 (2023), arXiv:2208.03282 [astro-ph.HE] .
  * Shakeri and Karkevandi [2022] S. Shakeri and D. R. Karkevandi, Bosonic Dark Matter in Light of the NICER Precise Mass-Radius Measurements (2022), arXiv:2210.17308 [astro-ph.HE] .
  * Diedrichs  _et al._ [2023] R. F. Diedrichs, N. Becker, C. Jockel, J.-E. Christian, L. Sagunski, and J. Schaffner-Bielich, Tidal Deformability of Fermion-Boson Stars: Neutron Stars Admixed with Ultra-Light Dark Matter (2023), arXiv:2303.04089 [gr-qc] .
  * Akiyama  _et al._ [2019] K. Akiyama _et al._ (Event Horizon Telescope), Astrophys. J. Lett. 875, L1 (2019), arXiv:1906.11238 [astro-ph.GA] .
  * Akiyama  _et al._ [2022] K. Akiyama _et al._ (Event Horizon Telescope), Astrophys. J. Lett. 930, L12 (2022).
  * Olivares  _et al._ [2020] H. Olivares, Z. Younsi, C. M. Fromm, M. De Laurentis, O. Porth, Y. Mizuno, H. Falcke, M. Kramer, and L. Rezzolla, Mon. Not. Roy. Astron. Soc. 497, 521 (2020), arXiv:1809.08682 [gr-qc] .
  * Bošković and Barausse [2022] M. Bošković and E. Barausse, JCAP 02 (02), 032, arXiv:2111.03870 [gr-qc] .
  * Schaffner-Bielich [2020] J. Schaffner-Bielich,  _Compact Star Physics_ (Cambridge University Press, 2020).
  * Claudel  _et al._ [2001] C.-M. Claudel, K. S. Virbhadra, and G. F. R. Ellis, J. Math. Phys. 42, 818 (2001), arXiv:gr-qc/0005050 .
  * Lattimer [2012] J. M. Lattimer, Ann. Rev. Nucl. Part. Sci. 62, 485 (2012), arXiv:1305.3510 [nucl-th] .