# Sidon-Ramsey and $B_{h}$-Ramsey numbers

Manuel A. Espinosa-García  Centro de Ciencias Matemáticas, UNAM Campus
Morelia, Morelia, Mexico<EMAIL_ADDRESS>,  Amanda Montejano  Facultad
de Ciencias, UNAM campus Juriquilla, Querétaro, Mexico
<EMAIL_ADDRESS>,  Edgardo Roldán-Pensado  Centro de Ciencias
Matemáticas, UNAM Campus Morelia, Morelia, Mexico<EMAIL_ADDRESS>and  J.
David Suárez  Facultad de Ciencias, UNAM campus Juriquilla, Querétaro, Mexico
<EMAIL_ADDRESS>
###### Abstract.

For a given positive integer $k$, the Sidon-Ramsey number
$\operatorname{SR}(k)$ is defined as the minimum value of $n$ such that, in
every partition of the set $[1,n]$ into $k$ parts, there exists a part that
contains two distinct pairs of numbers with the same sum. In other words,
there is a part that is not a Sidon set. In this paper, we investigate the
asymptotic behavior of this parameter and two generalizations of it. The first
generalization involves replacing pairs of numbers with $h$-tuples, such that
in every partition of $[1,n]$ into $k$ parts, there exists a part that
contains two distinct $h$-tuples with the same sum. Alternatively, there is a
part that is not a $B_{h}$ set. The second generalization considers the
scenario where the interval $[1,n]$ is substituted with a non-necessarily
symmetric $d$-dimensional box of the form $\prod_{i=1}^{d}[1,n_{i}]$. For the
general case of $h\geq 3$ and non-symmetric boxes, before applying our method
to obtain the Ramsey-type result, we needed to establish an upper bound for
the corresponding density parameter.

This research was supported by CONACyT project 282280 and PAPIIT project
IG100822.

##  1\. Introduction

A subset $S$ of an additive group $G$ is called a _Sidon set_ if the sums of
any two elements (possibly equal) of $S$ are distinct. In other words, if
$x,y,z,w\in S$ satisfy

$x+y=z+w,$

then $\\{x,y\\}=\\{z,w\\}$, which means that the equation above has only
trivial solutions in $S$. For a given $X\subset G$, an important problem is to
determine the maximum size of a Sidon set contained in $X$. This problem has
been mainly studied when $G=\mathbb{Z}$ and $X=[1,n]$. We use
$\operatorname{F_{2}}(n)$ to denote the size of the largest Sidon set
contained in $[1,n]$. It is known that

(1) $n^{1/2}(1-o(1))\leq\operatorname{F_{2}}(n)\leq n^{1/2}+0.998\,n^{1/4}.$

The upper bound has been progressively improved [ET41, Lin69, Cil10], the best
being recently established by Balogh, Füredi and Roy [BFR23]. The lower bound
may be inferred from several known constructions; in particular the one
provided by Singer concerning maximal Sidon sets in $X=G=\mathbb{Z}_{n}$,
where $n=q^{2}+q+1$ and $q$ is a prime power [Sin38]. For more information
about problems related with Sidon sets the reader may consult the survey paper
of O’Bryant [O’B04].

As with many density theorems, there is a Ramsey version of the problem of
maximizing the size of a Sidon set. For a given positive integer $k$, a _Sidon
$k$-partition_ of $X\subset G$ is a partition of $X$ into $k$ parts, all of
which are Sidon sets. Let $\operatorname{SR}(k)$ be the _Sidon-Ramsey number_
, defined as the minimum $n$ such that there is no Sidon $k$-partition of
$[1,n]$. This parameter can be found in different contexts under different
names. For instance, the existence of the Sidon-Ramsey numbers is a
consequence of a theorem of Rado [Rad43] with the matrix

$\begin{pmatrix}1&1&-1&-1&0&0\\\ 1&0&-1&0&1&0\\\ 1&0&0&-1&0&1\end{pmatrix},$

where two new variables and two new equations are introduced to assure that no
trivial solutions to the equation $x+y=z+w$ are considered. Liang et al.
[LLXX13] and Xu et al. [XLL18], using computer assistance, found the exact
values of $\operatorname{SR}(k)$ for $k\leq 5$ and gave specific bounds for
$k\leq 19$. They use the fact that $\operatorname{SR}(k)$ can be bounded from
above using the pigeonhole principle, that is,

(2) $\operatorname{SR}(k)\leq(t-1)k+1\text{ for any }k,t\geq 2\text{
satisfying }(\operatorname{F_{2}}(t)-1)/(t-1)>k.$

By computing the upper bound derived from (2) and (1), along with the well-
known Singer construction of Sidon sets, it can be easily deduced that the
Sidon-Ramsey number $\operatorname{SR}(k)$ behave asymptotically as $k^{2}$
(Theorem 1). Due to the interesting connection between Sidon sets and
$C_{4}$-free graphs, the proof of the lower bound of Theorem 1 closely
resembles the proof presented in [CG75] for Theorem 3. However, we have
included this proof in Section 2 to ensure completeness.

######  Theorem 1.

Let $k$ be a positive integer, then

$k^{2}-O(k^{c})\leq\operatorname{SR}(k)\leq k^{2}+Ck^{3/2}+O(k),$

where $c\leq 1.525$ depends on the distribution of the prime numbers and $C$
can be taken close to $1.996$ and depends on the best upper bound for Sidon
numbers.

In this paper, we investigate analogous results to Theorem 1 corresponding to
some generalizations of $F_{2}(n)$. Specifically, we study the Ramsey-type
parameter for $B_{h}$ sets in the interval $X=[1,n]\subset\mathbb{Z}$ (precise
definitions will be provided in Section 2). Additionally, we explore the
problem in higher dimensions by studying the size of the largest $B_{h}$-set
within a non-necessarily symmetric $d$-dimensional box
$X=\prod_{j=1}^{d}[1,n_{j}]\subset\mathbb{Z}^{d}$.

The paper is organized as follows. In Section 2 we deal with everything
related to the one dimensional case. First, we define a $B_{h}$ set in
$\mathbb{Z}$, its density parameter $\operatorname{F_{h}}(n)$, and its
corresponding Ramsey-type parameter $\operatorname{BR}_{h}(k)$. If $h=2$ this
coincides with the definition of a Sidon set, $\operatorname{F_{2}}(n)$, and
$\operatorname{SR}(k)$. Then, we present the proof of Theorem 1 which states
the asymptotic behavior of $\operatorname{SR}(k)$ (corresponding to the case
$h=2$) separated from the proof of Theorem 2, which states the asymptotic
behavior of $\operatorname{BR}_{h}(k)$ for $h\geq 3$.

Section 3 is devoted to the work in higher dimensions. Given positive integers
$n_{1}\leq\dots\leq n_{d}$, we use $\operatorname{F_{h}}(n_{1},\dots,n_{d})$
to denote the size of the largest $B_{h}$-set contained in a $d$-dimensional
box $[1,n_{1}]\times\dots\times[1,n_{d}]$. If $n_{1}=\dots=n_{d}=n$, we said
that the box is symmetric. Previous research by Lindström [Lin72], Cilleruelo
[Cil10], and Rackham and Šarka [RŠ10] has studied this parameter in specific
scenarios: symmetric boxes with $h=2$, asymmetrical boxes with $h=2$, and
symmetric boxes for $h\geq 2$, respectively. To complete this investigation,
we employ techniques developed by Cilleruelo [Cil10] and Rackham and Šarka
[RŠ10] to establish an upper bound for
$\operatorname{F_{h}}(n_{1},\dots,n_{d})$ in the most general case, where the
box is not necessarily symmetric and $h\geq 3$ (refer to Theorem 5). For the
Ramsey-type version, it is convenient to define
$\operatorname{br}_{h}(n_{1},\dots,n_{d})$ to be the largest positive integer
$k$ such that there is no $k$-partition of $\prod_{i\leq d}[1,n_{i}]$ in which
all its parts are $B_{h}$ sets. Note that, in dimension $d=1$,
$\operatorname{br}_{h}(n)$ is the counterpart of $\operatorname{BR}_{h}(k)$ in
the sense that $\operatorname{br}_{h}(\operatorname{BR}_{h}(k))=k$.

Tables 1 and 2 serve to identify the notation that we will use and locate the
bounds of each parameter in the manuscript.

| $h=2$ | $h\geq 3$
---|---|---
$\operatorname{F_{h}}(n)$ | (1) | (3)
$\operatorname{BR}_{h}(k)$ | Theorem 1 | Theorem 2
Table 1. We use $\operatorname{F_{h}}(n)$ to denote the size of the largest $B_{h}$-set contained in $[1,n]$, and so $\operatorname{F_{2}}(n)$ is the classical Sidon number. The corresponding Ramsey-type parameter is denoted by $\operatorname{BR}_{h}(k)$, and so $\operatorname{SR}(k)=\operatorname{BR}_{2}(k)$. |  | $h=2$ | $h\geq 3$
---|---|---|---
$\operatorname{F_{h}}(n_{1},\dots,n_{d})$ | $n_{1}=\dots=n_{d}=n$ | (6) | (10)
| $N=n_{1}n_{2}\dots n_{d}$ | (8) | Theorem 5
$\operatorname{br}_{h}(n_{1},\dots,n_{d})$ | $n_{1}=\dots=n_{d}=n$ | (13) | (15)
| $N=n_{1}n_{2}\dots n_{d}$ | (14) | Theorem 6
Table 2. We use $\operatorname{F_{h}}(n_{1},\dots,n_{d})$ to denote the size
of the largest $B_{h}$-set contained in a $d$-dimensional box
$[1,n_{1}]\times[1,n_{2}]\times\dots\times[1,n_{d}]$. The corresponding
Ramsey-type parameter is denoted by
$\operatorname{br}_{h}(n_{1},\dots,n_{d})$. When $n_{1}=n_{2}=\dots=n_{d}=n$
we said that the box is symmetric. In other case, the results are expressed in
terms of the product $N=n_{1}n_{2}\dots n_{d}$.

##  2\. $B_{h}$ and $B_{h}$-Ramsey numbers in intervals

A subset $S$ of an additive group $G$ is called a _$B_{h}$ -set_ if all sums
of the form $a_{1}+\dots+a_{h}$, where $a_{1},\dots,a_{h}\in A$, are distinct;
note that a $B_{2}$-set is a Sidon set. We use $\operatorname{F_{h}}(n)$ to
denote the size of the largest $B_{h}$-set contained in $[1,n]$. It is known
that

(3) $(1+o(1))n^{1/h}\leq\operatorname{F_{h}}(n)\leq c_{h}n^{1/h}.$

The lower bound was proved by Bose and Chowla [BC63], while the constant
$c_{h}$ in the upper bound has been successively improved [DR84, Jia93, Che94,
Cil01]. Currently, the best bounds are due to Green [Gre01], who proved that
$c_{3}<1.519$, $c_{4}<1.627$ and
$c_{h}\leq\frac{1}{2e}(h+(\frac{3}{2}+o(1))ln(h))$.

As with Sidon sets, we can define the Ramsey version of this problem. For a
given positive integer $k$ and a subset $X$ of an additive group $G$, a
_$B_{h}$ $k$-partition_ of $X$ is a partition of $X$ into $k$ parts, all of
which are $B_{h}$-sets. Let $\operatorname{BR}_{h}(k)$ be the _$B_{h}$ -Ramsey
number_, defined as the minimum $n$ such that there is no $B_{h}$
$k$-partition of $[1,n]$. We shall note that
$\operatorname{BR}_{2}(k)=\operatorname{SR}(k)$.

###### Proof of Theorem 1.

From the upper bound in (1) and using the pigeonhole principle we deduce that
the interval $[1,n]$ cannot be partitioned into less than

$k\geq\frac{n}{\operatorname{F_{2}}(n)}\geq\frac{n}{n^{1/2}+0.998\,n^{1/4}}$

Sidon sets. This can be simplified, using the geometric series, to

$\displaystyle k$ $\displaystyle\geq\frac{n^{1/2}}{1+0.998\,n^{-1/4}}$
$\displaystyle=n^{1/2}\left(1-0.998\,n^{-1/4}+O\left(n^{-1/2}\right)\right)$
$\displaystyle=\left(n^{1/4}-0.499\right)^{2}+O(1).$

Solving this for $n$ gives

$\displaystyle n$ $\displaystyle\leq\left((k+O(1))^{1/2}+0.499\right)^{4}$
$\displaystyle=(k+O(1))^{2}+1.996(k+O(1))^{3/2}+O(k)$
$\displaystyle=k^{2}+1.996\,k^{3/2}+O(k),$

which provides the upper bound.

To prove the lower bound, notice that a Sidon set in $\mathbb{Z}_{n}$ is also
a Sidon set in $[1,n]\subset\mathbb{Z}$. Therefore, any Sidon-Ramsey
$k$-partition of $\mathbb{Z}_{n}$ induces a Sidon-Ramsey $k$-partition of
$[1,n]$. Singer proved that, if $q$ is a prime power, there exists a Sidon set
$S$ in $\mathbb{Z}_{q^{2}+q+1}$ of size $q+1$ [Sin38]. Assume that
$S=\\{s_{1},s_{2},\dots,s_{q+1}\\}$ and consider the Sidon sets
$S_{i}=S-s_{i}$. Since $|S_{i}|=|S|=q+1$ and $S_{i}\cap S_{j}=\\{0\\}$, with
$i\neq j$, we have that

$\big{|}\bigcup_{i=1}^{q+1}S_{i}\big{|}=(q+1)(q+1)-(q+1)+1=q^{2}+q+1,$

then $\bigcup_{i=1}^{q+1}S_{i}=\mathbb{Z}_{q^{2}+q+1}$ and so
$\bigcup_{i=1}^{q+1}S_{i}$ covers $\mathbb{Z}_{q^{2}+q+1}$. Therefore, we may
construct a Sidon-Ramsey $(q+1)$-partition by taking subsets of the elements
of this cover. Let $k\in\mathbb{N}$ and let $p$ be the largest prime less than
$k$. It is known that $p=k-O(k^{0.525})$ (see e.g. [BHP01]). As there exists a
Sidon-Ramsey $(p+1)$-partition of $[1,p^{2}+p+1]$, we conclude that
$\operatorname{SR}(k)\geq\operatorname{SR}(p+1)\geq
p^{2}+p+1=k^{2}-O(k^{1.525})$, which completes the proof. ∎

For $h\geq 2$ we have a result analogous to the Theorem 1. We can derive an
upper bound of $\operatorname{BR}_{h}(k)$ from Green’s bound [Gre01] and the
pigeonhole principle. The lower bound comes from considering translates of the
$B_{h}$-set constructed by Ruzsa (for $h=2$) [Ruz93] and by Gómez-Trujillo
(for $h>2$) [GT11]. Note that the construction used for the $h=2$ case here,
is not the same as the one used in the proof of Theorem 1, although it gives
the same bound.

######  Theorem 2.

Let $k,h$ be positive integers, $h\geq 2$, then

$k^{\frac{h}{h-1}}-O\left(k^{1+\frac{c}{h-1}}\right)\leq\operatorname{BR}_{h}(k)\leq
C_{h}k^{\frac{h}{h-1}},$

where $c\leq 0.525$ depends on the distribution of the prime numbers and
$C_{h}$ depends on the best upper bound for $\operatorname{F_{h}}$ numbers.

###### Proof.

First, we work the upper bound for the $B_{h}-Ramsey$ numbers. The case $h=2$
follows from Theorem 1. If $h>2$, the best known upper bound for $B_{h}$ sets,
established by Green in 2001[Gre01], is

$F_{h}(n)\leq c_{h}n^{1/h},$

where $c_{h}\leq\frac{1}{2e}\left(h+\left(3/2+o(1)\right)\log(h)\right)$. From
this and using the pigeonhole principle we deduce that the interval $[1,n]$
cannot be partitioned into less than

$k\geq\frac{n}{F_{h}(n)}\geq\frac{n}{c_{h}n^{1/h}}$

$B_{h}$ sets. Solving this for $n$ gives

$n\leq C_{h}k^{\frac{h}{h-1}},$

where $C_{h}=c_{h}^{\frac{h}{h{\color[rgb]{1,0,0}-}1}}$. Which provides the
upper bound in Theorem 2 for $h>2$.

To prove the lower bound in Theorem 2, notice that a $B_{h}$ set in
$\mathbb{Z}_{n}$ is also a $B_{h}$ set in $[1,n]\subset\mathbb{Z}$. We
conclude that any $B_{h}$-Ramsey $k$-partition of $\mathbb{Z}_{n}$ induces a
$B_{h}$-Ramsey $k$-partition of $[1,n]$. The main idea to get the lower
bounds. In general, for any positive integers $a\leq b$ and any function
$f:\mathbb{Z}_{a}\to\mathbb{Z}_{b}$, the set

$C=\\{(t,f(t)):t\in\mathbb{Z}_{a}\\}$

and their translations $C+(0,1)$, $C+(0,2)$, $\dots$, $C+(0,b-1)$ give a
partition of $\mathbb{Z}_{a}\times\mathbb{Z}_{b}$. Ruzsa proved that for any
$p$ prime number, there is a function $f:\mathbb{Z}_{p-1}\to\mathbb{Z}_{p}$,
such that the set

$\\{(a,f(a)):a\in\mathbb{Z}_{p-1}\\}\subseteq\mathbb{Z}_{p-1}\times\mathbb{Z}_{p}$

is a $B_{2}$-set (see [Ruz93]). For $h>2$, Gómez and Trujillo proved that for
any $p$ prime number there exist a function
$f:\mathbb{Z}_{p}\to\mathbb{Z}_{p^{h-1}-1}$ such that the set

$\\{(a,f(a)):a\in\mathbb{Z}_{p}\\}\subseteq\mathbb{Z}_{p}\times\mathbb{Z}_{p^{h-1}-1}$

is a $B_{h}$-set (see [GT11]). Then, for any prime number $p$, there is a
$B_{2}$ Ramsey partition of $[1,p^{2}-p]$ in $p$ parts, and for $h>2$ there is
a $B_{h}$-Ramsey partition of $[1,p^{h}-p]$ in $p^{h-1}-1$ parts.

Let $k\in\mathbb{N}$ and, for $h>2$, let $p$ be the largest prime number such
that $p^{h-1}-1\leq k$ (for $h=2$ we consider $p$ to be the largest prime
number less than $k$). It is known that
$p=(k+1)^{1/(h-1)}-O((k+1)^{0.525/(h-1)})$ (for $h=2$, we have
$p=k-O(k^{0.525})$) (see [Dus99]). In the case that $h>2$ we have that

$\displaystyle BR_{h}(k)$
$\displaystyle\geq\operatorname{BR}_{h}(p^{h-1}-1)\geq p^{h}-p$
$\displaystyle=((k+1)^{1/(h-1)}-O((k+1)^{0.525/(h-1)}))^{h}-(k+1)^{1/(h-1)}+O((k+1)^{0.525/(h-1)})$
$\displaystyle=(k+1)^{h/(h-1)}-O((k+1)^{1+0.525/(h-1)}).$

In the case $h=2$ we have that

$BR_{2}(k)\geq\operatorname{BR}_{2}(p)\geq
p^{2}-p\geq(k-O(k^{0.525}))^{2}-(k-O(k^{0.525}))=k^{2}-O(k^{1.525}).\qed$

##  3\. $B_{h}$ and $B_{h}$-Ramsey numbers in $d$-dimensional boxes

It is an interesting problem to study $B_{h}$-sets in higher dimensions. For a
fixed positive integer $d$ and positive integers $n_{1}\leq\dots\leq n_{d}$,
we seek to bound the largest cardinality
$\operatorname{F_{h}}(n_{1},\dots,n_{d})$ of a $B_{h}$-set in the
$d$-dimensional (not necessarily symmetric) box

$X=[1,n_{1}]\times[1,n_{2}]\times\dots\times[1,n_{d}]=\prod_{j=1}^{d}[1,n_{j}]\subset\mathbb{Z}^{d}.$

Regarding the lower bounds, there is a natural way to map one-dimensional
$B_{h}$-sets to $d$-dimensional $B_{h}$-sets. Any integer $0\leq a\leq
n_{1}\cdots n_{d}-1$ has exactly one representation of the form
$a_{1}+a_{2}n_{1}+a_{3}n_{1}n_{2}+\cdots+a_{d-1}n_{1}n_{2}\dots n_{d-1}$,
where $0\leq a_{i}<n_{i}$. Let $\varphi:\mathbb{Z}\to\mathbb{Z}^{d}$ be such
that $\varphi(a)=(a_{1},a_{2},\dots,a_{d})$. Note that any $B_{h}$-set in
$[0,n_{1}\cdots n_{d}-1]$ gets sent into a $B_{h}$-set in
$[0,n_{1}-1]\times\cdots\times[0,n_{d}-1]$, since
$\varphi(x_{1})+\cdots+\varphi(x_{h})=\varphi(y_{1})+\cdots+\varphi(y_{h})$
implies that $x_{1}+\cdots+x_{h}=y_{1}+\cdots+y_{h}$. Using this property of
$\varphi$, we immediately have that

(4) $\operatorname{F_{h}}(n_{1}\cdots
n_{d})\leq\operatorname{F_{h}}(n_{1},\dots,n_{d}).$

So lower bounds in the $d$-dimensional case can be obtained from lower bounds
in the one-dimensional case. This was observed for Sidon sets by Cilleruelo
[Cil10].

As for the upper bounds, the first result was given by Lindström for the case
$n_{1}=\dots=n_{d}=n$ and $h=2$ [Lin72], i.e. Sidon sets in symmetrical boxes
in high dimensions. The bound obtained was

(5) $\operatorname{F_{2}}(n,\dots,n)\leq
n^{d/2}+O\left(n^{\frac{d^{2}}{2d+2}}\right),$

which together with (4) and (3) gives

(6) $n^{d/2}(1-o(1))\leq\operatorname{F_{2}}(n,\dots,n)\leq
n^{d/2}+O\left(n^{\frac{d^{2}}{2d+2}}\right).$

The best result obtained for asymmetrical boxes with $h=2$ is due to
Cilleruelo [Cil10], he proved that

(7) $\operatorname{F_{2}}(n_{1},\dots,n_{d})\leq
N^{1/2}\left(1+O\left(\left(\frac{N_{s-1}}{N^{1/2}}\right)^{\frac{1}{d-s+2}}\right)\right),$

where $N_{0}=1$, $N_{i}=\prod_{j=1}^{i}n_{j}$ for $1\leq i\leq d$, $N=N_{d}$
and $s$ is the least integer such that $N^{1/2}\leq n_{s}^{d-s+2}N_{s-1}$.
This match the lower bound obtained by (4) and (3),

(8) $N^{1/2}(1-o(1))\leq\operatorname{F_{2}}(n_{1},\dots,n_{d})\leq
N^{1/2}\left(1+O\left(\left(\frac{N_{s-1}}{N^{1/2}}\right)^{\frac{1}{d-s+2}}\right)\right).$

For $h\geq 2$, the best upper bound so far for symmetrical boxes is due to
Rackham and Šarka [RŠ10] who showed that

(9)
$\operatorname{F_{h}}(n,\dots,n)\leq\begin{cases}t^{\frac{d}{h}}(t!)^{\frac{1}{t}}n^{\frac{d}{h}}+O\left(n^{\frac{d^{2}}{h(d+1)}}\right)&\text{if
$h=2t$,}\\\
t^{\frac{d-1}{h}}(t!)^{\frac{2}{h}}n^{\frac{d}{h}}+O\left(n^{\frac{d^{2}}{h(d+1)}}\right)&\text{if
$h=2t-1$,}\\\ \end{cases}$

which again match the lower bound obtained from (4) and (3),

(10)
$\left(1+o(1)\right)n^{d/h}\leq\operatorname{F_{h}}(n,\dots,n)\leq\begin{cases}t^{\frac{d}{h}}(t!)^{\frac{1}{t}}n^{\frac{d}{h}}+O\left(n^{\frac{d^{2}}{h(d+1)}}\right)&\text{if
$h=2t$,}\\\
t^{\frac{d-1}{h}}(t!)^{\frac{2}{h}}n^{\frac{d}{h}}+O\left(n^{\frac{d^{2}}{h(d+1)}}\right)&\text{if
$h=2t-1$.}\\\ \end{cases}$

For large enough $h$, Rackham and Šarka gave the improvement

(11) $\operatorname{F_{h}}(n,\dots,n)\leq\begin{cases}(\pi
d)^{\frac{d}{2h}}(1+\epsilon(h))t^{\frac{d}{2h}}(t!)^{\frac{2}{h}}n^{\frac{d}{h}}+O\left(n^{\frac{d^{2}}{h(d+1)}}\right)&\text{if
$h=2t$,}\\\ (\pi
d)^{\frac{d}{2h}}(1+\epsilon(h))t^{\frac{d-2}{2h}}(t!)^{\frac{2}{h}}n^{\frac{d}{h}}+O\left(n^{\frac{d^{2}}{h(d+1)}}\right)&\text{if
$h=2t-1$,}\\\ \end{cases}$

where $\epsilon(h)$ is a function that approach $0$ as $h\to\infty$.

Note that the bound (7) by Cilleruelo in the symmetrical case and the bound
(9) by Rackham and Šarka for $h=2$ coincide with the bound (5) by Lindström.
In Theorem 5 we provide an upper bound for
$\operatorname{F_{h}}(n_{1},\dots,n_{d})$ in the most general case, where the
box is not necessarily symmetric and $h\geq 2$. We use a mix of techniques
used by Cilleruelo and Rackham and Šarka. Before continuing we need a couple
of definitions and lemmas.

Given subsets $A$ and $B$ of an additive group $G$, and $g\in G$, we define
the _sumset_ of $A$ and $B$ as

$A+B=\\{a+b:a\in A,b\in B\\},$

and the _additive energy_ between $A$ and $B$ as

$\sum_{g\in G}d_{A}(g)d_{B}(g),$

where

$d_{X}(g)=\lvert\\{(x,x^{\prime}):x,x^{\prime}\in X,x-x^{\prime}=g\\}\rvert.$

######  Lemma 3 ([Cil10]).

Let $G$ be an additive group and let $A,B\subset G$. Then

$\lvert A\rvert^{2}\leq\frac{\lvert A+B\rvert}{\lvert B\rvert^{2}}\sum_{g\in
G}d_{A}(g)d_{B}(g).$

The following lemma is a generalization of Lemmas 4.2.1 and 4.3.1 from [RŠ10].
The two lemmas from [RŠ10] require that $B$ is a symmetrical box, but this
fact is not actually used in the proof. We use $tA$ to denote the set of sums
of $t$ (not necessary different) elements of $A$, and $t\ast A$ to denote the
set of sums of $t$ distinct elements of $A$.

######  Lemma 4.

Let $h\geq 2$, $A$ be a $B_{h}$ set in $\mathbb{N}^{d}$ and
$B=[0,i_{1}-1]\times[0,i_{2}-1]\times\cdots\times[0,i_{d}-1]$.

  1. (1)

If $h=2t$ then

$\sum_{z\in\mathbb{Z}^{d}}d_{tA}(z)d_{B}(z)\leq\lvert B\rvert^{2}+O(\lvert
B\rvert\lvert A\rvert^{h-1}).$

  2. (2)

If $h=2t-1$ then

$\sum_{z\in\mathbb{Z}^{d}}d_{t\ast A}(z)d_{B}(z)\leq\frac{\lvert
A\rvert}{t}\lvert B|^{2}+O(\lvert B|\lvert A|^{h}).$

Using Lemma 3 and Lemma 4 together with techniques previously used by
Cilleruelo [Cil10] and Rackham and Šarka [RŠ10] we obtain the following
theorem.

######  Theorem 5.

Let $n_{1}\leq n_{2}\leq\dots\leq n_{d}$ be positive integers, set $N_{0}=1$,
$N_{i}=\prod_{j=1}^{i}n_{j}$ for $1\leq i\leq d$, $N=N_{d}$, and let $s$ be
the least integer such that $N^{1/h}\leq n_{s}^{d-s+2}N_{s-1}$. Then, for
$h\geq 2$,

$\operatorname{F_{h}}(n_{1},\dots,n_{d})\leq\begin{cases}(t!)^{\frac{2}{h}}t^{\frac{d}{h}}N^{\frac{1}{h}}\left(1+O\left(\left(\frac{N_{s-1}}{N^{\frac{1}{h}}}\right)^{\frac{1}{d-s+2}}\right)\right)&\text{if
$h=2t$,}\\\
(t!)^{\frac{2}{h}}t^{\frac{d-1}{h}}N^{\frac{1}{h}}\left(1+O\left(\left(\frac{N_{s-1}}{N^{\frac{1}{h}}}\right)^{\frac{1}{d-s+2}}\right)\right)&\text{if
$h=2t-1$.}\\\ \end{cases}$

###### Proof.

First, we will work the case $h=2t$. Let $A$ be a $B_{h}$ set in
$X=\prod_{j=1}^{d}[1,n_{j}]\subset\mathbb{Z}^{d}$, $1\leq s\leq d$ and let
$r_{j}=0$ for $j<s$ and $r_{j}=\lfloor tn_{j}M\rfloor$ for $j\geq s$ and some
$0<M<1$ fixed. Applying Lemma 3 with $tA$ and
$B=[0,r_{1}]\times\cdots\times[0,r_{d}]$, and using Lemma 4 we have that

(12) $\frac{\lvert tA\rvert^{2}\lvert B\rvert^{2}}{\lvert
tA+B\rvert}\leq\sum_{z\in\mathbb{Z}^{d}}d_{tA}(z)d_{B}(z)\leq\lvert
B\rvert^{2}+O(\lvert B\rvert\lvert A\rvert^{2t-1}).$

The number of elements in $tA$ is $\binom{\lvert A\rvert-1+t}{t}$. A basic
lower estimation to this binomial coefficient is $\frac{\lvert
A\rvert^{t}}{t!}\leq\binom{\lvert A\rvert-1+t}{t}$, then it follows that
$\frac{\lvert A\rvert^{t}}{t!}\leq\lvert tA\rvert$. Similarly, we get that
$\frac{\lvert A\rvert^{h}}{h!}\leq\lvert hA\rvert$ and $hA\subset hX$, which
implies $\frac{\lvert A\rvert^{h}}{h!}\leq\lvert hX\rvert=h^{d}N$, i.e.
$\lvert A\rvert=O(N^{\frac{1}{h}})$. Then we can estimate (12) as

$\frac{\lvert A\rvert^{2t}\lvert B\rvert^{2}}{(t!)^{2}\lvert
tA+B\rvert}\leq\lvert B\rvert^{2}+O(\lvert B\rvert\lvert
A\rvert^{2t-1})\leq\lvert B\rvert^{2}+O(\lvert B\rvert\lvert
N\rvert^{\frac{h-1}{h}})$

or equivalently

$\lvert A\rvert^{2t}\leq(t!)^{2}\lvert
tA+B\rvert\left(1+O\left(\frac{N^{\frac{h-1}{h}}}{\lvert
B\rvert}\right)\right).$

Notice that $tA+B\subset[1,tn_{1}+r_{1}]\times\cdots\times[1,tn_{d}+r_{d}]$.
As $r_{j}\leq tn_{j}M<r_{j}+1$ then

$\displaystyle\lvert tA+B\rvert$
$\displaystyle\leq\prod_{j=1}^{d}(tn_{j}+r_{j})=t^{d}N\prod_{j=1}^{d}\left(1+\frac{r_{j}}{tn_{j}}\right)=t^{d}N\prod_{j=s}^{d}\left(1+\frac{r_{j}}{tn_{j}}\right)$
$\displaystyle\leq t^{d}N(1+M)^{d-s+1}=t^{d}N(1+O(M)),$

and

$\lvert
B\rvert=\prod_{j=s}^{d}(1+r_{j})\geq\prod_{j=s}^{d}tn_{j}M=t^{d-s+1}\frac{N}{N_{s-1}}M^{d-s+1}.$

We conclude that

$\lvert
A\rvert^{2t}\leq(t!)^{2}t^{d}N(1+O(M))\left(1+O\left(\frac{N^{\frac{h-1}{h}}N_{s-1}}{NM^{d-s+1}}\right)\right).$

To minimize the order of the last expression, take
$M=(N_{s-1}/N^{1/h})^{1/(d-s+2)}$ and $s$ as the least integer such that
$N^{1/h}\leq n_{s}^{d-s+2}N_{s-1}$. Then

$\lvert
A\rvert^{2t}\leq(t!)^{2}t^{d}N\left(1+O\left(\left(\frac{N_{s-1}}{N^{\frac{1}{h}}}\right)^{\frac{1}{d-s+2}}\right)\right).$

We conclude the proof of this case by taking the $2t^{\text{th}}$ root.

The proof in the case $h=2t-1$ is pretty similar to the proof in the case
$h=2t$. Let $A$ be a $B_{h}$-set in $X$, $1\leq s\leq d$ and $0<M<1$. We take
$M$, $r_{j}$ and $B$ as the even case proof. Applying Lemma 3 with $t\ast A$
and $B$ and using Lemma 4 we have that

$\frac{\lvert t\ast A\rvert^{2}\lvert B\rvert^{2}}{\lvert t\ast
A+B\rvert}\leq\sum_{z\in\mathbb{Z}^{d}}d_{t\ast A}(z)d_{B}(z)\leq\frac{\lvert
A\rvert}{t}\lvert B\rvert^{2}+O(\lvert B\rvert\lvert A\rvert^{2t-1}).$

The number of elements in $t\ast A$ is $\binom{A}{t}$. A basic lower
estimation of the binomial coefficient is $\frac{\lvert
A\rvert^{t}}{t!}(1-\frac{c_{t}}{\lvert A\rvert})\leq\lvert t\ast A\rvert$
where $c_{t}$ depends only on $t$. Using that, and the facts that $\lvert
B\rvert=O(N)$ and $\lvert A\rvert=\Omega(N^{\frac{1}{h}})$ (this last follows
from (3) and (4), and from $\frac{\lvert A\rvert^{h}}{h!}\leq h^{d}N$), we
have that

$\displaystyle\frac{\lvert A\rvert^{2t}\left(1-\frac{c_{t}}{\lvert
A\rvert}\right)^{2}\lvert B\rvert^{2}}{(t!)^{2}\lvert t\ast
A+B\rvert}\leq\frac{\lvert A\rvert}{t}\left(\lvert B\rvert^{2}+O(\lvert
B\rvert\lvert A\rvert^{2t-2})\right),$

or equivalently

$\displaystyle\lvert A\rvert^{2t}\leq$ $\displaystyle\frac{(t!)^{2}\lvert
A\rvert}{t}\lvert t\ast A+B\rvert\left(1+O\left(\frac{\lvert
A\rvert^{2t-2}}{\lvert
B\rvert}\right)\right)\left(\frac{1}{1-\frac{c_{t}}{\lvert A\rvert}}\right)$
$\displaystyle\leq$ $\displaystyle\frac{(t!)^{2}\lvert A\rvert}{t}\lvert t\ast
A+B\rvert\left(1+O\left(\frac{\lvert A\rvert^{2t-2}}{\lvert
B\rvert}\right)\right)\left(1+O\left(\frac{1}{\lvert A\rvert}\right)\right)$
$\displaystyle\leq$ $\displaystyle\frac{(t!)^{2}\lvert A\rvert}{t}\lvert t\ast
A+B\rvert\left(1+O\left(\frac{N^{\frac{h-1}{h}}}{\lvert
B\rvert}\right)\right)$

Similarly as in the even case, we get that $\lvert B\rvert\geq
t^{d-s+1}\frac{N}{N_{s-1}}M^{d-s+1}$ and $\lvert t\ast A+B\rvert\leq
t^{d}N{(1+O(M))}$. Then

$\lvert
A\rvert^{2t-1}\leq\frac{(t!)^{2}}{t}t^{d}N(1+O(M))\left(1+O\left(\frac{N^{\frac{2t-2}{2t-1}}N_{s-1}}{NM^{d-s+1}}\right)\right).$

In order to minimize the last expression we take
$M=(N_{s-1}/N^{1/h})^{1/(d-s+2)}$ and $s$ as the least integer such that
$N^{1/h}\leq n_{s}^{d-s+2}N_{s-1}$. Then

$\lvert
A\rvert^{2t-1}\leq(t!)^{2}t^{d-1}N\left(1+O\left(\left(\frac{N_{s-1}}{N^{\frac{1}{h}}}\right)^{\frac{1}{d-s+2}}\right)\right).$

We conclude the proof taking the $(2t-1)^{\text{th}}$ root. ∎

When $h$ is large enough, it is sometimes possible to improve Theorem 5. If
the box $X=\prod_{j=1}^{d}[1,n_{j}]$ is sufficiently symmetric, we may
consider a slightly larger symmetrical box $X^{\prime}$ that contains $X$,
then it is possible to use (11) instead of (9) to obtain the improvement.

For the Ramsey version it is convenient to bound the size of the partition in
terms of the dimensions of the box. Let
$\operatorname{br}_{h}(n_{1},\dots,n_{d})$ be the largest positive integer $k$
such that there is no $B_{h}$ $k$-partition of $\prod_{i\leq d}[1,n_{i}]$. In
dimension $d=1$, $\operatorname{br}_{h}(n)$ is the counterpart of
$\operatorname{BR}_{h}(k)$ in the sense that
$\operatorname{br}_{h}(\operatorname{BR}_{h}(k))=k$.

The second result of this section is a lower and upper bound for
$\operatorname{br}_{h}(n_{1},\dots,n_{d})$.

######  Theorem 6.

Given positive integers $n_{1}\leq n_{2}\leq\dots\leq n_{d}$, let $N_{0}=1$,
$N_{i}=\prod_{j\leq i}n_{j}$ for $1\leq i\leq d$, $N=N_{d}$, and let $s$ be
the least index such that $N^{1/h}\leq n_{s}^{d-s+2}N_{s-1}$. Define

$\gamma_{h}(n_{1},\dots,n_{d})=\begin{cases}\frac{N^{\frac{h-1}{h}}}{\left(t!\right)^{\frac{2}{h}}t^{\frac{d}{h}}}\left(1-O\left(\left(\frac{N_{s-1}}{N^{\frac{1}{h}}}\right)^{\frac{1}{d-s+2}}\right)\right)&\text{if
}h=2t,\\\
\frac{N^{\frac{h-1}{h}}}{\left(t!\right)^{\frac{2}{h}}t^{\frac{d-1}{h}}}\left(1-O\left(\left(\frac{N_{s-1}}{N^{\frac{1}{h}}}\right)^{\frac{1}{d-s+2}}\right)\right)&\text{if
}h=2t-1,\end{cases}$

then

$\gamma_{h}(n_{1},\dots,n_{d})\leq\operatorname{br}_{h}(n_{1},\dots,n_{d})\leq
N^{\frac{h-1}{h}}\left(1+O\left(N\right)^{\frac{-1+c}{h}}\right),$

where $c\leq 0.525$ depends on the distribution of the prime numbers.

###### Proof.

The lower bound is obtained by using the pigeonhole principle with the bounds
given in Theorem 5. The upper bound follows from the lower bound in Theorem 2
and the fact that any $B_{h}$ $k$-partition in $[1,N]$ can be mapped to a
$B_{h}$ $k$-partition in $[1,n_{1}]\times\cdots\times[1,n_{d}]$. ∎

Therefore, the $B_{h}$-Ramsey number
$\operatorname{br}_{h}(n_{1},\dots,n_{d})$ behaves asymptotically as
$N^{(h-1)/h}$.

Note that the lower bound in Theorem 6 depends on the best upper bound for
$F_{h}$. In particular, for large enough $h$, it is sometimes possible to use
(11) instead of (9) to obtain a better lower bound for
$\operatorname{br}_{h}$. To finish the paper we evaluate some particular cases
of Theorem 6. When $h=2$ and the box is symmetric,

(13)
$n^{d/2}-O\left(n^{\frac{d^{2}}{2d+2}}\right)\leq\operatorname{br}_{2}(n,\dots,n)\leq
n^{d/2}+O\left(n^{dc/2}\right).$

When $h=2$ and the box is not necessarily symmetric,

(14)
$N^{1/2}\left(1-O\left(\frac{N_{s-1}}{N^{1/2}}\right)^{\frac{1}{d-s+2}}\right)\leq\operatorname{br}_{2}(n_{1},\dots,n_{d})\leq
N^{1/2}+O\left(N^{c/2}\right).$

Finally, when $h\geq 3$ and the box is symmetric

$\gamma_{h}(n,\dots,n)=\begin{cases}\frac{n^{d\frac{(h-1)}{h}}-O\left(n^{\frac{d(h-1)}{h}-\frac{d}{h(d+1)}}\right)}{\left(t!\right)^{\frac{2}{h}}t^{\frac{d}{h}}}&\text{if
}h=2t,\\\
\frac{n^{d\frac{(h-1)}{h}}-O\left(n^{\frac{d(h-1)}{h}-\frac{d}{h(d+1)}}\right)}{\left(t!\right)^{\frac{2}{h}}t^{\frac{d-1}{h}}}&\text{if
}h=2t-1,\end{cases}$

and then,

(15) $\gamma_{h}(n,\dots,n)\leq\operatorname{br}_{h}(n,\dots,n)\leq
n^{\frac{d(h-1)}{h}}+O\left(n^{\frac{d(c+h-2)}{h}}\right).$

## References

  * [BC63] R. C. Bose and S. Chowla, _Theorems in the additive theory of numbers._ , Comment. Math. Helv. (1962/63).
  * [BFR23] J. Balogh, Z. Füredi, and S. Roy, _An upper bound on the size of Sidon sets_ , Amer. Math. Monthly (2023).
  * [BHP01] R. C. Baker, G. Harman, and J. Pintz, _The difference between consecutive primes II_ , Proc. London Math. Soc. (3) 83 (2001), no. 3, 532–562.
  * [CG75] F. R. K. Chung and R. L. Graham, _On multicolor Ramsey numbers for complete bipartite graphs_ , J. Combinatorial Theory Ser. B 18 (1975), 164–169. MR 360329
  * [Che94] S. Chen, _On the size of finite Sidon sequences_ , Proc. Amer. Math. Soc. 121 (1994), no. 2, 353–356. MR 1196162
  * [Cil01] J. Cilleruelo, _New upper bounds for finite $B_{h}$ sequences_, Adv. Math. 159 (2001), no. 1, 1–17. MR 1823838
  * [Cil10] by same author, _Sidon sets in $\mathbb{N}^{d}$_, J. Combin. Theory Ser. A 117 (2010), no. 7, 857–871.
  * [DR84] A. G. D’yachkov and V. V. Rykov, _$B_{s}$ -sequences_, Mat. Zametki 36 (1984), no. 4, 593–601. MR 771237
  * [Dus99] P. Dusart, _The $k$th prime is greater than $k(\ln k+\ln\ln k-1)$ for $k\geq 2$_, Math. Comp. 68 (1999), no. 225, 411–415.
  * [ET41] P. Erdős and P. Turán, _On a problem of Sidon in additive number theory, and on some related problems_ , J. London Math. Soc. 16 (1941), 212–215.
  * [Gre01] B. Green, _The number of squares and $B_{h}[g]$ sets_, Acta Arith. 100 (2001), no. 4, 365–390. MR 1862059
  * [GT11] C. A. Gómez Ruiz and C. A. Trujillo Solarte, _A new construction of modular $B_{h}$-sequences_, Mat. Enseñ. Univ. (N. S.) 19 (2011), no. 1, 53–62. MR 2819036
  * [Jia93] X. D. Jia, _On finite Sidon sequences_ , J. Number Theory 44 (1993), no. 1, 84–92. MR 1219489
  * [Lin69] B. Lindström, _An inequality for $B_{2}$-sequences_, J. Combinatorial Theory 6 (1969), 211–212.
  * [Lin72] by same author, _On $B_{2}$-sequences of vectors_, J. Number Theory 4 (1972), 261–265.
  * [LLXX13] M. Liang, X. Li, B. Xiu, and X. Xu, _On Sidon-Ramsey numbers_ , Journal of Computational and Theoretical Nanoscience 10 (2013), no. 4, 884–887.
  * [O’B04] K. O’Bryant, _A complete annotated bibliography of work related to Sidon sequences_ , The Electronic Journal of Combinatorics (2004), DS11:Jul.
  * [Rad43] R. Rado, _Note on combinatorial analysis_ , Proc. London Math. Soc. (2) 48 (1943), 122–160.
  * [RŠ10] L. Rackham and P. Šarka, _$B_{h}$ sequences in higher dimensions_, Electron. J. Combin. 17 (2010), no. 1, Research Paper 35, 15. MR 2595495
  * [Ruz93] I. Z. Ruzsa, _Solving a linear equation in a set of integers I_ , Acta Arith. 65 (1993), no. 3, 259–282.
  * [Sin38] J. Singer, _A theorem in finite projective geometry and some applications to number theory_ , Trans. Amer. Math. Soc. 43 (1938), no. 3, 377–385.
  * [XLL18] X. Xu, M. Liang, and H. Luo, _Ramsey theory: unsolved problems and results_ , Berlin, Boston: De Gruyter, 2018.


Šrámková et al

*

# Oscillations of Fluid Tori around Neutron Stars

Eva Šrámková  Monika Matuszková  Kateřina Klimovičová  Jiři Horák  Odele
Straub  Gabriela Urbancová  Martin Urbanec  Vladimír Karas  Gabriel Török
Debora Lančová  Institute of Physics, Silesian University in Opava, Czech
Republic  Astronomical Institute, Czech Academy of Sciences, Prague, Czech
Republic  ORIGINS Excellence Cluster, Garching, Germany  Max-Planck-Institut
für extraterrestrische Physik, Garching, Germany<EMAIL_ADDRESS>Šrámková
E  M. Matuszková  K. Klimovičová  J. Horák  O. Straub  G. Urbancová  M.
Urbanec  V. Karas  G. Török  D. Lančová

###### Abstract

We examine the influence of quadrupole moment of a slowly rotating neutron
star (NS) on the oscillations of a fluid accretion disk (torus) orbiting a
compact object the spacetime around which is described by the Hartle-Thorne
geometry. Explicit formulae for non-geodesic orbital epicyclic and precession
frequencies, as well as their simplified practical versions that allow for an
expeditious application of the universal relations determining the NS
properties, are obtained and examined. We demonstrate that the difference in
the accretion disk precession frequencies for NSs of the same mass and angular
momentum, but different oblateness, can reach up to tens of percent. Even
higher differences can arise when NSs with the same mass and rotational
frequency, but different equations of state (EoS), are considered. In
particular, the Lense-Thirring precession frequency defined in the innermost
parts of the accretion region can differ by more than one order of magnitude
across NSs with different EoS. Our results have clear implications for models
of the LMXBs variability.

###### keywords:

Gravitation, Accretion, accretion disks, Dense matter, Equation of state

††articletype: Article Type

, , , , , , , , , and , Oscillations and Precession of Fluid Accretion Disks
around Neutron Stars, , .

##  1 Motivation

Oscillations of accretion disks have been studied in various astrophysical
contexts. Among others, the analysis of disk oscillation modes has been
motivated by the phenomenon of high-frequency quasi-periodic oscillations (HF
QPOs) observed in the light curves of Galactic low-mass X-ray binaries (see,
e.g., [Kotrlová . (2020)] and references therein for a review). Despite the
fact that a large number of these observations have been gathered over the
past three decades, there is yet no commonly accepted model for either black
hole (BH) or neutron star (NS) HF QPOs. Based on several strong arguments,
however, it is generally believed these oscillations originate in orbital
motion in the vicinity of the compact object, and numerous models have been
proposed in this context, often assuming some form of accretion disk
oscillations. Over the years, different scenarios have been assumed in
different models incorporating different oscillatory modes and different types
of accretion flow, either geodesic or non-geodesic. In this work, we focus on
the possibility of epicyclic oscillations in a non-geodesic flow, which is
described by the model of an equilibrium perfect fluid torus. We briefly
report on some results that have been obtained within our large and detailed
presently completed study.

##  2 Model of an oscillating torus

To model the torus, we assume a stationary, purely azimuthal flow of perfect
fluid with constant specific angular momentum $l=-u_{\phi}/u_{t}$. The fluid
has a four-velocity $u^{\mu}$ with only two non-zero components,

$u^{\mu}=A(1,0,0,\Omega),$ (1)

where $A$ is the time component $u^{t}$ and $\Omega$ is the orbital velocity
with respect to a distant observer. One may write

$\displaystyle A=u^{t}=(-g_{tt}-2\Omega
g_{t\varphi}-\Omega^{2}g_{\varphi\varphi})^{-1/2},$ (2)
$\displaystyle\Omega=\frac{u^{\varphi}}{u^{t}}=\frac{g^{t\varphi}-lg^{\varphi\varphi}}{g^{tt}-lg^{t\varphi}}.$
(3)

The perfect fluid with rest-mass density $\rho$, pressure $p$ and the energy
density $e$ is characterized by the stress-energy tensor

$T^{\mu\nu}=(p+e)u^{\mu}u^{\nu}+pg^{\mu\nu}.$ (4)

For a polytropic fluid, we may write

$\displaystyle p=K\rho^{\frac{n+1}{n}},\,\,\,\,\,\,\,e=np+\rho,$ (5)

where $K$ and $n$ denote the polytropic constant and the polytropic index,
respectively. In this work, we use $n=3$, which describes a radiation-
pressure-dominated torus.

The Euler equation is obtained from the energy–momentum conservation law,
$\nabla_{\mu}T^{\mu}_{\;\ \;\nu}=0$,

$\nabla_{\mu}(\ln\epsilon)=-\frac{\nabla_{\mu}p}{p+e},$ (6)

with $\epsilon$ being the specific energy

$\epsilon=-u_{t}=\left(-g^{tt}+2lg^{t\varphi}-l^{2}g^{\varphi\varphi}\right)^{-1/2}.$
(7)

By integrating (6), we obtain the Bernoulli equation in the form

$H\epsilon=\mathrm{const}.,$ (8)

where $H=(p+e)/\rho$ denotes the specific enthalpy. From relation (8), one can
easily derive equations describing the structure and the shape of the torus:

$\displaystyle\frac{p}{\rho}$
$\displaystyle=\frac{p_{0}}{\rho_{0}}f(r,\theta),,\,\,\,\,f(r,\theta)=\frac{1}{nc_{\mathrm{s},0}^{2}}\left[\left(1+nc_{\mathrm{s},0}^{2}\right)\frac{\epsilon_{0}}{\epsilon}-1\right],$
(9)

where $c_{\mathrm{s}}$ is the sound speed in the fluid. The surface of the
torus, which coincides with the surface of zero pressure, is given by the
condition $f(r,\theta)=0$.

###  2.1 Epicyclic oscillations

The concept of epicyclic oscillations within test particle motion is a well
explored topic, see, e.g., [Aliev Galtsov (1981)]. In this work, we are
interested in how the properties of epicyclic oscillations, in particular the
oscillation frequencies, differ for an oscillating perfect fluid torus. Such
problem has been investigated in the past, first in a pseudo-Newtonian study
(Blaes ., 2007) and later it was generalized to Kerr geometry (Straub
Šrámková, 2009). Here we assume fluid tori around rotating NSs and extend the
previous calculations for the case of Hartle-Thorne (HT) metric that well
approximates the spacetime around such objects.

##  3 Hartle-Thorne geometry and its parameters range relevant to rotating NSs

The exterior solution of the HT metric is characterized by three parameters:
the gravitational mass $M$, angular momentum $J$ and the quadrupole moment $Q$
of the star. In this work, we use the HT metric assuming dimensionless forms
of the angular momentum and the quadrupole moment, $j=J/M^{2}$ and
$q=Q/M^{3}$.

The explicit form of the metric, which is too long to be presented here, may
be found in Urbancová . (2019) where there is also a thorough discussion of
the relevance of the HT geometry for the calculations of the geodesic orbital
motion. Here we just briefly summarize the appropriate ranges of the
individual parameters that are implied by the present NS equations of state:
The maximum value of the specific angular momentum of a NS is about
$j_{\mathrm{max}}~{}\sim~{}0.7$, the so-called oblateness parameter defined as
$q/j^{2}$ takes values from $q/j^{2}~{}\sim~{}1.5$ for a very massive
(compact) NS up to $q/j^{2}~{}\sim~{}10$ for a low-mass NS, and the
conservative expectations of the NS mass values are about
$1.4-2.5\,M_{\odot}$.

##  4 Using perturbation method to calculate epicyclic frequencies in thicker
tori

As in the previous studies, we use a perturbation method to calculate the
epicyclic frequencies in fluid tori. We start with the so-called Papaloizou-
Pringle equation, which is a partial differential equation derived by
Papaloizou Pringle (1984) who studied global linear stability of Newtonian
fluid tori with constant specific angular momentum with respect to non-
axisymmetric perturbations. The derived equation represents an eigenvalue
problem that in general case has no analytic solution. It however does have
analytic solution in the specific limit of an infinitely slender torus.

This was utilized by Blaes (1985) who calculated the full spectrum of modes of
a Newtonian infinitely slender torus with constant specific angular momentum.
The epicyclic modes are one of these modes. Later, Abramowicz . (2006) derived
a relativistic version of this equation, which can be written as

$\displaystyle\frac{1}{\sqrt{-g}}\partial_{\mu}\frac{\sqrt{-g}g^{\mu\nu}f^{n}\partial_{\nu}W}{nc_{\mathrm{s},0}^{2}f+1}+\left(l_{0}\omega-m\right)^{2}\frac{\Omega
g^{t\phi}-g^{\phi\phi}}{1-\Omega
l_{0}}\frac{f^{n}}{nc_{\mathrm{s},0}^{2}f+1}W$
$\displaystyle=-\frac{2n\aaa^{2}\left(\overline{\omega}-m\overline{\Omega}\right)^{2}}{\beta^{2}r_{0}^{2}}f^{n-1}W,$
(10)

where $\\{\mu,\nu\\}\in\\{r,\theta\\}$, $\aaa\equiv A/A_{0}$,
$\overline{\Omega}\equiv\Omega/\Omega_{0}$,
$\overline{\omega}\equiv\omega/\Omega_{0}$, $g$ is the determinant of the
metric tensor and $W$ equals to

$W=-\frac{\delta p}{A\rho\left(\omega-m\Omega\right)}.$ (11)

The subscript zero refers to quantities defined at the torus centre,
$r=r_{0}$.

We now expand the perturbation equation around its slender torus limit. In
order to do that, we introduce a dimenssionless parameter (Abramowicz .,
2006),

$\beta^{2}=\frac{2nc_{\mathrm{S},0}^{2}}{r_{0}^{2}\Omega_{0}^{2}A_{0}^{2}},$
(12)

which determines the thickness of the torus. Its values range from zero to
one, $\beta->0$ corresponding to the limit of an infinitely slender torus and
$\beta->1$ corresponding to a torus whose outer edge extends to infinity. We
expand the formulae for relevant quantities in $\beta$ by writing

$\displaystyle Q=Q^{(0)}+\beta Q^{(1)}+\beta^{2}Q^{(2)}+\cdots,$ (13)
$\displaystyle
Q\in\left\\{\overline{\omega},W,\aaa,\overline{\Omega},f\right\\}.$ (14)

Substituting that into the relativistic Papaloizou-Pringle equation and
comparing terms of the same order in $\beta$, we obtain analytic expressions
for the corrections to eigenfunctions and eigenfrequencies of the desired
modes. Doing so leads to very complicated and long expressions, to which we
have found shorter approximations, but here we choose to present our results
in terms of figures.

##  5 Results and Conclusions

We present here the results for the $m=-1$ nonaxisymmetric vertical epicyclic
mode, which corresponds to the often discussed and astrophysically relevant
Lense-Thirring precession. This mode in our study possesses the most
interesting features when studied from the point of view of the quadrupole
moment influence. Figures 1 and 2 display the behaviour of its frequency
calculated for different torus thickness ($\beta$) and NS angular momentum
($j$). Figure 1 shows the frequency behaviour in the case of a small
rotationally induced NS oblateness, while Figure 2 displays the same, but for
a large NS oblateness. One may see from Figure 1 that with growing angular
momentum of the star, the frequency behaviour looks qualitatively similar, but
the frequencies themselves tend to increase. As we move towards higher NS
oblateness (Figure 2), however, the effect of increased frequencies begins to
fade away. This means that the presence of the quadrupole moment suppresses
the effect of rotation, and for highly rotating stars we thus get similar
frequency behaviour as we would for the same star with small or no rotation.
Furthermore, for high NS oblateness, qualitative changes of frequency
behaviour begin to arise with growing values of $j$, up to a point, where the
frequencies take on negative values. This corresponds to a change of sign in
the precession.

Our results indicate that, while the shape of a (non-oscillating) torus is not
much sensitive to the NS quadrupole moment, the frequencies of the epicyclic
modes of tori oscillations are affected significantly. In particular, the
difference of frequencies of oscillations of tori around NSs of the same mass
and angular momentum, but different NS oblateness, can reach up to tens of
percent. Consequently, the large differences can also arise when NSs with the
same mass and rotational frequency, but different EoS, are considered. For
instance, the Lense-Thirring precession frequency in the innermost parts of
the accretion region can differ by more than one order of magnitude among
stars with different EoS.

The presented results have a clear relevance for modeling of the high-
frequency quasi-periodic oscillations, in particular for models that directly
deal with epicyclic modes of tori oscillations, see the work of Török . (2022)
and references therein. The properties of epicyclic oscillations along with
related astrophysical consequences will be explored in detail within our
larger study that is presently being completed.

Figure 1: The Lense-Thirring precession frequency ($m=-1$ non-axisymmetric
vertical epicyclic frequency) calculated for different torus thickness
($\beta$) and NS angular momentum ($j$) in the case of a small rotationally
induced neutron star oblateness ($q/j^{2}=3$).

Figure 2: The same as Figure 1 but for a large rotationally induced NS
oblateness ($q/j^{2}=10$).

## Acknowledgments

We acknowledge the Czech Science Foundation (GAČR) grant No. 21-02685X. We
express our thanks to the INTEREXCELLENCE project No. LTI17018. We furthermore
acknowledge two internal grants of the Silesian University in Opava,
SGS/12,13/2019. KK thanks to the INTER-EXCELLENCE project No. LTT17003. MM and
MU thank to the INTER-EXCELLENCE project No. LTC18058. DL thanks the Student
Grant Foundation of the Silesian University in Opava, Grant No. SGF/1/2020,
which has been carried out within the EU OPSRE project entitled “Improving the
quality of the internal grant scheme of the Silesian University in Opava”,
reg. number: CZ.02.2.69/0.0/0.0/19_073/0016951.

## References

  * Abramowicz . (2006) Abr-etal:2006Abramowicz, MA., Blaes, OM., Horák, J., Kluzniak, W. Rebusco, P.  200603,  Classical and Quantum Gravity2351689-1696.  10.1088/0264-9381/23/5/014
  * Aliev Galtsov (1981) Ali+Gal-etal:1981Aliev, AN. Galtsov, DV.  198110,  General Relativity and Gravitation1310899-912.  10.1007/BF00756068
  * Blaes (1985) blaes:1985Blaes, OM.  198510,  MNRAS216553-563.  10.1093/mnras/216.3.553
  * Blaes . (2007) bla-etal:2007Blaes, OM., Šrámková, E., Abramowicz, MA., Kluźniak, W. Torkelsson, U.  200708,  ApJ6651642-653.  10.1086/519782
  * Kotrlová . (2020) kot-etal:2020Kotrlová, A., Šrámková, E., Török, G. et al.  202011,  A&A643A31.  10.1051/0004-6361/201937097
  * Papaloizou Pringle (1984) pap+pri-etal:1984Papaloizou, JCB. Pringle, JE.  198406,  MNRAS208721-750.  10.1093/mnras/208.4.721
  * Straub Šrámková (2009) str+sram-etal:2009Straub, O. Šrámková, E.  200903,  Classical and Quantum Gravity265055011.  10.1088/0264-9381/26/5/055011
  * Török . (2022) tor-etal:2022Török, G., Kotrlová, A., Matuszková, M., Klimovičová, K., Lančová, D., Urbancová, G. Šrámková, E.  202204,  ApJ929128.  10.3847/1538-4357/ac5ab6
  * Urbancová . (2019) urb-etal:2019Urbancová, G., Urbanec, M., Török, G., Stuchlík, Z., Blaschke, M. Miller, JC.  201906,  ApJ877266.  10.3847/1538-4357/ab1b4c


# Tensor spectrum of turbulence-sourced gravitational waves as a constraint on
graviton mass

Yutong He  Axel Brandenburg  and Aditya Sinha

###### Abstract

We consider a generic dispersive massive gravity theory and numerically study
its resulting modified energy and strain spectra of tensor gravitational waves
(GWs) sourced by (i) fully developed turbulence during the electroweak phase
transition (EWPT) and (ii) forced hydromagnetic turbulence during the QCD
phase transition (QCDPT). The GW spectra are then computed in both spatial and
temporal Fourier domains. We find, from the spatial spectra, that the slope
modifications are weakly dependent on the eddy size at QCDPT, and, from the
temporal spectra, that the modifications are pronounced in the $1$–$10\,{\rm
nHz}$ range – the sensitivity range of the North American Nanohertz
Observatory for Gravitational Waves (NANOGrav) – for a graviton mass $m_{\rm
g}$ in the range $2\times 10^{-23}\,{\rm eV}\lesssim m_{\rm g}c^{2}\lesssim
7\times 10^{-22}\,{\rm eV}$.

##  1 Introduction

The history of gravity theories alternative to general relativity (GR) is
almost as long as that of GR itself. GR propagates massless gravitons
described by a linear dispersion relation $\omega=ck$, where $c$ is the speed
of light. One possibility of modifying GR is by having a nonlinear dispersion
relation of the form $\omega=\sqrt{c^{2}k^{2}+\omega_{\rm 0}^{2}}$, where
$\omega_{\rm 0}^{2}=(m_{\rm
g}c^{2}/\hbar)^{2}+\mathcal{A}c^{\alpha}k^{\alpha}$ is a frequently adopted
form of a nonlinear modification [1, 2]. Here $m_{\rm g}c^{2}/\hbar$ is an
effective mass term for a nonzero graviton mass $m_{\rm g}$ [3], and
$\mathcal{A}$ and $\alpha$ are two Lorentz-violating parameters that do not
contain $m_{\rm g}$ [1, 2]. In this paper, we focus on having an effective
massive graviton term only, i.e., we only consider $\omega_{\rm cut}=m_{\rm
g}c^{2}/\hbar$, and numerically studying the resulting modified gravitational
wave (GW) energy spectra. However, we believe the idea of using the GW energy
spectra as a constraint can be applied for a general, Lorentz-violating,
modification.

So far, the graviton mass $m_{\rm g}$ has been constrained to be $m_{\rm
g}\lesssim 7.7\times 10^{-23}\,{\rm eV}/c^{2}$ by LIGO [5, 4], which is
already tighter than the constraint on the photon mass $m_{\gamma}\lesssim
10^{-18}\,{\rm eV}/c^{2}$ [6]. However, motivations to develop massive gravity
theories remain. One such motivation is that massive gravity can explain the
accelerated expansion of the universe more “naturally” than dark energy, i.e.,
a Yukawa-type gravitational potential of the form $\propto r^{-1}e^{-m_{\rm
g}rc/\hbar}$ [8, 3, 7] arises naturally from many massive gravity theories and
thus dilutes the gravitational strength at large distances without the need
for a cosmological constant $\Lambda$ [9, 10]. For recent reviews, see ref.
[11] on massive gravity in the context of $\Lambda$-related topics, and ref.
[12] on massive gravity theories more comprehensively.

Phenomenologically, GWs offer a clean and direct constraint on the graviton
mass, albeit not the tightest ones [13]. Some of the methods using GWs include
orbital decay of pulsars [8], modified dispersion relation and alternative
polarization modes using pulsar timing arrays [3], waveforms of extreme-mass-
ratio inspirals [14] and black hole ringdowns [15], modified dispersion
relation [16] and standard sirens [17] using LISA. Multi-messenger detection
of GW and electromagnetic (EM) waves from binary neutron star mergers provide
another direct constraint on the propagation speed difference between GW and
EM waves. However, so far the bound using this method is no better than
waveform measurements of GWs alone [18]; see ref. [19] on various graviton
mass bounds.

Although modified GW energy spectra have been studied in the context of
inflation [20], it is not well explored with a turbulent source. In
particular, it is then not clear whether the spectral GW energy is enhanced or
decreased in the presence of a finite graviton mass. One particular point of
uncertainty arises from the fact that the equivalence between temporal and
spatial spectra is now broken. In turbulence theory and numerical simulations,
one computes spatial spectra, but the measurements in wind tunnels is almost
always based on temporal spectra. In that case, the approximate equivalence
between both spectra is accomplished by the fact that a chunk of turbulence
passes by the detector at a certain mean speed $\overline{\bm{u}}$, so the
fluctuations as a function of time $t$ can be translated into a spatial
dependence on $\bm{x}$ through $\bm{x}=\bm{x}_{0}-\overline{\bm{u}}t$,
relative to some reference point $\bm{x}_{0}$. Relic GWs, on the other hand,
come from all directions, so the equivalence between spatial and temporal
spectra is not that obvious—especially when there is dispersion.

Meanwhile, significant progress has been made in the numerical solution of
relic GWs from turbulent sources [21, 22, 23, 24, 25]. A useful tool is the
Pencil Code, a massively parallel public domain code developed by the
community of users for a broad range of applications [26]. It comes with a GW
solver, where the modification to dispersive GWs is straightforward.
Therefore, in this work we consider turbulence-sourced GWs in a generic
massive gravity theory by adding a nonzero graviton mass term $m_{\rm g}$ to
the otherwise massless GW equation and explore its effect on the resulting GW
spectra. Specifically, we consider GWs sourced by fully developed turbulence
with an initial Kolmogorov scaling during the electroweak phase transition
(EWPT) [27] and more realistic hydromagnetic turbulence that may have been
present during the QCD phase transition (QCDPT) [25]. We then compute the
resulting energy and strain spectra in spatial and temporal Fourier domains,
which are now different due to the dispersion relation being nonlinear.

We begin by briefly recalling the basic phenomenology of massive GWs (section
2) and the relevant parameters for GWs produced during EWPTs and QCDPTs
(section 3). We then present the governing equations solved in this paper
(section 4), and turn then to the discussion of our results (section 5). We
conclude in section 6. We use the $(-+++)$ metric signature and set $c=1$,
unless specifically noted otherwise. We also normalize the critical energy
density at the time of GW generation $t_{*}$ to be unity, i.e., $\rho_{\rm
crit}(t_{*})=1$.

##  2 Massive gravity and its phenomenology

For a metric $g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu}$, where $\eta_{\mu\nu}=\rm
diag(-+++)$ is the Minkowski background and $|h_{\mu\nu}|\ll|\eta_{\mu\nu}|$
is some small perturbation, the action of a generic massive gravity theory
then can be written as

$S=\int d^{4}x(\mathcal{L}_{\rm EH}+\mathcal{L}_{m_{\rm g}}+\mathcal{L}_{\rm
mat}),$ (2.1)

where $\mathcal{L}_{\rm EH}$ is the usual Einstein-Hilbert Lagrangian in GR,
$\mathcal{L}_{m_{\rm g}}$ is a Lagrangian containing the graviton mass term
$m_{\rm g}$, and $\mathcal{L}_{\rm mat}$ is a Lagrangian for matter-energy.
They take the forms

$\mathcal{L}_{\rm EH}=\sqrt{-g}R,\;\;\mathcal{L}_{m_{\rm g}}=\frac{1}{4}m_{\rm
g}^{2}\sqrt{-g}\Big{(}h_{\mu\nu}h^{\mu\nu}-\frac{1}{2}h^{2}\Big{)},\;\;\mathcal{L}_{\rm
mat}=-8\pi G\sqrt{-g}h^{\mu\nu}T_{\mu\nu},$ (2.2)

where $R$ is the Ricci scalar and $T_{\mu\nu}$ is the stress-energy tensor.

Note that the modified action we are working with here is the Fierz-Pauli (FP)
action [28], which is the simplest massive gravity theory. Generalizing the FP
action leads to the full dRGT gravity [29, 30], for which the massive graviton
term reads $\mathcal{L}_{m_{\rm g}}\propto\mathcal{U}(g,H)$, where
$\mathcal{U}(g,H)$ is the sum of up to quintic order interaction terms in the
perturbation metric $H_{\mu\nu}$. Conversely, the dRGT gravity reduces to FP
gravity at the leading order. Therefore, since our main objective here is to
deliver a qualitatively modified GW spectrum instead of the exacts of massive
gravity theory, we choose the simpler FP action here.

For a conserved source, the action in equation (2.1) can be minimized by
setting its variation with respect to the metric to be zero, which reads

$0=\frac{\delta S}{\delta g^{\mu\nu}}=\int
d^{4}x\sqrt{-g}\Big{(}G_{\mu\nu}+\frac{1}{2}m_{\rm g}^{2}\bar{h}_{\mu\nu}-8\pi
GT_{\mu\nu}\Big{)},$ (2.3)

where $G_{\mu\nu}\equiv R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R$ is the Einstein
tensor and $\bar{h}_{\mu\nu}\equiv h_{\mu\nu}-\frac{1}{2}\eta_{\mu\nu}h$ is
the trace-reversed perturbation. Keeping only the leading order terms in the
perturbation $h_{\mu\nu}$ and applying the harmonic gauge
$\partial^{\mu}\bar{h}_{\mu\nu}=0$, the Einstein tensor becomes

$G_{\mu\nu}=\frac{1}{2}\Big{(}\partial_{\mu}\partial_{\rho}\bar{h}^{\rho}_{\nu}+\partial_{\nu}\partial_{\rho}\bar{h}^{\rho}_{\mu}-\partial_{\rho}\partial^{\rho}\bar{h}_{\mu\nu}-\eta_{\mu\nu}\partial_{\rho}\partial_{\sigma}\bar{h}^{\rho\sigma}\Big{)}=-\frac{1}{2}\Box\bar{h}_{\mu\nu},$
(2.4)

where
$\Box\equiv\partial_{\rho}\partial^{\rho}=-\partial_{t}^{2}+\bm{\nabla}^{2}$
is the d’Alembert operator.

Next, inserting equation (2.4) into equation (2.3) gives the linearized
equation for massive gravity [8],

$(\Box-m_{\rm g}^{2})\bar{h}_{\mu\nu}=-16\pi GT_{\mu\nu},$ (2.5)

which is the modified GW equation. Note that we have not yet taken into
account the scale factor $a$ characterizing the expansion of the universe –
this will be presented in section 3. We drop the overbar of $\bar{h}_{\mu\nu}$
from here on.

In momentum space, equation (2.5) becomes

$\ddot{h}_{\mu\nu}(\bm{k},t)+(\bm{k}^{2}+m_{\rm
g}^{2})h_{\mu\nu}(\bm{k},t)=16\pi GT_{\mu\nu}(\bm{k},t),$ (2.6)

where the double dots denote a second derivative in time, i.e.,
$\ddot{h}_{\mu\nu}(\bm{k},t)\equiv\partial_{t}^{2}h_{\mu\nu}(\bm{k},t)$. A
direct consequence of a generic massive gravity theory is the modified
dispersion relation, which now takes the form

$\omega=\sqrt{k^{2}+\omega_{\rm cut}^{2}},$ (2.7)

where $\omega_{\rm cut}=m_{\rm g}c^{2}/\hbar$ is an effective mass term for
the graviton [3].

Another important implication of massive gravity is that the FP action given
by equations (2.1) and (2.2) propagates five degrees of freedom, namely two
tensor, two vector, and one scalar mode, whereas standard GR only contains two
tensor modes, $+$ and $\times$. For more detailed studies on the polarization
states of massive gravity, see ref. [31]. In this paper, however, we only
study the spectral modifications of the tensor mode for the following reasons.
First, detection of any extra polarization modes would automatically indicate
a modified gravity theory, whereas modifications of the tensor modes that are
also present in GR are more subtle to discern. Second, although the
construction is different from the FP gravity here, there exists a Lorentz-
violating minimal theory of massive gravity (MTMG) that only carries two
tensor modes of GWs [32, 33]. These highlight the significance of studying
tensor GWs separately. And third, the nonlinear dispersion in equation (2.7)
holds for all modes, which means the spectral modifications of tensor modes
should have qualitatively valid features for additional modes, too. However,
perhaps for a future project, it is also important to investigate the extra
polarization modes in the context of low-frequency GWs, since below around
$10^{-7}\,{\rm Hz}$, these modes can have amplitudes similar to those of the
tensor modes in GR [34].

##  3 Tensor mode gravitational waves from the early universe

Equation (2.5) gives the linearized equation for massive gravity. For the
early universe, we adopt the normalization such that the scale factor at the
time of GW generation, $t_{*}$, is set to be unity, i.e., $a_{*}\equiv
a(t_{*})=1$. The Hubble parameter at that time is $H_{\ast}$, defined to be
$H_{*}=\sqrt{8\pi G\rho_{\rm crit}(t_{*})/3}$, where $\rho_{\rm crit}(t_{*})$
is the critical energy density at time $t_{*}$. We denote $\omega_{0}^{H}=2\pi
H_{*}/a_{0}$ as the Hubble frequency scaled to the present, $m_{\rm
g}^{H}=\hbar\omega_{0}^{H}$ as the graviton mass corresponding to the Hubble
frequency and $\lambdabar_{0}^{H}=1/\omega_{0}^{H}$ as the reduced Compton
wavelength of such a graviton scaled to the present.

Furthermore, we assume an adiabatic expansion of the universe during the
radiation-dominated era. This means $g_{S}(T)\,T^{3}a^{3}(T)$ stays constant
by entropy conservation, where $g_{S}(T)$ is the number of adiabatic degrees
of freedom at temperature $T$. Using $T_{*}^{\rm EW}\approx 100\,{\rm GeV}$
and $g_{S}(T_{*}^{\rm EW})\approx 100$ at EWPT, $T_{*}^{\rm QCD}\approx
150\,{\rm MeV}$ and $g_{S}(T_{*}^{\rm QCD})\approx 15$ at QCDPT, combined with
$T_{0}=2.7\,{\rm K}$ and $g_{S}(T_{0})=3.9$ at the present, we obtain
$a_{0}/a_{*}^{\rm EW}\sim 1.3\times 10^{15}$ and $a_{0}/a_{*}^{\rm QCD}\sim
10^{12}$.

To obtain the Hubble parameter, we need the critical energy density in
physical units. During the radiation era, the critical energy density of the
universe is approximately the radiation energy density, i.e., $\rho_{\rm
crit}(t_{*})\simeq\mathcal{E}_{\rm
rad}(t_{*})=\pi^{2}g_{S}(T_{*})k_{B}^{4}T_{*}^{4}/(30\hbar^{3}c^{3})$, where
$k_{\rm B}$ is the Boltzmann constant. We note that, in agreement with some
earlier work [23, 24, 25], we use now the symbol $\mathcal{E}$ for mean energy
densities, but keep the commonly used symbol $\rho$ for the critical density,
which includes the rest mass density.

For EWPT and QCDPT, we get $H_{*}^{\rm EW}\approx 2.1\times 10^{10}\rm s^{-1}$
and $H_{*}^{\rm QCD}\approx 1.8\times 10^{4}\rm s^{-1}$, respectively. Given
$H_{0}\approx 3.2\times 10^{-18}\rm s^{-1}$ today, we obtain $H_{0}/H_{*}^{\rm
EW}\approx 1.5\times 10^{-28}$ and $H_{0}/H_{*}^{\rm QCD}\approx 1.8\times
10^{-22}$. Using these, we can obtain the energy dilution factor
$(a_{*}/a_{0})^{4}(H_{*}/H_{0})^{2}$, which we should multiply the generated
energy by in order to obtain its present day value in the form
$h_{0}^{2}\Omega_{\rm GW}$ [21, 22], which is independent of the uncertainty
in the normalized present day Hubble parameter $h_{0}\approx 0.7$. Table 1
summarizes the aforementioned scaling factors, together with some useful
parameters for EWPT and QCDPT.

We adopt the conformal scaling such that $h_{\mu\nu}=ah_{\mu\nu}^{\rm phys}$
and $T_{\mu\nu}^{\rm TT}=a^{4}T_{\mu\nu,{\rm phys}}^{\rm TT}$, and the
relation $a(t)\propto t$ during the radiation dominated era, so the GW
equation (2.5) becomes

$\Big{(}\partial_{t}^{2}+m_{\rm
g}^{2}-\bm{\nabla}^{2}\Big{)}h_{\mu\nu}=\frac{16\pi G}{a}T_{\mu\nu}^{\rm TT},$
(3.1)

where the superscript TT denotes the transverse-traceless (TT) projection. We
solve this equation in Fourier space and work with the stress projected onto
the linear polarization basis, where $+$ and $\times$ denote the plus and
cross polarizations; see [35, 36] for details and [21] for the implementation
in the Pencil Code.

Using normalized conformal time $\bar{t}=t/t_{*}$, scaled wave vector
$\bar{\bm{k}}=\bm{k}/H_{*}$, and scaled normalized stress $\bar{T}^{\rm
TT}_{+/\times}=T^{\rm TT}_{+/\times}/\mathcal{E}_{\rm rad}^{*}$, the GW
equation in momentum space (equation (2.6)) becomes

$\Big{(}\partial_{\bar{t}}^{2}+m_{\rm
g}^{2}+\bar{\bm{k}}^{2}\Big{)}h_{+/\times}(\bm{k},t)=\frac{6}{\bar{t}}\,\bar{T}^{\rm
TT}_{+/\times}(\bm{k},t),$ (3.2)

but we will drop the overbars from now on,111Note that the overbars here
indicate scaled quantities, which means they are different from the overbars
that appeared, and were subsequently dropped, in equations (2.3)–(2.5), where
$\bar{h}_{\mu\nu}$ denoted the trace-reversed perturbation. except for one
case where we explicitly compare with $\bm{x}$ and $t$ in physical space. Note
that this is modified upon equation (13) in ref. [21].

Finally, in terms of observables, we define the characteristic strain
amplitude as $h_{\rm rms}=\langle h^{2}\rangle^{1/2}$, where $h^{2}\equiv
h_{+}^{2}+h_{\times}^{2}=h_{\mu\nu}^{2}/2$, as well as the scaled GW energy
$\mathcal{E}_{\rm GW}=\langle\dot{h}^{2}\rangle/6$, where
$\dot{h}^{2}=\dot{h}_{+}+\dot{h}_{\times}$ and
$\dot{h}_{+/\times}=\partial_{t}h_{+/\times}$, but see ref. [21] for a small
correction term that will here be neglected.

Event | $k_{\rm B}T$ | $\omega_{0}^{H}[\rm Hz]$ | $m_{\rm g}^{H}[\,{\rm eV}/c^{2}]$ | $\lambdabar_{0}^{H}[\rm m]$ | $(a_{*}/a_{0})^{4}(H_{*}/H_{0})^{2}$
---|---|---|---|---|---
EWPT | $100\,{\rm GeV}$ | $1.0\times 10^{-4}$ | $6.8\times 10^{-20}$ | $1.9\times 10^{13}$ | $1.6\times 10^{-5}$
QCDPT | $150\,{\rm MeV}$ | $1.1\times 10^{-7}$ | $7.4\times 10^{-23}$ | $1.7\times 10^{16}$ | $3.1\times 10^{-5}$
Table 1: Scaling factors and useful parameters.

##  4 Hydromagnetic turbulent sources

The full set of governing equations for the density $\rho$, velocity field
$\bm{u}$, and magnetic field $\bm{B}$ with $\bm{\nabla}\cdot\bm{B}=0$ in
conformal time and comoving variables [37, 27] are

$\displaystyle\frac{\partial\ln\rho}{\partial t}=$
$\displaystyle-\frac{4}{3}(\bm{\nabla}\cdot\bm{u}+\bm{u}\cdot\bm{\nabla}\ln\rho)+\frac{1}{\rho}[\bm{u}\cdot(\bm{J}\cdot\bm{B})+\eta\bm{J}^{2}]$
(4.1)
$\displaystyle\begin{split}\frac{D\bm{u}}{Dt}=&\frac{\bm{u}}{3}(\bm{\nabla}\cdot\bm{u}+\bm{u}\cdot\bm{\nabla}\ln\rho)-\frac{\bm{u}}{\rho}[\bm{u}\cdot(\bm{J}\times\bm{B})+\eta\bm{J}^{2}]\\\
&-\frac{1}{4}\bm{\nabla}\ln\rho+\frac{3}{4\rho}\bm{J}\times\bm{B}+\frac{2}{\rho}\bm{\nabla}\cdot(\rho\nu\bm{\mathsf{S}})\end{split}$
(4.2) $\displaystyle\frac{\partial\bm{B}}{\partial t}=$
$\displaystyle\bm{\nabla}\times(\bm{u}\times\bm{B}-\eta\bm{J}+\bm{\mathcal{F}}),$
(4.3)

where $\bm{J}=\bm{\nabla}\times\bm{B}$ is the magnetic current,
$D/Dt=\partial/\partial t+\bm{u}\cdot\bm{\nabla}$ is the advective derivative,
$\nu$ is the kinematic viscosity, $\eta$ is the magnetic diffusivity, and
$\mathsf{S}_{ij}=\frac{1}{2}(u_{i,j}+u_{j,i})-\frac{1}{3}\delta_{ij}u_{k,k}$
are the components of the traceless strain tensor $\bm{\mathsf{S}}$.

For the EWPT, we start with a turbulence spectrum that has Kolmogorov scaling
in $k$ space, i.e., the initial condition of the evolving turbulence
considered in ref. [27]. We have a magnetic field of the form

$\bm{B}_{i}(\bm{k})=\bm{B}_{*}[P_{ij}(\bm{k})-i\sigma_{\rm
M}\epsilon_{ijl}\hat{k}_{l}]g_{j}(\bm{k})S(k),$ (4.4)

where $P_{ij}=\delta_{ij}-\hat{k}_{i}\hat{k}_{j}$ is the projection operator,
$\sigma_{\rm M}$ indicates helicity and is set to 1 in our runs,
$\epsilon_{ijl}$ is the Levi-Civita symbol, $g_{j}(\bm{k})$ is the Fourier
transform of a random $\delta$-correlated vector field with Gaussian
fluctuations, i.e.,
$g_{i}(\bm{x})g_{j}(\bm{x}^{\prime})=\delta_{ij}\delta^{3}(\bm{x}-\bm{x}^{\prime})$,
and $S(k)$ determines the spectral shape with

$S(k)=\frac{k_{\rm f}^{-3/2}(k/k_{\rm f})^{\alpha/2-1}}{[1+(k/k_{\rm
f})^{2(\alpha+5/3)}]^{1/4}},$ (4.5)

where $k_{\rm f}$ is the wave number of the energy-carrying eddies and
$\alpha=4$ for a causal spectrum, such that $S(k)\sim k$ for small $k$ and
$S(k)\sim k^{-5/3}$ for large $k$.

For QCDPT, we start with zero magnetic field and, as in ref. [25], apply
instead a forcing function $\bm{\mathcal{F}}$ with

$\bm{\mathcal{F}}(\bm{x},t)={\rm
Re}[\mathcal{N}\tilde{\bm{f}}(\bm{k})\exp(i\bm{k}\cdot\bm{x}+i\varphi)],$
(4.6)

where the wave vector $\bm{k}(t)$ and the phase $\varphi(t)$ change randomly
at each time step. This forcing function is therefore white noise in time and
consists of plane waves with average wave number $k_{\rm f}$ such that
$|\bm{k}|$ lies in an interval $k_{\rm f}-\delta k/2\leq|\bm{k}|<k_{\rm
f}+\delta k/2$ of width $\delta k$. $\mathcal{N}={\cal F}_{0}/\delta t^{1/2}$
is a normalization factor, where $\delta t$ is the time step and ${\cal
F}_{0}$ is varied to achieve a certain magnetic field strength after a certain
time, and
$\tilde{\mathcal{F}}(\bm{k})=(\bm{k}\times\bm{e})/[\bm{k}^{2}-(\bm{k}\cdot\bm{e})^{2}]^{1/2}$
is a nonhelical forcing function. Here $\bm{e}$ is an arbitrary unit vector
that is not aligned with $\bm{k}$. Note that $|\mathcal{F}|^{2}=1$. As in ref.
[25], this forcing is only enabled during the time interval $1\leq t\leq 2$.
The kinetic and magnetic energy densities are defined as $\mathcal{E}_{\rm
K}(t)=\langle\rho\bm{u}^{2}\rangle/2$ and $\mathcal{E}_{\rm
M}(t)=\langle\bm{B}^{2}\rangle/2$.

##  5 Energy spectra from numerical simulations

We use the Pencil Code [26] to solve equation (3.2) together with a forced
magnetic field given by equations (4.4) and (4.5) at EWPT, and together with
equations (4.1)–(4.3) for turbulence at QCDPT. The runs discussed in this
paper are summarized in table 2. The numerical data for our spectra are
publicly available [38]. We recall that in the code, our nondimensional wave
numbers and frequencies correspond, at the present time, to
$1/\lambdabar_{0}^{H}$ and $\omega_{0}^{H}$, respectively. The numerical
resolution for the runs is arranged as follows: all EWPT runs (E1 through F5)
have $1152^{3}$ mesh points each, except for F1, which has $1024^{3}$ points;
and all QCDPT runs (Q1 through P5) have $512^{3}$ mesh points. We also set
$\nu=\eta=5\times 10^{-5}$ for all runs. In table 2 we list four groups of
runs: two that are applied to the EWPT (E1–E4 with $\mathcal{E}_{\rm
M}=7.8\times 10^{-3}$ and F1–F5 with $\mathcal{E}_{\rm M}=5.6\times 10^{-3}$)
and two that are applied to the QCDPT (Q1–Q4 with $\mathcal{E}_{\rm
M}=3.9\times 10^{-2}$ and P1–P5 with $\mathcal{E}_{\rm M}=3.8\times 10^{-2}$).
For each group, we vary the value of $\omega_{\rm cut}$, which does not affect
the values of $\mathcal{E}_{\rm M}$, which are therefore not listed in table
2.

Runs | $k_{1}$ | $k_{\rm f}$ | $\omega_{\rm cut}$ | $\mathcal{E}_{\rm GW}^{\rm sat}$ | $h_{\rm rms}^{\rm sat}$ | $h_{0}^{2}\Omega_{\rm GW}$ | $h_{\rm c}$
---|---|---|---|---|---|---|---
E1 | $100$ | $600$ | 0 | $1.15\times 10^{-10}$ | $9.36\times 10^{-8}$ | $1.89\times 10^{-15}$ | $7.46\times 10^{-23}$
E2 | $100$ | $600$ | 10 | $1.13\times 10^{-10}$ | $1.15\times 10^{-7}$ | $1.85\times 10^{-15}$ | $9.19\times 10^{-23}$
E3 | $100$ | $600$ | 50 | $1.12\times 10^{-10}$ | $9.95\times 10^{-8}$ | $1.85\times 10^{-15}$ | $7.93\times 10^{-23}$
E4 | $100$ | $600$ | 200 | $1.03\times 10^{-10}$ | $6.44\times 10^{-8}$ | $1.69\times 10^{-15}$ | $5.13\times 10^{-23}$
F1 | $1$ | $100$ | 0 | $3.77\times 10^{-9}$ | $8.93\times 10^{-6}$ | $6.20\times 10^{-14}$ | $7.12\times 10^{-21}$
F2 | $1$ | $100$ | 0.3 | $3.84\times 10^{-9}$ | $8.34\times 10^{-6}$ | $6.31\times 10^{-14}$ | $6.65\times 10^{-21}$
F3 | $1$ | $100$ | 1 | $3.83\times 10^{-9}$ | $7.06\times 10^{-6}$ | $6.29\times 10^{-14}$ | $5.63\times 10^{-21}$
F4 | $1$ | $100$ | 3 | $3.74\times 10^{-9}$ | $3.91\times 10^{-6}$ | $6.14\times 10^{-14}$ | $3.12\times 10^{-21}$
F5 | $1$ | $100$ | 10 | $3.68\times 10^{-9}$ | $2.01\times 10^{-6}$ | $6.06\times 10^{-14}$ | $1.60\times 10^{-21}$
Q1 | 0.3 | 2 | 0 | $4.86\times 10^{-4}$ | $4.16\times 10^{-2}$ | $1.50\times 10^{-8}$ | $4.17\times 10^{-14}$
Q2 | 0.3 | 2 | 0.3 | $4.70\times 10^{-4}$ | $3.44\times 10^{-2}$ | $1.45\times 10^{-8}$ | $3.44\times 10^{-14}$
Q3 | 0.3 | 2 | 1 | $3.91\times 10^{-4}$ | $2.05\times 10^{-2}$ | $1.21\times 10^{-8}$ | $2.06\times 10^{-14}$
Q4 | 0.3 | 2 | 3 | $1.99\times 10^{-4}$ | $8.48\times 10^{-3}$ | $6.15\times 10^{-9}$ | $8.49\times 10^{-15}$
Q5 | 0.3 | 2 | 10 | $3.41\times 10^{-5}$ | $1.37\times 10^{-3}$ | $1.05\times 10^{-9}$ | $1.37\times 10^{-15}$
P1 | 1 | 6 | 0 | $5.05\times 10^{-5}$ | $4.62\times 10^{-3}$ | $1.56\times 10^{-9}$ | $4.63\times 10^{-15}$
P2 | 1 | 6 | 0.3 | $5.02\times 10^{-5}$ | $4.48\times 10^{-3}$ | $1.55\times 10^{-9}$ | $4.48\times 10^{-15}$
P3 | 1 | 6 | 1 | $4.81\times 10^{-5}$ | $3.54\times 10^{-3}$ | $1.49\times 10^{-9}$ | $3.55\times 10^{-15}$
P4 | 1 | 6 | 3 | $3.97\times 10^{-5}$ | $2.07\times 10^{-3}$ | $1.23\times 10^{-9}$ | $2.08\times 10^{-15}$
P5 | 1 | 6 | 10 | $1.85\times 10^{-5}$ | $7.71\times 10^{-4}$ | $5.72\times 10^{-10}$ | $7.71\times 10^{-16}$
Table 2: Summary of runs shown in the paper. $k_{1}$ is the smallest wave
number in the simulation domain, and $k_{\rm f}$ characterizes the peak
magnetic energy. $h_{0}^{2}\Omega_{\rm GW}$ and $h_{\rm c}$ are the values of
$\mathcal{E}_{\rm GW}^{\rm sat}$ and $h_{\rm rms}^{\rm sat}$ scaled to the
present day, respectively.

###  5.1 Spatial Fourier spectra

Following ref. [21], we define spatial Fourier spectra as integrals over
concentric shells in wave number space (indicated now by a single tilde),
i.e.,

${\rm
Sp}_{k}(h)=\int_{4\pi}\left[\left|\tilde{h}_{+}(\bm{k},t)\right|^{2}+\left|\tilde{h}_{\times}(\bm{k},t)\right|^{2}\right]k^{2}\,{\rm
d}\Omega_{k},$ (5.1)

where ${\rm d}\Omega_{k}$ is the differential over the solid angle in $\bm{k}$
space. An analogous definition applies also to $\dot{h}$, which is used for
calculating $E_{\rm GW}(k,t)={\rm Sp}_{k}(\dot{h})/6$, but see equation (B.36)
of ref. [21] for lower order correction terms that are here neglected.

For EWPT, we first choose the same parameters as in ref. [22], where the
smallest wave number in the simulation domain is $k_{1}=100$ and the peak wave
number is $k_{\rm f}=600$. The resulting GW energy density today,
$h_{0}^{2}\Omega_{\rm GW}(k)$, and the strain today, $h_{\rm c}(k)$, as
functions of wave number $k$ are shown in figure 1(a). We see that for all the
effective mass terms $\omega_{\rm cut}\in\\{10,50,200\\}$ corresponding to
$m_{\rm g}\in\\{6.8\times 10^{-19},3.4\times 10^{-18},1.4\times
10^{-17}\\}\,{\rm eV}/c^{2}$, the spectral modifications are not significant,
even though $\omega_{\rm cut}=200$ already corresponds to an unrealistically
large graviton mass $m_{\rm g}=1.4\times 10^{-17}\,{\rm eV}/c^{2}$.

Next, we explore the same EWPT era but with $k_{1}=1$ and $k_{\rm f}=100$, and
the effective mass term $\omega_{\rm cut}\in\\{0.3,1,3,10\\}$, corresponding
to a range of graviton masses $m_{\rm g}\in\\{2.0\times 10^{-20},6.8\times
10^{-20},2.0\times 10^{-19},6.8\times 10^{-19}\\}\,{\rm eV}/c^{2}$. The
resulting spectra are shown in figure 1(b). Now the spectral differences are
more significant than before, especially towards the lower wave number at
around $k_{1}$. For the spectral energy density $h_{0}^{2}\Omega_{\rm GW}(k)$,
the left panel of figure 1(b) shows spectral shape changes by about a factor
of $k^{1.2}$, from $k^{6/5}$ to $k^{12/5}$. For the strain $h_{\rm c}(k)$, the
right-hand panel of figure 1(b) shows slope changes correspondingly by about a
factor $k^{1.6}$, from $k^{-2/5}$ to $k^{6/5}$ as we increase the graviton
mass. However, the constraint on the graviton mass provided by GWs from EWPT
is not significant overall, as can be seen from table 1, compared to the
existing constraint of $m_{\rm g}\lesssim 7.7\times 10^{-23}\,{\rm eV}/c^{2}$
[5, 4] that we quoted in the introduction.

Since the QCDPT era could provide a constraint on the graviton mass tighter
than EWPT by about three orders of magnitude (table 1), we would like to
explore the spectral behaviors of GWs from the QCDPT era. For this we adopt
the simulation setup of previous work [25], since the only change here is
adding an effective graviton mass term $\omega_{\rm cut}$. We have two series
of runs, the first with $k_{1}=0.3$ and $k_{\rm f}=2$ (runs Q2 to Q5), and the
second with $k_{1}=1$ and $k_{\rm f}=6$ (runs P2 to P5). For both series, we
vary $\omega_{\rm cut}\in\\{0.3,1,3,10\\}$, corresponding to a range of
graviton masses $m_{\rm g}\in\\{2.2\times 10^{-23},7.4\times
10^{-23},2.2\times 10^{-22},7.4\times 10^{-22}\\}\,{\rm eV}/c^{2}$. The
resulting energy and strain spectra can be seen in figures 1(c) and 1(d),
which show remarkably consistent modifications for the two series. We see that
for both cases, as the graviton mass increases, the spectral shapes at lower
wave numbers, around the corresponding value of $k_{1}$, become steeper. In
particular, figure 1(c) shows that for $k_{\rm f}=2$, the energy density
$h_{0}^{2}\Omega_{\rm GW}$ goes from $k^{8/5}$ to $k^{3}$, and the strain
$h_{\rm c}$ goes from $k^{-1/5}$ to $k^{3/2}$; and figure 1(d) shows that for
$k_{\rm f}=6$, the energy density changes from $k^{1}$ to $k^{3}$, and the
strain from $k^{-1/2}$ to $k^{3/2}$, respectively. Therefore, the spectral
modifications are only weakly dependent on the driving wave number.

Note that in figures 1(a)–d, we observe the expected slope correspondence
between GW energy density and strain, i.e., in the small graviton mass limit,
we have $h_{0}^{2}\Omega_{\rm GW}(k)\propto k^{2}h_{c}^{2}(k)$, and in the
large graviton mass limit, the relation changes to $h_{0}^{2}\Omega_{\rm
GW}(k)\propto h_{c}^{2}(k)$.

(a) Runs E1 to E4: $\omega_{\rm cut}=0$ (blue), 10 (green), 50 (orange), and
200 (red).

(b) Runs F2 to F5: $\omega_{\rm cut}=0.3$ (blue), 1 (green), 3 (orange), and
10 (red).

(c) Runs Q2 to Q5: $\omega_{\rm cut}=0.3$ (blue), 1 (green), 3 (orange), and
10 (red).

(d) Runs P2 to P5: $\omega_{\rm cut}=0.3$ (blue), 1 (green), 3 (orange), and
10 (red).

Figure 1: Spectra of GW energy density $h_{0}^{2}\Omega_{\rm GW}(k)$ and
strain $h_{\rm c}(k)$ from EWPT and QCDPT scaled to the present time.

To understand the systematic change from $h_{0}^{2}\Omega_{\rm GW}(k)\propto
k^{3}$ to $h_{\rm c}(k)\propto k^{3/2}$ in the QCD runs, it is important to
recall that for a causal spectrum of the magnetic field with ${\rm
Sp}(\bm{B})\propto k^{4}$ [39], the stress only has a spectrum proportion to
$k^{2}$ [40]. This is because a blue spectrum convolved with itself can only
become white noise. It implies that ${\rm Sp}([k^{2}+\omega_{\rm
cut}^{2}]h)\propto k^{2}$, and therefore

$h_{\rm c}=\sqrt{k\,{\rm Sp}(h)}\propto k^{3/2}/(k^{2}+\omega_{\rm
cut}^{2})\propto k^{3/2}$ (5.2)

for $k\ll\omega_{\rm cut}$. The spectrum of $\dot{h}$, on the other hand, is
given by $\omega^{2}{\rm Sp}(h)\propto k^{2}/(k^{2}+\omega_{\rm cut}^{2})$,
and therefore

$h_{0}^{2}\Omega_{\rm GW}(k)=k\,{\rm Sp}({\dot{h})}\propto
k^{3}/(k^{2}+\omega_{\rm cut}^{2})\propto k^{3}$ (5.3)

for $k\ll\omega_{\rm cut}$. These considerations clearly demonstrate the
obtained trend in figures 1(c) and 1(d) from the $h_{0}^{2}\Omega_{\rm
GW}(k)\propto k^{3}$ scaling to a $h_{\rm c}(k)\propto k^{3/2}$ scaling. We
also see that, relative to the case $\omega_{\rm cut}=0$, both spectra
decrease in amplitude by a factor $1/\omega_{\rm cut}^{2}$. This is well borne
out by the simulations, where we see a drop by $1/100$ in both panels of
figure 1(b) for $\omega_{\rm cut}=10$.

For all the graviton masses considered, the patterns at higher wave numbers
are approximately the same, including a sharp drop by many orders of magnitude
found previously in the massless case discussed in ref. [25]. We return to
this feature further below.

In order to make closer connections to potential detections, we would like to
inspect the temporal Fourier spectra next. In particular, we will show that
both EWPT and QCDPT could induce observable features in the $\sim 10\,{\rm
nHz}$ range, accessible by NANOGrav.

###  5.2 Temporal Fourier spectra

As alluded to in the introduction, spatial and temporal spectra can be
different from each other if there is nonlinear dispersion. To the best of our
knowledge, it is the first time that GW spectra have been obtained from
turbulence simulations in terms of temporal frequency. To compute temporal
spectra, we first Fourier transform $\tilde{h}_{+}(\bm{k},t)$ and
$\tilde{h}_{\times}(\bm{k},t)$ back into real space to obtain time series

$h_{+/\times}(\bm{x},t)=\int\tilde{h}_{+/\times}(\bm{k},t)\,e^{{\rm
i}\bm{k}\cdot\bm{x}}\,{\rm d}^{3}\bm{k}/(2\pi)^{3}.$ (5.4)

We then compute their Fourier transforms to $\omega$ space as
$\tilde{\tilde{h}}_{+/\times}(\bm{x},\omega)=\int
h_{+/\times}(\bm{x},t)\,e^{{\rm i}\omega t}\,{\rm d}t$ at several points
$\bm{x}_{i}$. We finally compute the mean spectrum as

${\rm
Sp}_{\omega}(h)=N^{-1}\sum_{i=1}^{N}\left[\left|\tilde{\tilde{h}}_{+}(\bm{x}_{i},\omega)\right|^{2}+\left|\tilde{\tilde{h}}_{\times}(\bm{x}_{i},\omega)\right|^{2}\right],$
(5.5)

which is an average over $N$ spatial points. In practice, we take $N=1024$,
which is the number of mesh points of the simulations in the $x$ direction.
The results for ${\rm Sp}_{\omega}(h)$ are shown in figure 2, where we also
compare with ${\rm Sp}_{k}(h)$ using both $\omega=k$ (ignoring dispersion,
which is only valid for $\omega_{\rm cut}=0$) and
$\omega=\sqrt{k^{2}+\omega_{\rm cut}^{2}}$ (valid also in the presence of
dispersion, with $\omega_{\rm cut}\neq 0$). Since

$\int{\rm Sp}_{\omega}(h)\,{\rm d}\omega=\int{\rm Sp}_{k}(h)\,{\rm d}k=h_{\rm
rms}^{2},$ (5.6)

we present in the following $(\omega/k)\,{\rm Sp}_{k(\omega)}(h)$, where the
factor ${\rm d}k/{\rm d}\omega=\omega/k$ (for $\omega>\omega_{\rm cut}$) has
been applied to take the effect of dispersion into account.

Figure 2:  Left: temporal strain spectra (solid lines) compared with spatial
strain spectra with (red dashed lines) and without (blue dotted lines) the
transformation to $\omega$ space applied, for $\omega_{\rm cut}=0$, 3, and 10.
Right: space–time diagrams of the strain for runs F1, F4, and F5 with
$\omega_{\rm cut}=0$, 3, and 10. The small white boxes on the lower left
corner of the last two strain diagrams with a size of $2\pi/\omega_{\rm cut}$
are given for orientation. The color bars give the present-day physical strain
multiplied by $10^{18}$.

In figure 2, we also show contour plots of $h(\bm{x},t)$ through an
arbitrarily selected section $\bm{x}=(x,0,0)$ for the same three cases with
$\omega_{\rm cut}=0$, 3, and 10. They would give us a direct impression of how
the GW field would vary in the proximity of the Sun in space and time.

The temporal spectra show a sharp drop near $\omega=\omega_{\rm cut}$, but
${\rm Sp}_{\omega}(h)$ only drops by about two orders of magnitude and does
not vanish completely. Moreover, towards smaller values of $\omega<\omega_{\rm
cut}$, ${\rm Sp}_{\omega}(h)$ begins to rise again. This is probably just a
consequence of having a finite length of the time series, which makes the
statistical error large at small frequencies.

The drop in spectral power below $\omega_{\rm cut}$ also has marked
consequences for the appearance of $h(x,t)$ in real space (right-hand panels
of figure 2, which show smaller scale patches in space and time whose size is
comparable to $2\pi/\omega_{\rm cut}$ in space, but only about half as long as
$2\pi/\omega_{\rm cut}$ in time. For $\omega_{\rm cut}=0$, on the other hand,
we see extended patches in $h(x,t)$ that are all comparable in size to the
scale of the domain. We note at this point that in all three plots of
$h(x,t)$, we have selected a time interval of the length $2\pi$, i.e., the
same interval as in $x$. We select the time interval to be near the end of the
run, but other time intervals look similar.

On the outer axes of figure 2, we have indicated the values scaled to a
physical time corresponding to the QCDPT, using $t=(a_{0}/H_{\ast})\,\bar{t}$
and $x=(ca_{0}/H_{\ast})\,\bar{x}$ with $a_{0}=10^{12}$, $H_{\ast}=1.8\times
10^{4}\,{\rm s}^{-1}$, $c=3\times 10^{8}\,{\rm m}\,{\rm s}^{-1}$ being the
physical dimensions used. We see that the patches can become as short as a
fraction of a year.

In figure 3, we show ${\rm Sp}_{\omega}(h)$ and contour plots of $h(\bm{x},t)$
for runs Q1 and Q3 with $\omega_{\rm cut}=0$ and 1. One of the major features
in the runs shown in figure 2 is the presence of a sharp drop of spectral
power for wave numbers above $2k_{\rm f}\approx 6$. This translates to a
similar drop also in the $\omega$ spectra, but now, when $\omega_{\rm cut}=1$,
GWs exist only over a narrow range of frequencies. This is probably also the
reason why the space–time diagrams of the strain show a somewhat more regular
pattern for $\omega_{\rm cut}=1$.

Finally, we note that in figure 3 between $\omega\sim 10$ and $\omega\sim
100$, the blue dotted lines sharply drop by several orders of magnitude,
exposing a divergence between the black solid lines. This disagreement between
the $\omega$ and $\bm{k}$ spectra is not related to the graviton mass, but it
is primarily related to the sharp drop in ${\rm Sp}(h)$ above
$\bar{\omega}=k_{\rm f}=3$. We have seen such a drop in runs with low $k_{\rm
f}$ for the QCD runs (figures 1(c) and 1(d)). One might think that it could be
related to the Reynolds number not being large enough. However, we have been
able to increase the Reynolds number by a factor of ten and have still not
seen any significant change in this drop. On the other hand, temporal spectra
are not as accurate as spatial ones, which take the entire data cube into
account. They are therefore more easily affected by statistical noise. Given
that the drop extends over several orders of magnitude, the differences
between the spatial and temporal spectra occur where the spatial spectrum is
already very weak and below the noise level of the temporal spectrum. Note
that this feature of a sharp drop is also present in a physically more
realistic turbulence source driven by inflationary magnetogenesis [41]. The
weak temporal signal beyond this drop is therefore most likely no longer
reliable.

Figure 3:  Similar to figure 2, but for runs Q1 and Q3 with $\omega_{\rm
cut}=0$ and 1, respectively. The color bars give the physical strain
multiplied by $10^{14}$.

##  6 Discussions and conclusions

We have seen that, depending on different values of graviton masses
considered, GW energy and strain can exhibit significant modifications, mainly
in terms of spectral shapes. The spectra of relic GWs from EWPT and QCDPT have
the potential to constrain the graviton mass and, in turn, massive gravity
theories. For GWs from fully developed hydromagnetic turbulence at EWPT, the
spectral changes are not so pronounced when $k_{\rm f}=600$ and $k_{1}=100$,
but become more significant when $k_{\rm f}=100$ and $k_{1}=1$, allowing us to
see changes at low wave numbers near $k_{1}$. Specifically, the changes here
are that the energy and strain spectra both become steeper by slightly more
than a $k^{1}$ factor as we increase the graviton mass by about 30 fold. For
QCDPT, we see that GWs induced by forced hydromagnetic turbulence also exhibit
steepening spectral slopes around $k_{1}$ by up to a $k^{2}$ factor as the
graviton mass increases by about 30 fold. Quite remarkably, the spectral
slopes are only weakly affected by different driving wave numbers, meaning
that, even if we are unsure about the exact number of turbulent eddies at the
time of QCDPT, we could still use the GW spectra to constrain the graviton
mass. For both EWPT and QCDPT, the relic GW energy in the low wave number
tails is expected to be lower in the presence of a nonzero graviton mass,
although the values at high wave numbers are roughly unchanged.

The frequency spectra can directly be obtained from wave number spectra
through a simple transformation, which takes $\omega_{\rm cut}\neq 0$ into
account. Our work has demonstrated that the agreement between both types of
spectra is rather good, but frequency spectra obtained from a single detector
can be rather noisy. As the name $\omega_{\rm cut}$ suggests, there is no
signal below the cutoff frequency. Determining this cutoff frequency can be a
sensitive means for constraining the value of the graviton mass. In wave
number space, by contrast, the graviton mass only affects the spectral slope
in the proximity of the associated cutoff wave number. Even in real space, a
finite graviton mass manifests itself through a striking absence of waves with
periods above the cutoff value.

#### Data availability.

The source code used for the simulations of this study, the Pencil Code, is
freely available [26]; see also ref. [38] for the numerical data of the
spectra.

## Acknowledgments

We thank the anonymous referee for useful remarks and suggestions. Support
through the grant 2019-04234 from the Swedish Research Council
(Vetenskapsrådet) is gratefully acknowledged. We acknowledge the allocation of
computing resources provided by the Swedish National Allocations Committee at
the Center for Parallel Computers at the Royal Institute of Technology in
Stockholm.

## References

  * [1] S. Mirshekari, N. Yunes and C. M. Will, _Constraining Generic Lorentz Violation and the Speed of the Graviton with Gravitational Waves_ , _Phys. Rev. D_. 85 (2012), 024041 doi:10.1103/PhysRevD.85.024041 [arXiv:1110.2720 [gr-qc]].
  * [2] J. A. R. Cembranos, M. Coma Díaz and P. Martín-Moruno, _Modified gravity as a diagravitational medium_ , _Phys. Lett. B_. 788 (2019), 336-340 doi:10.1016/j.physletb.2018.10.068 [arXiv:1805.09629 [gr-qc]].
  * [3] K. J. Lee, _Pulsar Timing Arrays and Gravity Tests in the Radiative Regime_ , _Class. Quant. Grav._ 30 (2013), 224016 doi:10.1088/0264-9381/30/22/224016 [arXiv:1404.2090 [astro-ph.CO]].
  * [4] B. P. Abbott et al. [LIGO Scientific and VIRGO], _GW170104: Observation of a 50-Solar-Mass Binary Black Hole Coalescence at Redshift 0.2_ , _Phys. Rev. Lett._ 118 (2017) no.22, 221101 [erratum: Phys. Rev. Lett. 121 (2018) no.12, 129901] doi:10.1103/PhysRevLett.118.221101 [arXiv:1706.01812 [gr-qc]].
  * [5] B. P. Abbott et al. [LIGO Scientific and Virgo], _Tests of general relativity with GW150914_ , _Phys. Rev. Lett._ 116 (2016) no.22, 221101 [erratum: Phys. Rev. Lett. 121 (2018) no.12, 129902] doi:10.1103/PhysRevLett.116.221101 [arXiv:1602.03841 [gr-qc]].
  * [6] P. A. Zyla et al. [Particle Data Group], _Review of Particle Physics_ , PTEP 2020 (2020) no.8, 083C01 doi:10.1093/ptep/ptaa104
  * [7] T. Kahniashvili, A. Kar, G. Lavrelashvili, N. Agarwal, L. Heisenberg and A. Kosowsky, _“Cosmic expansion in extended quasidilaton massive gravity_ , _Phys. Rev. D_. 91 (2015) no.4, 041301 [erratum: Phys. Rev. D 100 (2019) no.8, 089902] doi:10.1103/PhysRevD.91.041301 [arXiv:1412.4300 [astro-ph.CO]].
  * [8] L. S. Finn and P. J. Sutton, _Bounding the mass of the graviton using binary pulsar observations_ , _Phys. Rev. D_. 65 (2002), 044022 doi:10.1103/PhysRevD.65.044022 [arXiv:gr-qc/0109049 [gr-qc]].
  * [9] G. Dvali, G. Gabadadze and M. Shifman, _Diluting cosmological constant via large distance modification of gravity_ , doi:10.1142/9789812776310_0034 [arXiv:hep-th/0208096 [hep-th]].
  * [10] S. V. Babak and L. P. Grishchuk, _Finite range gravity and its role in gravitational waves, black holes and cosmology_ , _Int. J. Mod. Phys. D_. 12 (2003), 1905-1960 doi:10.1142/S0218271803004250 [arXiv:gr-qc/0209006 [gr-qc]].
  * [11] M. Maggiore, _Nonlocal Infrared Modifications of Gravity. A Review_ , _Fundam. Theor. Phys._ 187 (2017), 221-281 doi:10.1007/978-3-319-51700-1_16 [arXiv:1606.08784 [hep-th]].
  * [12] C. de Rham, _Massive Gravity_ , _Living Rev. Rel._ 17 (2014), 7 doi:10.12942/lrr-2014-7 [arXiv:1401.4173 [hep-th]].
  * [13] S. Desai, _Limit on graviton mass from galaxy cluster Abell 1689_ , _Phys. Lett. B_. 778 (2018), 325-331 doi:10.1016/j.physletb.2018.01.052 [arXiv:1708.06502 [astro-ph.CO]].
  * [14] V. Cardoso, G. Castro and A. Maselli, _Gravitational waves in massive gravity theories: waveforms, fluxes and constraints from extreme-mass-ratio mergers_ , _Phys. Rev. Lett._ 121 (2018) no.25, 251103 doi:10.1103/PhysRevLett.121.251103 [arXiv:1809.00673 [gr-qc]].
  * [15] R. Dong and D. Stojkovic, _Gravitational wave echoes from black holes in massive gravity_ , _Phys. Rev. D_. 103 (2021) no.2, 024058 doi:10.1103/PhysRevD.103.024058 [arXiv:2011.04032 [gr-qc]].
  * [16] E. Barausse, E. Berti, T. Hertog, S. A. Hughes, P. Jetzer, P. Pani, T. P. Sotiriou, N. Tamanini, H. Witek and K. Yagi, et al. _Prospects for Fundamental Physics with LISA_ , _Gen. Rel. Grav._ 52 (2020) no.8, 81 doi:10.1007/s10714-020-02691-1 [arXiv:2001.09793 [gr-qc]].
  * [17] E. Belgacem et al. [LISA Cosmology Working Group], _Testing modified gravity at cosmological distances with LISA standard sirens_ , _JCAP_ 07 (2019), 024 doi:10.1088/1475-7516/2019/07/024 [arXiv:1906.01593 [astro-ph.CO]].
  * [18] I. M. Shoemaker and K. Murase, _Constraints from the time lag between gravitational waves and gamma rays: Implications of GW170817 and GRB 170817A_ , _Phys. Rev. D_ 97 (2018) no.8, 083013 doi:10.1103/PhysRevD.97.083013 [arXiv:1710.06427 [astro-ph.HE]].
  * [19] C. de Rham, J. T. Deskins, A. J. Tolley and S. Y. Zhou, _Graviton Mass Bounds_ , _Rev. Mod. Phys._ 89 (2017) no.2, 025004 doi:10.1103/RevModPhys.89.025004 [arXiv:1606.08462 [astro-ph.CO]].
  * [20] T. Fujita, S. Kuroyanagi, S. Mizuno and S. Mukohyama, _Blue-tilted Primordial Gravitational Waves from Massive Gravity_ , _Phys. Lett. B_. 789 (2019), 215-219 doi:10.1016/j.physletb.2018.12.025 [arXiv:1808.02381 [gr-qc]].
  * [21] A. Roper Pol, A. Brandenburg, T. Kahniashvili, A. Kosowsky and S. Mandal, _The timestep constraint in solving the gravitational wave equations sourced by hydromagnetic turbulence_ , _Geophys. Astrophys. Fluid Dynamics_. 114 (2020), 130-161 doi:10.1080/03091929.2019.1653460 [arXiv:1807.05479 [physics.flu-dyn]].
  * [22] A. Roper Pol, S. Mandal, A. Brandenburg, T. Kahniashvili and A. Kosowsky, _Numerical simulations of gravitational waves from early-universe turbulence_ , _Phys. Rev. D_. 102 (2020) no.8, 083512 doi:10.1103/PhysRevD.102.083512 [arXiv:1903.08585 [astro-ph.CO]].
  * [23] T. Kahniashvili, A. Brandenburg, G. Gogoberidze, S. Mandal and A. Roper Pol, _Circular polarization of gravitational waves from early-Universe helical turbulence_ , _Phys. Rev. Res._ 3 (2021) no.1, 013193 doi:10.1103/PhysRevResearch.3.013193 [arXiv:2011.05556 [astro-ph.CO]].
  * [24] A. Brandenburg, Y. He, T. Kahniashvili, M. Rheinhardt and J. Schober, _Relic gravitational waves from the chiral magnetic effect_ , Astrophys. J. 911 (2021) no.1, 110 doi:10.3847/1538-4357/abe4d7 [arXiv:2101.08178 [astro-ph.CO]].
  * [25] A. Brandenburg, E. Clarke, Y. He and T. Kahniashvili, _Can we observe the QCD phase transition-generated gravitational waves through pulsar timing arrays?_ [arXiv:2102.12428 [astro-ph.CO]].
  * [26] Pencil Code Collaboration et al. _The Pencil Code, a modular MPI code for partial differential equations and particles: multipurpose and multiuser-maintained_ , _Journal of Open Source Software_ , 6 (2021), 2807, https://doi.org/10.21105/joss.02807 [arXiv:2009.08231 [astro-ph.IM]].
  * [27] A. Brandenburg, T. Kahniashvili, S. Mandal, A. Roper Pol, A. G. Tevzadze and T. Vachaspati, _Evolution of hydromagnetic turbulence from the electroweak phase transition_ , _Phys. Rev. D_. 96 (2017) no.12, 123528 doi:10.1103/PhysRevD.96.123528 [arXiv:1711.03804 [astro-ph.CO]].
  * [28] M. Fierz and W. Pauli, _On relativistic wave equations for particles of arbitrary spin in an electromagnetic field_ , _Proc. Roy. Soc. Lond. A_. 173 (1939), 211-232 doi:10.1098/rspa.1939.0140.
  * [29] C. de Rham and G. Gabadadze, _Generalization of the Fierz-Pauli Action_ , _Phys. Rev. D_. 82 (2010), 044020 doi:10.1103/PhysRevD.82.044020 [arXiv:1007.0443 [hep-th]].
  * [30] C. de Rham, G. Gabadadze and A. J. Tolley, _Resummation of Massive Gravity_ , _Phys. Rev. Lett._ 106 (2011), 231101 doi:10.1103/PhysRevLett.106.231101 [arXiv:1011.1232 [hep-th]].
  * [31] T. Tachinami, S. Tonosaki and Y. Sendouda, _Gravitational-wave polarizations in generic linear massive gravity and generic higher-curvature gravity_ , _Phys. Rev. D_ 103 (2021) no.10, 104037 doi:10.1103/PhysRevD.103.104037 [arXiv:2102.05540 [gr-qc]].
  * [32] A. De Felice and S. Mukohyama, _Minimal theory of massive gravity_ , _Phys. Lett. B_ 752 (2016), 302-305 doi:10.1016/j.physletb.2015.11.050 [arXiv:1506.01594 [hep-th]].
  * [33] A. De Felice and S. Mukohyama, _Phenomenology in minimal theory of massive gravity_ , _JCAP_ 04 (2016), 028 doi:10.1088/1475-7516/2016/04/028 [arXiv:1512.04008 [hep-th]].
  * [34] W. L. S. de Paula, O. D. Miranda and R. M. Marinho, _Polarization states of gravitational waves with a massive graviton_ , _Class. Quant. Grav._ 21 (2004), 4595-4606 doi:10.1088/0264-9381/21/19/008 [arXiv:gr-qc/0409041 [gr-qc]].
  * [35] M. Maggiore, _Gravitational Waves. Vol. 1: Theory and Experiments_ , Oxford University Press, Oxford U.K. (2007).
  * [36] C. Caprini, R. Durrer and T. Kahniashvili, _The Cosmic microwave background and helical magnetic fields: The Tensor mode_ , _Phys. Rev. D_. 69 (2004), 063006 doi:10.1103/PhysRevD.69.063006 [arXiv:astro-ph/0304556 [astro-ph]].
  * [37] A. Brandenburg, K. Enqvist and P. Olesen, _Large scale magnetic fields from hydromagnetic turbulence in the very early universe_ , _Phys. Rev. D_. 54 (1996), 1291-1300 doi:10.1103/PhysRevD.54.1291 [arXiv:astro-ph/9602031 [astro-ph]].
  * [38] Y. He, A. Brandenburg, and A. Sinha Datasets for “Spectrum of turbulence-sourced gravitational waves as a constraint on graviton mass” (v2021.04.06), doi:10.5281/zenodo.4666074; see also http://www.nordita.org/~brandenb/projects/GravitonGW/ for easier access.
  * [39] R. Durrer and C. Caprini, _Primordial magnetic fields and causality_ , _JCAP_ 11 (2003), 010 doi:10.1088/1475-7516/2003/11/010 [arXiv:astro-ph/0305059 [astro-ph]].
  * [40] A. Brandenburg and S. Boldyrev _The turbulent stress spectrum in the inertial and subinertial ranges_ , Astrophys. J. 892 (2020), 80 [arXiv:1912.07499 [astro-ph.CO]].
  * [41] A. Brandenburg and R. Sharma, _Simulating relic gravitational waves from inflationary magnetogenesis_ , [arXiv:2106.03857 [astro-ph.CO]].