Analysis of Large Scale Web Experiments Using Sequences of Estimators Ian E.
Fellows (Streamlet Data)

###### Abstract

Experimental testing is vital in the optimization of web applications, and as
such A/B testing has been widely adopted as a methodology for determining
optimal content for many web applications. While some testing platforms
provide sequentially valid inferences, a large proportion of online tests
still utilize traditional statistical tests that do not allow for interim
“peeking” at the data or extending the test past its proposed sample size.

In this paper we develop results useful for the sequential analysis of large
scale experiments. In particular, the properties of sequences of maximum
likelihood and generalized method of moments estimators are examined. This
leads to new tests of odds ratios and relative risks for binary outcomes. For
continuous and ordinal outcome we develop a test of mean difference and a non-
parametric test of Area Under the Curve (AUC). Additionally, multivariate
versions of these tests are proposed.

##  1 Introduction

Experimental testing is vital in the optimization of web applications, and as
such A/B testing has been widely adopted as a methodology for determining
optimal content for many web applications. A/B testing, from a statistical
perspective, is a randomized controlled trial (RCT), where users are
randomized to one of many user experiences. The goal of these trials is to
drive users into particular behaviors such as signing up to a service (a
conversion), or purchase of a product.

Despite the widespread adoption of A/B testing, the methodology used to make
decisions often does not fit the analyst’s needs or actions. Traditional
hypothesis tests require the collection of a sample with fixed size, with no
interim analysis and no post-hoc study continuation. However, in the context
of application optimization, data streams in over the course of days or weeks
and discussion making needs to be flexible and dynamic. For this reason,
sequentially valid hypothesis tests are required, where significance can be
evaluated at any point during the run of the experiment, and type I error is
controlled regardless of the number of interim analyses and length of the
experiment.

In Section 2 we review the mixture sequential likelihood ratio (mSPRT) and
define asymptotic equivalence. Section 3 provides results examining the
asymptotic behavior of mSPRT type ratios for sequences of maximum likelihood
or generalized method of moment estimators (GMM). Section 4 proposed an
efficient family of mixture distributions that can be calculated in closed
form. Section 5 applies the theoretical results of Section 4 to the risk
ratio, odds ratio and difference of proportions metrics for binary outcomes
and the difference of means and AUC metrics for non-binary outcomes.

##  2 Preliminaries: Sequential Likelihood Ratio Tests

Suppose that we wish to test the hypothesis $H_{0}:\theta_{1}=\theta_{0}$ for
some family of probability distributions $f(x\mid\theta)$, where $\theta$ is a
vector of parameters. In a sequential experiment, we observe a sequence of
observations from this distribution $X_{1},X_{2},...X_{\infty}$, and wish to
determine a stopping rule $T$, which may at any point in the sequence reject
the null hypothesis and terminate the experiment. The sequential likelihood
ratio of $\theta_{1}$ versus $\theta_{0}$ is defined as

$L_{n}=\prod_{i=1}^{n}\frac{f(X_{1},...,X_{n}\mid\theta_{1})}{f(X_{1},...,X_{n}\mid\theta_{0})}$

Analogous to the likelihood ratio test in classical statistics, Wald (1945)
suggested rejecting $H_{0}$ when the sequential likelihood ratio rose above a
certain level. Based on the fact that $L_{n}$ is a martingale, Wald (1945)
derived type I error probabilities for this test based on the identity

$P_{\theta_{0}}(\textrm{max}\ L_{n}\geq\alpha^{-1})\leq\alpha.$

This identity guarantees that rejecting the null hypothesis when the
likelihood ratio attains a value of $1/\alpha$ provides an error rate less
than or equal to $\alpha$.

For the majority of real world analyses, the hypotheses are composite rather
than simple, which complicates the problem considerably. Researchers have
considered a number of different generalizations to the Wald sequential
likelihood ratio, including the adaptive likelihood ratio (Robbins et al.,
1970) and the sequential generalized likelihood ratio (Schwarz et al., 1962)
(see Lai (2004) for a review). For this work we consider the mixture
likelihood ratio test (mSPRT) (Robbins, 1970), which averages the numerator
over a specified prior distribution ($g$) for $\theta$

$\Lambda_{n}^{g}=\frac{\int\prod_{i}^{n}f(X_{i}\mid\theta)g(\theta)d\theta}{\prod_{i}^{n}f(X_{i}\mid\theta_{0})}.$

The mixture likelihood ratio then rejects the null hypothesis when

$T(\alpha)=\Lambda_{n}^{g}>\frac{1}{\alpha},$

and terminates the experiment and rejects the null hypothesis at sample size

$\tau_{G}=\textrm{inf}\bigg{\\{}n\geq
1:\Lambda_{n}^{g}>\frac{1}{\alpha}\bigg{\\}}.$

Like $L_{n}$, $\Lambda_{n}^{g}$ is a martingale, and thus the mixture
likelihood ratio test is guaranteed to be a level alpha test, in that under
the null hypothesis the probability of terminating is

$P(\tau_{G}<\infty)\leq\alpha.$

Note that the inequality here is generally fairly tight as the “overshoot” of
sequential tests is typically small (Lai, 2004).

For composite null hypothesis, Wald (1945) suggested integrating over the
denominator as well as the numerator, yielding a likelihood ratio of the form

$\Psi_{n}^{g_{1}g_{2}}=\frac{\int\prod_{i}^{n}f(X_{i}\mid\theta)g_{2}(\theta)d\theta}{\int\prod_{i}^{n}f(X_{i}\mid\theta)g_{1}(\theta)d\theta},$

where $g_{2}$ is the distribution of the parameters under the alternate
hypothesis and $g_{1}$ is the distribution under the null hypothesis.
Rejecting the null when $\textrm{max}\ \Psi^{g_{1}g_{2}}_{n}\geq\alpha^{-1}$
limits the averaged type one error over the null distribution
$\int\alpha(\theta)g_{1}(\theta)d\theta<\alpha$, where $\alpha(\theta)$ is the
type I error for the particular paramter configuration $\theta$ (Lai, 2004).

For large scale experiments we are particularly interested in the large sample
behavior of the mSPRT and asymptotic approximations of it. Examination of the
large sample characteristics of sequential tests has a long history, starting
with Bartlett (1946) and Cox (1963) with early work summarized by Joanes
(1972). We now introduce two definitions relevant for large sample
approximations to the mSPRT

###### Definition 1

Let $\Lambda_{n}^{g}$ be the mSPRT testing the null hypothesis
$H_{0}:\theta=\theta^{0}$. A sequence $\Lambda_{n}^{\prime}$ is asymptotically
valid for the mSPRT $\Lambda_{n}$ if

$\Lambda_{n}^{\prime}\overset{d}{\Rightarrow}\Lambda_{n}^{g}$

under $H_{0}$.

###### Definition 2

Let $\Lambda_{n}^{g}$ be the mSPRT testing the null hypothesis
$H_{0}:\theta=\theta^{0}$. A sequence $\Lambda_{n}^{\prime}$ is asymptotically
equivalent to the mSPRT $\Lambda_{n}$ if

$\Lambda_{n}^{\prime}\overset{d}{\Rightarrow}\Lambda_{n}^{g}$

under any value of $\theta$.

Tests which are asymptotically equivalent to an mSPRT perform similarly at
large sample sizes to the underlying mSPRT, while asymptotically valid
approximations only perform similarly under the null hypothesis. Obviously,
asymptotically equivalent test statistics are also asymptotically valid.

Tests which are asymptotically valid do not necessarily keep the same type I
error characteristics of the mSPRT. At low sample sizes, the asymptotic test
may reject too much or not enough. There are two approaches that alleviate
this behavior. Firstly, it may be prudent to ignore rejections at low sample
sizes and only trigger a termination of the test once a threshold sample size
has been reached. Secondly, if $g$ is heavily concentrated around $\theta^{0}$
then rejections at low sample size are unlikely even under an approximation of
$\Lambda^{g}_{n}$, and thus inflated type I errors are of little practical
concern.

##  3 Sequences of Estimators

Let $\hat{\theta}_{n}(X_{1},...,X_{n})$ be a consistent estimator for the
parameter vector $\theta$, which conditional upon $X_{1},...,X_{n-1}$ is a
one-to-one function of $X_{n}$. Given a known sampling distribution for
$\hat{\theta}_{n}$, it is desirable to create a valid sequential hypothesis
test based on this distribution.

###### Theorem 1

Let $H_{0}:\theta=\theta_{0}$ be a statistical hypothesis and
$\hat{\theta}_{n}$ be sufficient for $\theta$, then the mSPRT is equal to

$\Lambda_{n}^{g}=\frac{\int
p(\hat{\theta}_{n}\mid\theta)g(\theta)d\theta}{p(\hat{\theta}_{n}\mid\theta_{0})}.$

Proof: Let $S_{n}=n\hat{\theta}_{n}-(n-1)\hat{\theta}_{n-1}$, then

$\displaystyle\prod_{i}^{n}f(X_{i}\mid\theta)$
$\displaystyle=f(X_{1},...,X_{n}\mid\theta)$
$\displaystyle=p(S_{1},...,S_{n}\mid\theta)$
$\displaystyle=p(\sum^{n}_{i}S_{i}\mid\theta)p(S_{1},...,S_{n}\mid\sum^{n}_{i}S_{i},\theta)$
$\displaystyle=p(\hat{\theta}_{n}\mid\theta)p(S_{1},...,S_{n}\mid\hat{\theta}_{n}).$

Because the estimator is sufficient for $\theta$, the second probability
cancels out in the likelihood ratio, leaving

$\displaystyle\Lambda_{n}^{g}$
$\displaystyle=\frac{\int\prod_{i}^{n}f(X_{i}\mid\theta)g(\theta)d\theta}{\prod_{i}^{n}f(X_{i}\mid\theta_{0})}$
$\displaystyle=\frac{\int
p(\hat{\theta}_{n}\mid\theta)g(\theta)d\theta}{p(\hat{\theta}_{n}\mid\theta_{0})}.$

$\blacksquare$

Theorem 1 allows for a significant reduction in computational complexity when
it can be applied, as there are often easy closed form approximations of the
distribution of $\hat{\theta}$.

###### Theorem 2

Let $H_{0}:\theta=\theta_{0}$ be a statistical hypothesis, and
$\hat{\theta}_{n}$ be the maximum likelihood estimators for $\theta$. Further,
let
$\sqrt{n}\hat{\Sigma}_{n}^{-\frac{1}{2}}(\hat{\theta}_{n}-\theta)\overset{d}{\Rightarrow}\mathcal{N}(0,I)$,
where $I$ is the identity matrix and $\hat{\Sigma}_{n}$ is a consistent
estimate of the limiting covariance when $H_{0}$ is true (
$\Sigma(\theta^{0})$ ).

$\frac{\int\phi(\hat{\theta}_{n}\mid\theta,n^{-1}\hat{\Sigma}_{n})g(\theta)d\theta}{\phi(\hat{\theta}_{n}\mid\theta^{0},n^{-1}\hat{\Sigma}_{n})}$

is asymptotically valid for $\Lambda_{n}^{g}$, and is asymptotically
equivalent if $\hat{\Sigma}_{n}$ is a consistent for all $\theta$.

Proof:

Let $\hat{\theta}_{n}$ represent the maximum likelihood estimators for $f$,
then the Taylor expansion around the MLE is

$\log f(\theta)=\log f(\hat{\theta})+(\theta-\hat{\theta}_{n})^{T}\nabla
f\big{|}_{\hat{\theta}_{n}}+\frac{1}{2}(\theta-\hat{\theta}_{n})^{T}\nabla^{2}f\big{|}_{\hat{\theta}_{n}}(\theta-\hat{\theta}_{n})+...,$

with higher order terms tending to 0 as $n$ increases. Using the fact that
$\nabla f\big{|}_{\hat{\theta}_{n}}=0$, we express the log likelihood ratio as

$\displaystyle\log f(\theta)-\log f(\theta^{0})$ $\displaystyle=\log
f(\theta)-\log f(\hat{\theta}_{n})-(\log f(\theta^{0})-\log
f(\hat{\theta}_{n}))$
$\displaystyle\Rightarrow-\frac{1}{2}(\theta-\hat{\theta}_{n})^{T}\nabla^{2}f\big{|}_{\hat{\theta}_{n}}(\theta-\hat{\theta}_{n})+\frac{1}{2}(\theta^{0}-\hat{\theta}_{n})^{T}\nabla^{2}f\big{|}_{\hat{\theta}_{n}}(\theta^{0}-\hat{\theta}_{n}).$

Exponentiating both sides yields

$\displaystyle\int\frac{f(X_{1},...,X_{n}|\theta)}{f(X_{1},...,X_{n}|\theta^{0})}g(\theta)d\theta$
$\displaystyle\Rightarrow\frac{\int\phi(\hat{\theta}_{n}\mid\theta,-(\nabla^{2}f\big{|}_{\hat{\theta}_{n}})^{-1})g(\theta)d\theta}{\phi(\hat{\theta}_{n}\mid\theta^{0},-(\nabla^{2}f\big{|}_{\hat{\theta}_{n}})^{-1})}$

If $\hat{\Sigma}_{n}\Rightarrow-n(\nabla^{2}f\big{|}_{\hat{\theta}_{n}})^{-1}$
equivalence is achieved by substituting $n^{-1}\hat{\Sigma}_{n}$ into the
ratio. If it is only consistent under the null hypothesis then the ratio is
valid for $\Lambda_{n}^{g}$.

$\blacksquare$

###### Theorem 3

Let $\theta=\begin{bmatrix}\beta\\\ \eta\end{bmatrix}$ be a partitioning of
the parameter space, $H_{0}:\beta=\beta^{0}$ be the hypothesis of interest and

$\Lambda_{n}=\frac{\int\phi(\hat{\theta}_{n}\mid\theta,n^{-1}\hat{\Sigma}_{n})g(\beta)d\beta}{\phi(\hat{\theta}_{n}\mid\theta^{0},n^{-1}\hat{\Sigma}_{n})}$

be an asymptotically valid mSPRT, where
$\hat{\Sigma}_{n}=\begin{bmatrix}\hat{\Sigma}^{\beta\beta}_{n}&\hat{\Sigma}^{\beta\eta}_{n}\\\
\hat{\Sigma}^{\eta\beta}_{n}&\hat{\Sigma}^{\eta\eta}_{n}\end{bmatrix}$ is a
consistent estimate of the limiting covariance matrix and
$\theta^{0}=\begin{bmatrix}\beta^{0}\\\ \eta\end{bmatrix}$. Further, let
$\hat{\Sigma}^{\eta\theta}=\hat{\Sigma}^{\theta\eta}=0$. Then

$\Lambda_{n}=\frac{\int\phi(\hat{\beta}_{n}\mid\beta,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})g(\beta)d\theta}{\phi(\hat{\beta}_{n}\mid\beta^{0},n^{-1}\hat{\Sigma}_{n}^{\beta\beta})}$

Proof:

$\displaystyle\frac{\int\phi(\hat{\theta}_{n}\mid\theta,n^{-1}\hat{\Sigma}_{n})g(\beta)d\beta}{\phi(\hat{\theta}_{n}\mid\theta^{0},n^{-1}\hat{\Sigma}_{n})}$
$\displaystyle=\frac{\int\phi(\hat{\eta}_{n}\mid\eta,n^{-1}\hat{\Sigma}_{n}^{\eta\eta})\phi(\hat{\beta}_{n}\mid\beta,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})g(\beta)d\theta}{\phi(\hat{\eta}_{n}\mid\eta,n^{-1}\hat{\Sigma}_{n}^{\eta\eta})\phi(\hat{\beta}_{n}\mid\beta^{0},n^{-1}\hat{\Sigma}_{n}^{\beta\beta})}$
$\displaystyle=\frac{\int\phi(\hat{\beta}_{n}\mid\beta,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})g(\beta)d\beta}{\phi(\hat{\beta}_{n}\mid\beta^{0},n^{-1}\hat{\Sigma}_{n}^{\beta\beta})}$

$\blacksquare$

In many cases, the assumption that $\hat{\Sigma}^{\theta\eta}=0$ is not met,
making inference in the case of nuisance parameters difficult. Let $h$ be the
true prior distribution of $\eta$. Averaging the null hypothesis over $h$
yields a likelihood ratio of

$\Psi_{n}^{h(gh)}=\frac{\int\int\prod_{i}^{n}f(X_{i}\mid\beta,\eta)g(\beta)h(\eta)d\eta
d\beta}{\int\prod_{i}^{n}f(X_{i}\mid\beta_{0},\eta)h(\eta)d\eta}.$

Using this ratio as the test statistic limits type I error averaged across all
experiments to the specified level. Of course $h$ is unknown, making direct
use of $\Psi$ problematic. Fortunately, the following theorem shows that
asymptotically, the ratio does not depend on $h$.

###### Theorem 4

Let $\theta=\begin{bmatrix}\beta\\\ \eta\end{bmatrix}$ be a partitioning of
the parameter space, $H_{0}:\beta=\beta^{0}$ be the hypothesis of interest and
$\hat{\theta}_{n}$ be the maximum likelihood estimators for $\theta$. Further,
let
$\sqrt{n}\hat{\Sigma}_{n}^{-\frac{1}{2}}(\hat{\theta}_{n}-\theta)\overset{d}{\Rightarrow}\mathcal{N}(0,I)$,
where $I$ is the identity matrix and
$\hat{\Sigma}_{n}=\begin{bmatrix}\hat{\Sigma}^{\beta\beta}_{n}&\hat{\Sigma}^{\beta\eta}_{n}\\\
\hat{\Sigma}^{\eta\beta}_{n}&\hat{\Sigma}^{\eta\eta}_{n}\end{bmatrix}$ is a
consistent estimate of the limiting covariance $\Sigma(\theta)$.

(1) If $h(\eta)\propto 1$, then

$\frac{\int\phi(\hat{\beta}_{n}\mid\beta,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})g(\beta)d\theta}{\phi(\hat{\beta}_{n}\mid\beta^{0},n^{-1}\hat{\Sigma}_{n}^{\beta\beta})}$

is asymptotically equivalent to $\Psi_{n}^{h(gh)}$.

(2) If
$h^{\prime}(\eta|\beta)=h(\eta+\Sigma^{\eta^{0}\beta^{0}}(\Sigma^{\beta^{0}\beta^{0}})^{-1}(\beta_{0}-\beta))$
then

$\frac{\int\phi(\hat{\beta}_{n}\mid\beta,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})g(\beta)d\theta}{\phi(\hat{\beta}_{n}\mid\beta^{0},n^{-1}\hat{\Sigma}_{n}^{\beta\beta})}$

is asymptotically valid for $\Psi_{n}^{h(gh^{\prime})}$ at
$\theta^{0}=\begin{bmatrix}\beta^{0}\\\ \eta^{0}\end{bmatrix}$.

Proof: As with Theorem 2, we use the fact that

$\log f(\theta)-\log
f(\hat{\theta})\Rightarrow\frac{1}{2}(\theta-\hat{\theta}_{n})^{T}\nabla^{2}f\big{|}_{\hat{\theta}_{n}}(\theta-\hat{\theta}_{n}),$

and thus

$\displaystyle\Psi_{n}^{h(gh)}$
$\displaystyle\Rightarrow\int\frac{\int\phi(\hat{\theta}_{n}\mid\theta,n^{-1}\hat{\Sigma}_{n})h(\eta)d\eta}{\int\phi(\hat{\theta}_{n}\mid\theta^{0},n^{-1}\hat{\Sigma}_{n})h(\eta)d\eta}g(\beta)d\beta$
$\displaystyle=\int\frac{\phi(\hat{\beta}_{n}\mid\beta,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})\int\phi(\hat{\eta}_{n}\mid\eta+\hat{\Sigma}_{n}^{\eta\beta}(\hat{\Sigma}_{n}^{\beta\beta})^{-1}(\hat{\beta}_{n}-\beta),V_{n})h(\eta)d\eta}{\phi(\hat{\beta}_{n}\mid\beta^{0},n^{-1}\hat{\Sigma}_{n}^{\beta\beta})\int\phi(\hat{\eta}_{n}\mid\eta+\hat{\Sigma}_{n}^{\eta\beta}(\hat{\Sigma}_{n}^{\beta\beta})^{-1}(\hat{\beta}_{n}-\beta^{0}),V_{n})h(\eta)d\eta}g(\beta)d\beta$

where
$V_{n}=n^{-1}\hat{\Sigma}_{n}^{\eta\eta}-n^{-2}\hat{\Sigma}_{n}^{\beta\eta}(\hat{\Sigma}_{n}^{\eta\eta})^{-1}\hat{\Sigma}_{n}^{\eta\beta}$.
Using the fact that $h$ is uniform, the inner integrals cancel out

$\displaystyle\Psi_{n}^{h(gh)}$
$\displaystyle\Rightarrow\int\frac{\phi(\hat{\beta}_{n}\mid\beta,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})\int\phi(\hat{\eta}_{n}\mid\eta+\hat{\Sigma}_{n}^{\eta\beta}(\hat{\Sigma}_{n}^{\beta\beta})^{-1}(\hat{\beta}_{n}-\beta),V_{n})d\eta}{\phi(\hat{\beta}_{n}\mid\beta^{0},n^{-1}\hat{\Sigma}_{n}^{\beta\beta})\int\phi(\hat{\eta}_{n}\mid\eta+\hat{\Sigma}_{n}^{\eta\beta}(\hat{\Sigma}_{n}^{\beta\beta})^{-1}(\hat{\beta}_{n}-\beta^{0}),V_{n})d\eta}g(\beta)d\beta$
$\displaystyle=\int\frac{\phi(\hat{\beta}_{n}\mid\beta,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})}{\phi(\hat{\beta}_{n}\mid\beta^{0},n^{-1}\hat{\Sigma}_{n}^{\beta\beta})}g(\beta)d\beta$

This proves the first part of the proposition. For the second proposition, we
have that

$\displaystyle\Psi_{n}^{h(gh^{\prime})}$
$\displaystyle\Rightarrow\int\frac{\phi(\hat{\beta}_{n}\mid\beta,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})\int\phi(\hat{\eta}_{n}\mid\eta+\hat{\Sigma}_{n}^{\eta\beta}(\hat{\Sigma}_{n}^{\beta\beta})^{-1}(\hat{\beta}_{n}-\beta),V_{n})h(\eta+\Sigma^{\eta^{0}\beta^{0}}(\Sigma^{\beta^{0}\beta^{0}})^{-1}(\beta_{0}-\beta))d\eta}{\phi(\hat{\beta}_{n}\mid\beta^{0},n^{-1}\hat{\Sigma}_{n}^{\beta\beta})\int\phi(\hat{\eta}_{n}\mid\eta+\hat{\Sigma}_{n}^{\eta\beta}(\hat{\Sigma}_{n}^{\beta\beta})^{-1}(\hat{\beta}_{n}-\beta^{0}),V_{n})h(\eta)d\eta}g(\beta)d\beta$

Consider the inner integrals of the numerator

$\phi(\hat{\eta}_{n}\mid\eta+\hat{\Sigma}_{n}^{\eta\beta}(\hat{\Sigma}_{n}^{\beta\beta})^{-1}(\hat{\beta}_{n}-\beta),V_{n})h(\eta+\Sigma^{\eta^{0}\beta^{0}}(\Sigma^{\beta^{0}\beta^{0}})^{-1}(\beta_{0}-\beta))$

and denominator

$\phi(\hat{\eta}_{n}\mid\eta+\hat{\Sigma}_{n}^{\eta\beta}(\hat{\Sigma}_{n}^{\beta\beta})^{-1}(\hat{\beta}_{n}-\beta^{0}),V_{n})h(\eta)$

define $\bar{\eta}_{1}(\beta)$ to be the maximizer of the numerator and
$\bar{\eta}_{0}(\beta)$ the maximizer of the denominator term. Let
$\eta^{*}(\beta)=\hat{\eta}_{n}-\hat{\Sigma}_{n}^{\eta\beta}(\hat{\Sigma}_{n}^{\beta\beta})^{-1}(\hat{\beta}_{n}-\beta)$
to be the maximum likelihood solution to
$\phi(\hat{\eta}_{n}\mid\eta+\hat{\Sigma}_{n}^{\eta\beta}(\hat{\Sigma}_{n}^{\beta\beta})^{-1}(\hat{\beta}_{n}-\beta),V_{n})$.
Applying the Laplace approximation to the numerator and denominator yields.

$\displaystyle\int\frac{\phi(\hat{\beta}_{n}\mid\beta,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})\phi(\hat{\eta}_{n}\mid\bar{\eta}_{1}(\beta)+\hat{\Sigma}_{n}^{\eta\beta}(\hat{\Sigma}_{n}^{\beta\beta})^{-1}(\hat{\beta}_{n}-\beta),V_{n})h(\bar{\eta}_{1}(\beta)+\Sigma^{\eta^{0}\beta^{0}}(\Sigma^{\beta^{0}\beta^{0}})^{-1}(\beta_{0}-\beta))|\Omega^{1}|^{\frac{1}{2}}}{\phi(\hat{\beta}_{n}\mid\beta^{0},n^{-1}\hat{\Sigma}_{n}^{\beta\beta})\phi(\hat{\eta}_{n}\mid\bar{\eta}_{0}(\beta^{0})+\hat{\Sigma}_{n}^{\eta\beta}(\hat{\Sigma}_{n}^{\beta\beta})^{-1}(\hat{\beta}_{n}-\beta^{0}),V_{n})h(\bar{\eta}_{0}(\beta^{0}))|\Omega^{0}|^{\frac{1}{2}}}g(\beta)d\beta,$

where the $\Omega^{i}$ are the inverted second derivative matrices at
$\bar{\eta}^{i}$. Because
$\sqrt{n}(\bar{\eta}(\beta)-\eta^{*}(\beta))\Rightarrow 0$ (Ghosh et al.,
2007) we may substitute in $\eta^{*}$ resulting in

$\displaystyle\Psi_{n}^{h(gh^{\prime})}$
$\displaystyle\Rightarrow\int\frac{\phi(\hat{\beta}_{n}\mid\beta,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})\phi(\hat{\eta}_{n}\mid\hat{\eta}_{n},V_{n})h(\hat{\eta}_{n}-\hat{\Sigma}_{n}^{\eta\beta}(\hat{\Sigma}_{n}^{\beta\beta})^{-1}(\hat{\beta}_{n}-\beta)+\Sigma^{\eta^{0}\beta^{0}}(\Sigma^{\beta^{0}\beta^{0}})^{-1}(\beta_{0}-\beta))|\Omega^{1}|^{\frac{1}{2}}}{\phi(\hat{\beta}_{n}\mid\beta^{0},n^{-1}\hat{\Sigma}_{n}^{\beta\beta})\phi(\hat{\eta}_{n}\mid\hat{\eta}_{n},V_{n})h(\hat{\eta}_{n}-\hat{\Sigma}_{n}^{\eta\beta}(\hat{\Sigma}_{n}^{\beta\beta})^{-1}(\hat{\beta}_{n}-\beta^{0}))|\Omega^{0}|^{\frac{1}{2}}}g(\beta)d\beta$
$\displaystyle=\int\frac{\phi(\hat{\beta}_{n}\mid\beta,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})h(\hat{\eta}_{n}-\hat{\Sigma}_{n}^{\eta\beta}(\hat{\Sigma}_{n}^{\beta\beta})^{-1}(\hat{\beta}_{n}-\beta)+\Sigma^{\eta^{0}\beta^{0}}(\Sigma^{\beta^{0}\beta^{0}})^{-1}(\beta_{0}-\beta))|\Omega^{1}|^{\frac{1}{2}}}{\phi(\hat{\beta}_{n}\mid\beta^{0},n^{-1}\hat{\Sigma}_{n}^{\beta\beta})h(\hat{\eta}_{n}-\hat{\Sigma}_{n}^{\eta\beta}(\hat{\Sigma}_{n}^{\beta\beta})^{-1}(\hat{\beta}_{n}-\beta^{0}))|\Omega^{0}|^{\frac{1}{2}}}g(\beta)d\beta.$

Under the null hypothesis of $\theta^{0}=\begin{bmatrix}\beta^{0}\\\
\eta^{0}\end{bmatrix}$, $\hat{\beta}_{n}\rightarrow\beta^{0}$ and
$\hat{\Sigma}\rightarrow\Sigma(\theta^{0})$ so we may substitute these in
within $h$.

$\Psi_{n}^{h(gh^{\prime})}\Rightarrow\int\frac{\phi(\hat{\beta}_{n}\mid\beta,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})h(\hat{\eta}_{n})|\Omega^{1}|^{\frac{1}{2}}}{\phi(\hat{\beta}_{n}\mid\beta^{0},n^{-1}\hat{\Sigma}_{n}^{\beta\beta})h(\hat{\eta}_{n})|\Omega^{0}|^{\frac{1}{2}}}g(\beta)d\beta.$

Similarly, under the null hypothesis of $\theta^{0}$,
$\Omega^{1}\rightarrow\Omega^{0}$. This can be seen by noting that the
quantities inside the inner integrals of the numerator and denominator are
asymptotically equivalent at $\eta^{*}(\beta)$ and $\eta^{*}(\beta^{0})$
respectively. Thus we then have the result that

$\Psi_{n}^{h(gh^{\prime})}\Rightarrow\int\frac{\phi(\hat{\beta}_{n}\mid\beta,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})}{\phi(\hat{\beta}_{n}\mid\beta^{0},n^{-1}\hat{\Sigma}_{n}^{\beta\beta})}g(\beta)d\beta.$

$\blacksquare$

###### Theorem 5

Let $\hat{\theta}_{nk}$ be a (sub-)set of generalized method of moments
estimators with limiting distribution
$\sqrt{nk}\hat{\Sigma}_{nk}^{-\frac{1}{2}}(\hat{\theta}_{nk}-\theta)\Rightarrow
N(0,I)$. The likelihood ratio

$\Lambda_{n}=\frac{\int\phi(\hat{\theta}_{nk}\mid\theta,(nk)^{-1}\hat{\Sigma})g(\theta)d\theta}{\phi(\hat{\theta}_{nk}\mid\theta_{0},(nk)^{-1}\hat{\Sigma})},$

is an asymptotically equivalent mSPRT as $n$ and $k>0$ become large.

Proof:

The sample moment conditions for $\hat{\theta}_{nk}$ are

$m_{n}(\theta)=\frac{1}{n}\sum_{i}^{n}\gamma(X_{i},\theta),$

and standard GMM asymptotic theory finds that

$\sqrt{nk}\Gamma^{-1}(\hat{\theta}_{nk}-\theta)\Rightarrow\sqrt{nk}m_{nk}(\theta),$

where $\Gamma$ is a constant matrix, and $m_{nk}(\theta)\Rightarrow
N(0,nk\Omega)$ by the central limit theorem. Consider observing the sequence

$S_{n}=\Gamma\sum_{i=(n-1)k+1}^{nk}\gamma(X_{i},\theta)+k\theta,$

then by construction

$\sqrt{nk}\hat{\theta}_{nk}\Rightarrow\frac{\sum_{i}^{n}S_{n}}{\sqrt{nk}}.$

By the central limit theorem, for large $k$,

$S_{n}\Rightarrow N(k\theta,k\Gamma\Omega\Gamma^{T}).$

The sample mean and sample covariance are the maximum likelihood estimates of
$k\theta$ and $k\Gamma\Omega\Gamma^{T}$, so we may apply Theorem 2 to obtain
an asymptotically valid mSPRT. Because the sample mean and sample covariance
are independent, we may then apply Theorem 3 setting $\beta=k\theta$ and
$\eta=k\Gamma\Omega\Gamma^{T}$. This yields an mSPRT of

$\Lambda_{n}=\frac{\int\phi(\sum S_{i}\mid
nk\theta,nk\hat{\Sigma})g(\theta)d\theta}{\phi(\sum S_{i}\mid
nk\theta_{0},nk\hat{\Sigma})}.$

Replacing the $\sum S_{i}$ by the GMM estimator we arrive at

$\Lambda_{n}=\frac{\int\phi(\hat{\theta}_{nk}\mid\theta,(nk)^{-1}\hat{\Sigma})g(\theta)d\theta}{\phi(\hat{\theta}_{nk}\mid\theta_{0},(nk)^{-1}\hat{\Sigma})}.$

$\blacksquare$

By Theorem 5 we can still construct valid mSPRT tests for a wide variety of
estimators by restricting our interim-analyses to occur every $k$
observations, provided $k$ is large enough. The result relies on the
convergence of the central limit theorem within each batch of size $k$. In the
case of online experiments, where tens of thousands of observations stream
into a test per day, this it is not an unreasonable limitation. While a fixed
$k$ interval is posited, the size of the interval may be allowed to grow or
shrink over the course of the sequential test provided that each interval is
long enough for the central limit theorem to apply to the batch.

###### Corollary 1

Let $\hat{\theta}_{nk}$ be a (sub-)set of maximum likelihood estimators with
limiting distribution
$\sqrt{nk}\hat{\Sigma}_{nk}^{-\frac{1}{2}}(\hat{\theta}_{nk}-\theta)\Rightarrow
N(0,I)$. The likelihood ratio

$\Lambda_{n}=\frac{\int\phi(\hat{\theta}_{nk}\mid\theta,(nk)^{-1}\hat{\Sigma})g(\theta)d\theta}{\phi(\hat{\theta}_{nk}\mid\theta_{0},(nk)^{-1}\hat{\Sigma})},$

is an asymptotically equivalent mSPRT as $n$ and $k>0$ become large.

Corollary 1 follows immediately from the fact that maximum likelihood
estimators are also generalized method of moments estimators. Thus we now have
three ways to justify the use of the asymptotic marginal distribution of a
subset of the MLEs. Firstly, if the subset of estimators is independent of the
rest of the estimators, then Theorem 3 may be used to remove the nuisance
parameters. Otherwise, Theorem 4 shows that mSPRT limits to $L_{n}^{hg}$,
which controls the averaged type I error of the test. Finally, Corollary 1
shows that (non-averaged) type I error is controlled so long as the mSPRT is
evaluated at intervals. The fact that the mSPRT is valid when evaluated at
intervals suggests that type I errors are likely to be well controlled even
when $k=1$, as the deviation of the test statistic within each interval is
bounded.

###### Theorem 6

Let $Y_{1},...,Y_{\infty}$ be an ancillary sequence of random variables such
that

$p(X_{1},...,X_{n},Y_{1},...,Y_{n})=\prod_{i=1}^{n}p(X_{i}\mid\theta,Y_{i})p(Y_{i}|Y_{1},...,Y_{i-1}).$

then

$\Lambda_{n}^{g}=\frac{\int
p(X_{1},...,X_{n}\mid\theta,Y_{1},...,Y_{n})g(\theta)d\theta}{p(X_{1},...,X_{n}\mid\theta_{0},Y_{1},...,Y_{n})},$

Proof:

$\displaystyle\Lambda_{n}^{g}$
$\displaystyle=\frac{\int\prod_{i=1}^{n}p(X_{i}\mid\theta,Y_{i})p(Y_{i}|Y_{1},...,Y_{i-1})g(\theta)d\theta}{\prod_{i=1}^{n}p(X_{i}\mid\theta_{0},Y_{i})p(Y_{i}|Y_{1},...,Y_{i-1})}$
$\displaystyle=\frac{\int\prod_{i=1}^{n}p(X_{i}\mid\theta,Y_{i})g(\theta)d\theta}{\prod_{i=1}^{n}p(X_{i}\mid\theta_{0},Y_{i})}$
$\displaystyle=\frac{\int
p(X_{1},...,X_{n}\mid\theta,Y_{1},...,Y_{n})g(\theta)d\theta}{p(X_{1},...,X_{n}\mid\theta_{0},Y_{1},...,Y_{n})}.$

$\blacksquare$

##  4 An Efficient Family of Priors

Thus far we have avoided putting a defined functional form on the prior under
the alternative hypothesis $g$. While in principle any distribution may be
selected, computing the mixture integral numerically for each observation can
be prohibitively computationally expensive, especially in the case of online
experiments, where the number of observations is typically greater than
10,000, and may scale up to the millions. Fortunately, we may specify a family
of distributions that is both flexible enough to approximate any arbitrary
distribution while providing a closed form solution to the integral.

Let $g^{\prime}$ be a multivariate mixture normal density with $r$ components

$g^{\prime}(\theta)=\sum_{i=1}^{r}\phi(\theta|\mu_{i},\Upsilon_{i})w_{i},$

where $w_{i}$ is the probability of selecting the $i$th component. Let
$\Lambda_{n}$ be a mSPRT with normal likelihoods of the form

$\frac{\int\phi(\hat{\theta}_{n}\mid\theta,\hat{\Sigma}_{n})g(\theta)d\theta}{\phi(\hat{\theta}_{k}\mid\theta_{0},\hat{\Sigma}_{n})}.$

Given $g=g^{\prime}$, the mSPRT simplifies to

$\frac{\sum_{i}^{r}\phi(\hat{\theta}_{n}\mid\mu_{i},\hat{\Sigma}_{n}+\Upsilon_{i})w_{i}}{\phi(\hat{\theta}_{k}\mid\theta_{0},\hat{\Sigma}_{n})},$

removing the need for numeric integration. The mixture normal distribution
allows us to efficiently represent most distributions with just a few
components and also allows us the flexibility to model any continuous
distribution by simply increasing the number of terms.

##  5 Applications in Online Testing

###  5.1 $m$-Sample Tests of a Binomial Outcome

One of the most common goals for an online A/B test is to determine whether
one arm of the trial leads to more “conversions” then the others. A conversion
might indicate signing up for a news letter, a purchase, clicking on an Ad, or
any other positive action by the user. The outcome is therefore a Bernoulli
random variable $X_{i}\sim\textrm{Ber}(p_{Y_{i}})$, where
$Y_{i}\in\\{1,...,m\\}$ is the arm assigned to the $i$th individual. The
maximum likelihood estimators of $p$ are simply the sample proportions within
each group $\hat{p}_{j}=\frac{1}{n_{j}}\sum_{i}X_{i}\mathds{1}(Y_{i}=j)$,
where $\mathds{1}$ is the indicator function and
$n_{j}=\sum_{i}\mathds{1}(Y_{i}=j)$.

####  5.1.1 Risk Ratio

Since $\log$ is a one-to-one function, the estimators $\log(\hat{p})$ are the
maximum likelihood estimators of $\log(p)$ and by the delta method are
independent asymptotically normal variables with variance

$v_{j}=\hat{\textrm{var}}(\log(\hat{p}_{j}))=\frac{1-\hat{p}_{j}}{\hat{p}_{j}n_{j}}.$

Defining $\beta^{i}=\log(p_{i+1})-\log(p_{1})$ to be the log risk ratio versus
baseline, $\eta=\sum_{j}\log(p_{j})$ to be a nuisance parameter,
$\hat{\theta}_{n}=\begin{bmatrix}\hat{\eta}_{n}\\\
\hat{\beta}_{n}\end{bmatrix}$ to be the maximum likelihood estimates of the
re-parameterized model obtained by replacing $p$ by $\hat{p}$. The covariance
between $\beta^{i}$ and $\eta$ is $v_{i}-v_{1}$. This covariance tends to 0
under the the null hypothesis of $H_{0}:\beta=0$ if the allocation rates are
balanced across all variants. Therefore a consistent estimate of the
covariance under $H_{0}$ and balanced allocation is

$n^{-1}\hat{\Sigma}_{n}=\begin{bmatrix}\sum_{j}v_{j}&0\\\
0&n^{-1}\hat{\Sigma}^{\beta\beta}_{n}\end{bmatrix},$

where $n^{-1}(\hat{\Sigma}^{\beta\beta}_{n})_{jj}=v_{j+1}+v_{1}$ and
$n^{-1}(\hat{\Sigma}^{\beta\beta}_{n})_{ij}=v_{1}\ \ \forall i\neq j$. We may
now apply Theorems 2 and 3 to obtain an asymptotically valid mSPRT

$\frac{\int\phi(\hat{\beta}_{n}\mid\beta,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})g(\beta)d\beta}{\phi(\hat{\beta}_{n}\mid
0,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})}.$

Because $\hat{\beta}^{i}_{n}$ are maximum likelihood estimators for the
reparametrized model, when allocations are unbalanced or adaptive allocation,
Theorem 4 and Corollary 1 may be used to justify the mSPRT.

####  5.1.2 Odds Ratio

Derivation of an mSPRT for the log odds ratio proceeds similarly to the risk
ratio case. Setting
$\beta^{i}=\log(\frac{p_{i+1}}{1-p_{i+1}})-\log(\frac{p_{1}}{1-p_{1}})$ to be
the log odds ratio versus baseline, $\eta=\sum_{j}\log(\frac{p_{j}}{1-p_{j}})$
to be a nuisance parameter we apply Theorems 2 and 3 to obtain an
asymptotically valid mSPRT in the case of balanced allocations and utilize
Theorem 3 and Corollary 1 to justify it in the case of unbalanced or adaptive
allocations.

$\frac{\int\phi(\hat{\beta}_{n}\mid\beta,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})g(\beta)d\beta}{\phi(\hat{\beta}_{n}\mid
0,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})}.$

where $n^{-1}(\hat{\Sigma}^{\beta\beta}_{n})_{jj}=v_{j+1}+v_{1}$ and
$n^{-1}(\hat{\Sigma}^{\beta\beta}_{n})_{ij}=v_{1}\ \ \forall i\neq j$ for
$v_{i}=\frac{1}{n_{i}\hat{p}_{i}}+\frac{1}{n_{i}(1-\hat{p}_{i})}$.

####  5.1.3 Difference in Proportions

We begin our mSPRT for difference in proportions similarly to the risk and
odds ratio cases, by defining $\beta^{i}=p_{i+1}-p_{1}$ and
$\eta=\sum_{i}p_{i}$. By an identical argument,

$\frac{\int\phi(\hat{\beta}_{n}\mid\beta,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})g(\beta)d\beta}{\phi(\hat{\beta}_{n}\mid
0,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})}$

is an asymptotically valid mSPRT under balanced allocation, where
$n^{-1}(\hat{\Sigma}^{\beta\beta}_{n})_{jj}=v_{j+1}+v_{1}$ and
$n^{-1}(\hat{\Sigma}^{\beta\beta}_{n})_{ij}=v_{1}\ \ \forall i\neq j$ for
$v_{i}=\frac{\hat{p}_{i}(1-\hat{p}_{i})}{n_{i}}$. Again, we utilize Theorem 4
and Corollary 1 to justify it in the case of unbalanced or adaptive
allocations.

When there are only two arms to the trial, this reduces to the difference in
proportions test described in Johari et al. (2015); However, the derivation of
the test in that paper is incorrect. Lemma 4 of Johari et al. (2015) first
asserts that $\hat{p}_{2}-\hat{p}_{1}$ is asymptotically independent of
$\hat{p}_{2}+\hat{p}_{1}$ for any $p$, whereas this is only true under the
null hypothesis of equality of proportions and equal allocation rates since
$\textrm{cov}(\hat{p}_{2}-\hat{p}_{1},\hat{p}_{2}+\hat{p}_{1})=v_{2}-v_{1}$.
Secondly, they assert that $\hat{\Sigma}^{\beta\beta}_{n}$ converges to both
$\Sigma^{\beta\beta}(\beta=0)$ and $\Sigma^{\beta\beta}(\beta)\ \
\forall\beta\neq 0$ in the same experiment, which of course it can not do
because nothing can converge to two different matrices.

Another problem with difference in proportion is that it is a poor measure of
the effect of an intervention. For example, increasing the conversion rate of
a page from 50% to 51% is a small change, likely to lead to a small increase
in the underlying profit. A change from a 1% conversion rate to a 2%
conversion rate on the other hand is a huge effect, which (if revenue is
linked to conversions) could double profit. For this reason, difference in
proportion is generally of less business use than risk ratio.

Because difference in proportion effect size varies based on the base rate,
identifying a $g$ that is useful for both high and low conversion experiments
is a challenge. Our approach is to model the prior of the Cohen’s D effect
size (Hedges and Olkin, 2014) instead of the raw proportion. Effect size is
defined as

$\frac{p_{j}-p_{1}}{\bar{\sigma}},$

where $\bar{\sigma}$ is a measure of the underlying variation, and could be
the baseline standard deviation $\sqrt{p_{1}(1-p_{1})}$, the average variation
across groups $\sqrt{\frac{1}{m}\sum_{i}p_{i}(1-p_{i})}$, or any other
population quantity that can be estimated consistently. Then

$\frac{\int\phi(\hat{\beta}_{n}\mid\beta,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})g(\frac{\beta}{\hat{\sigma}_{n}})d\beta}{\phi(\hat{\beta}_{n}\mid
0,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})}\Rightarrow\frac{\int\phi(\hat{\beta}_{n}\mid\beta,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})g(\frac{\beta}{\bar{\sigma}})d\beta}{\phi(\hat{\beta}_{n}\mid
0,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})},$

where $\hat{\sigma}_{n}$ is a consistent approximation of $\bar{\sigma}$ for
example $\sqrt{\hat{p}_{1}(1-\hat{p}_{1})}$ or
$\sqrt{\frac{1}{m}\sum_{i}\hat{p}_{i}(1-\hat{p}_{i})}$.

###  5.2 $m$-Sample Tests of a Numeric Outcome

For numeric outcomes, instead of being distributed binomially, the outcome is
distributed according to $X_{i}\sim f_{Y_{i}}$, where $f_{Y_{i}}$ is the
distribution of the $i$th arm of the study. Of particular interest is the case
when each group is normally distributed
$f_{Y_{i}}(X_{i})=\phi(X_{i}|\mu_{Y_{i}},\sigma^{2}_{Y_{i}})$.

####  5.2.1 Difference in Means

Let us begin by considering the simplified case of
$f_{Y_{i}}(X_{i})=\phi(X_{i}|\mu_{Y_{i}},\sigma^{2}_{Y_{i}})$. The maximum
likelihood estimates are the group means and standard deviations
$\hat{\mu}_{j}=\frac{1}{n_{j}}\sum_{i}X_{i}\mathds{1}(Y_{i}=j)$ and
$\hat{\sigma}^{2}_{j}=\frac{1}{n_{j}-1}\sum_{i}(X_{i}-\hat{\mu}_{j})\mathds{1}(Y_{i}=j)$.
We begin by reparameterizing the model as $\beta^{i}=\mu_{i+1}-\mu_{1}$,
$\eta=(\sum_{i}\mu_{i},\sigma^{2}_{1},...,\sigma^{2}_{n})$.

The MLEs for $\mu$ and $\sigma^{2}$ are the sample means and variances, which
under the normal model are independent. Under the null hypothesis of equality
of means, if the allocations are balance and $\sigma^{2}_{i}$ are all equal
then the asymptotic covariance between $\hat{\beta}_{n}$ and $\hat{\eta}_{n}$
are all zero. Therefore we may apply Theorems 2 and 3 to obtain

$\frac{\int\phi(\hat{\beta}_{n}\mid\beta,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})g(\beta)d\beta}{\phi(\hat{\beta}_{n}\mid
0,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})},$

where $n^{-1}(\hat{\Sigma}^{\beta\beta}_{n})_{jj}=v_{j+1}+v_{1}$ and
$n^{-1}(\hat{\Sigma}^{\beta\beta}_{n})_{ij}=v_{1}\ \ \forall i\neq j$ for
$v_{i}=\frac{\hat{\sigma}^{2}_{i}}{n_{i}}$.

Because $\hat{\beta}^{i}_{n}$ are maximum likelihood estimators for the
reparametrized normal model, when allocations are unbalanced and/or the true
variances are unequal, Theorem 4 and Corollary 1 may be used to justify the
mSPRT.

In reality, very few outcomes in online testing applications are even
approximately normally distributed. Many have heavy tails and high skew.
Fortunately, Theorem 5 provides a foundation for inference in the case of non-
normal data. Firstly note that from the central limit theorem we have that
$\hat{\Sigma}^{\beta\beta}_{n}$ remains a consistent estimator of the
asymptotic covariance. Let $Y_{n}$ be distributed multinomial with
probabilities $q$, and consider

$\gamma(X_{i},Y_{i},\mu_{j}-\mu_{1})=\frac{X_{i}\mathds{1}(Y_{i}=j)}{q_{j}}-\frac{X_{i}\mathds{1}(Y_{i}=1)}{q_{1}}-(\mu_{j}-\mu_{1}).$

The expectation of $\gamma$ is 0, so

$\frac{1}{n}\sum_{i}\gamma(X_{i},Y_{i},\mu_{j}-\mu_{1})=0$

is a (generalized) method of moments estimator with solution

$(\hat{\beta}^{j}_{n})_{\textrm{gmm}}=\frac{1}{n}\sum_{i}\big{(}\frac{X_{i}\mathds{1}(Y_{i}=j)}{q_{j}}-\frac{X_{i}\mathds{1}(Y_{i}=1)}{q_{1}}\big{)}=\mu_{j}-\mu_{1}.$

Because $n_{j}\Rightarrow nq_{j}$, we have that
$\sqrt{n}\hat{\beta}^{j}_{n}\Rightarrow\sqrt{n}(\hat{\beta}^{j}_{n})_{\textrm{gmm}}$,
and so the maximum likelihood estimators under the normal model are
asymptotically equivalent to GMM estimators under the general distribution
model. Thus, we may use Theorem 5 as a basis for inference.

As with differences in proportions, the scale of $X$ is an important factor in
determining the likely prior distribution effect sizes and again the solution
is to model the prior effect size using Cohen’s D type statistics

$\frac{\mu_{j}-\mu_{1}}{\bar{\sigma}},$

where $\bar{\sigma}$ here equals the population quantity that may be estimated
consistently; for example, the variance of the baseline variant
($\textrm{var}_{f_{1}}(X)$) or the average variance across variants
($\frac{1}{m}\sum_{j}^{m}\textrm{var}_{f_{j}}(X)$). Then

$\frac{\int\phi(\hat{\beta}_{n}\mid\beta,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})g(\frac{\beta}{\hat{\sigma}_{n}})d\beta}{\phi(\hat{\beta}_{n}\mid
0,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})}\Rightarrow\frac{\int\phi(\hat{\beta}_{n}\mid\beta,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})g(\frac{\beta}{\bar{\sigma}})d\beta}{\phi(\hat{\beta}_{n}\mid
0,n^{-1}\hat{\Sigma}_{n}^{\beta\beta})},$

where $\hat{\sigma}$ is the sample estimate of $\bar{\sigma}$, for example
$\hat{\sigma}_{n}=\sqrt{\frac{1}{n_{1}-1}\sum_{i}(X_{i}-\hat{\mu}_{n})^{2}\mathds{1}(Y_{i}=1)}$.

###  5.3 Non-parametric Superiority

As mentioned in the previous section. Many applications in online testing
involve outcomes with heavy tails and high skew. The presence of outliers due
either to exceptional users, or data collection errors may factor into
choosing appropriate methodologies for analysis.

For heavy tailed distributions with outliers, it is well known that the mean,
as a measure of central tendency, is a questionable choice. The mean is
heavily influenced by the tail behavior of a distribution, and thus any test
based on mean differences will require large sample sizes to reach
significance. Further, the result of that test may be dominated by the
behavior of a minority of exceptional users rather than representing the
effects of the experiment on the majority.

Another important use case is ordinal data, which, while ordered, does not
have an intrinsic unit of measurement. Examples from online testing might be
the number of steps a user took through the registration process, or a user
selected product rating from “Very good” to “Very Poor.” Using means to
measure the central tendency of an ordinal variable imposes an arbitrary unit
on the variable, which may not be appropriate.

Addressing both the continuous and ordinal case in the comparison of two
samples is known as the Nonparametric Behrens-Fisher Problem (Brunner and
Munzel, 2000). In the non-sequential context Brunner and Munzel (2000)
developed a two-sample test that shows good small sample characteristics.

Let $g_{i}$ be the indexes of $X$ belonging to group $i$, then the treatment
effect of group $i$ over group $1$ is defined as

$p_{i}=P(X_{(g_{1})_{1}}<X_{(g_{i})_{1}})+\frac{1}{2}P(X_{(g_{1})_{1}}=X_{(g_{i})_{1}}).$

The interpretation of this treatment effect is that $p_{i}$ is the probability
that a random chosen member of group $i$ has a higher value of $X$ than a
randomly chosen member of group 1, plus 0.5 times the probability that they
tie. If the two distributions are equal, then $p_{i}=0.5$. If group $i$ tends
to have higher values, then $0.5<p_{i}\leq 1$ and if it tends to have lower
values then $0\leq p_{i}<0.5$. $p_{i}$ is also known as the area under the
curve (AUC) (Mason and Graham, 2002).

Following Brunner and Munzel (2000), we define the normalized distribution
function $F_{i}(x)=\frac{1}{2}\big{(}F_{i}^{-}(x)+F_{i}^{+}(x)\big{)}$, where
$F_{i}^{-}(x)=P(X_{(g_{i})_{1}}<x)$ is the left continuous distribution
function and $F_{i}^{+}(x)=P(X_{(g_{i})_{1}}\leq x)$ is the right continuous
version. Empirical approximations $\hat{F}$ of $F$ may be estimated by
replacing the probabilities by their sample analogs. Further, we define the
mid-rank of each $X_{j}$ as $R_{j}$, and
$\bar{R}_{i}=\frac{1}{n_{i}}\sum_{j\in g_{i}}R_{j}$ to be the observed mean
rank of group $i$. An unbiased estimate of $p_{i}$ is then

$\hat{p}_{i}=\int\hat{F}_{1}d\hat{F}_{i}=\frac{1}{n_{1}}\Bigg{(}\bar{R}_{i}-\frac{n_{i}-1}{2}\Bigg{)}.$

Brunner and Munzel (2000) then make large sample inference possible by showing
that asymptotically,

$\sqrt{n}(\hat{p}_{i}-p_{i})\Rightarrow
U_{n}=\sqrt{n}\Bigg{(}\frac{1}{n_{i}}\sum_{j\in
g_{i}}F_{1}(X_{j})-\frac{1}{n_{1}}\sum_{j\in
g_{1}}F_{i}(X_{j})+1-2p_{i}\Bigg{)}.$ (1)

The right hand side of Equation 1 is the difference of two sums of independent
variables, and thus the central limit theorem may be invoked for asymptotic
normality. The variance may be expressed as

$\textrm{var}(U_{n})=n\Bigg{(}\frac{v_{i}^{2}}{n_{1}}+\frac{\sigma_{i}^{2}}{n_{i}}\Bigg{)},$

where $\sigma_{i}^{2}=\textrm{var}(F_{i}(X_{(g_{1})_{1}}))$ and
$v^{2}=\textrm{var}(F_{1}(X_{(g_{i})_{1}}))$. $\sigma$ and $v$ may be
consistently approximated by
$\hat{\sigma}_{i}^{2}=\hat{\textrm{var}}(\hat{F}_{i}(X_{(g_{1})_{1}}))$ and
$\hat{v}^{2}=\hat{\textrm{var}}(\hat{F}_{1}(X_{(g_{i})_{1}}))$, where var is
the sample variance. So, the asymptotic distribution of $\hat{p}_{i}$ may be
approximated as

$\hat{p}_{i}\sim
N\Big{(}p_{i},\frac{\hat{v}_{i}^{2}}{n_{1}}+\frac{\hat{\sigma}_{i}^{2}}{n_{i}}\Big{)}$

####  5.3.1 Extending the Test to Sequential Data

Unlike the discrete, or normal distribution cases, $\hat{p}$ are not maximum
likelihood estimators. The asymptotic distribution of $\hat{p}$ is normal, and
the diagonal terms of the covariance matrix are

$n^{-1}\hat{\Sigma}_{ii}=\frac{\hat{v}_{i}^{2}}{n_{1}}+\frac{\hat{\sigma}_{i}^{2}}{n_{i}}.$

The off-diagonal terms may be estimated noting that

$\displaystyle\textrm{cov}(\hat{p}_{l},\hat{p}_{m})$
$\displaystyle\Rightarrow\textrm{cov}\Bigg{(}\frac{1}{n_{l}}\sum_{j\in
g_{l}}F_{1}(X_{j})-\frac{1}{n_{1}}\sum_{j\in
g_{1}}F_{l}(X_{j}),\frac{1}{n_{m}}\sum_{j\in
g_{m}}F_{1}(X_{j})-\frac{1}{n_{1}}\sum_{j\in g_{1}}F_{m}(X_{j})\Bigg{)}$
$\displaystyle=\textrm{cov}\Bigg{(}\frac{1}{n_{1}}\sum_{j\in
g_{1}}F_{l}(X_{j}),\frac{1}{n_{1}}\sum_{j\in g_{1}}F_{m}(X_{j})\Bigg{)}$
$\displaystyle=\frac{1}{n_{1}}\textrm{cov}\Bigg{(}F_{l}(X_{(g_{1})_{1}}),F_{m}(X_{(g_{1})_{1}})\Bigg{)}$
$\displaystyle\approx
n^{-1}\hat{\Sigma}_{lm}=\frac{1}{n_{1}}\hat{\textrm{cov}}\Bigg{(}\hat{F}_{l}(X_{(g_{1})_{1}}),\hat{F}_{m}(X_{(g_{1})_{1}})\Bigg{)},$

where $\hat{\textrm{cov}}$ is the sample covariance. In order for the
asymptotic mSPRT

$\frac{\int\phi(\hat{p}\mid
p,n^{-1}\hat{\Sigma})g(p)dp}{\phi(\hat{p}\mid\frac{1}{2},n^{-1}\hat{\Sigma})},$
(2)

to be justified, we must show that $\hat{p}$ is asymptotically equivalent to a
generalized method of moments estimator. As with the derivation with regard to
difference in means, we let $Y_{n}$ be distributed multinomial with
probabilities $q$. Consider the function

$\gamma_{j}(X_{i},Y_{i},p_{j})=\frac{F_{1}(X_{i})\mathds{1}(Y_{i}=j)}{q_{j}}-\frac{F_{j}(X_{i})\mathds{1}(Y_{i}=1)}{q_{1}}+1-2p_{j}.$

From Brunner and Munzel (2000) we know that
$E(\gamma_{j}(X_{i},Y_{i},p_{j}))=0$, and thus the solution to

$\frac{1}{n}\sum_{i}\gamma_{j}(X_{i},Y_{i},p_{j})=0$

is a generalized method of moments estimator. Because
$\frac{n_{i}}{n}\Rightarrow q_{i}$,

$\displaystyle\sqrt{n}(\hat{p}_{j}-p_{j})$ $\displaystyle\Rightarrow U_{n}$
$\displaystyle\Rightarrow\sqrt{n}\frac{1}{n}\sum_{i}\gamma_{j}(X_{i},Y_{i},p_{j}),$

and thus the solution $\hat{p}_{j}$ is asymptotically equivalent to a
generalized method of moments estimator. Theorem 5 may then be applied to
justify Equation 2 as a sequential mSPRT.

##  6 Discussion

The field of A/B testing is just recently realizing the importance of
sequentially valid inference. Utilizing traditional statistical equations in
an environment where the analyst is checking for significance daily leads to
wildly inflated type I errors. It is thus important to understand how
sequentially valid inference may be applied to the performance indicators and
metrics of interest to the community.

The methodology developed here allows for any performance indicator that can
be seen as either a maximum likelihood estimator, or a generalized method of
moments estimator to be formulated into a sequentially valid test provided a
estimate of it’s standard error is available.

## References

  * Bartlett [1946] MS Bartlett.  The large-sample theory of sequential tests.  In _Mathematical Proceedings of the Cambridge Philosophical Society_ , volume 42, pages 239–244. Cambridge Univ Press, 1946.
  * Brunner and Munzel [2000] Edgar Brunner and Ullrich Munzel.  The nonparametric behrens-fisher problem: Asymptotic theory and a small-sample approximation.  _Biometrical Journal_ , 42(1):17–25, 2000.
  * Cox [1963] David Roxbee Cox.  Large sample sequential tests for composite hypotheses.  _Sankhyā: The Indian Journal of Statistics, Series A_ , pages 5–12, 1963.
  * Ghosh et al. [2007] Jayanta K Ghosh, Mohan Delampady, and Tapas Samanta.  _An introduction to Bayesian analysis: theory and methods_.  Springer Science & Business Media, 2007.
  * Hedges and Olkin [2014] Larry V Hedges and Ingram Olkin.  _Statistical methods for meta-analysis_.  Academic press, 2014.
  * Joanes [1972] DN Joanes.  Sequential tests of composite hypotheses.  _Biometrika_ , 59(3):633–637, 1972.
  * Johari et al. [2015] Ramesh Johari, Leo Pekelis, and David J Walsh.  Always valid inference: Bringing sequential analysis to a/b testing.  _arXiv preprint arXiv:1512.04922_ , 2015.
  * Lai [2004] Tze Leung Lai.  Likelihood ratio identities and their applications to sequential analysis.  _Sequential Analysis_ , 23(4):467–497, 2004.
  * Mason and Graham [2002] Simon J Mason and Nicholas E Graham.  Areas beneath the relative operating characteristics (roc) and relative operating levels (rol) curves: Statistical significance and interpretation.  _Quarterly Journal of the Royal Meteorological Society_ , 128(584):2145–2166, 2002.
  * Robbins et al. [1970] H Robbins, D Siegmund, LM Le Cam, J Neyman, and EL Scott.  A class of stopping rules for testing parameter hypotheses.  In _Proceedings of the Sixth Berkeley Symposium on Mathematical Statistics and Probability_ , volume 4, pages 37–41, 1970.
  * Robbins [1970] Herbert Robbins.  Statistical methods related to the law of the iterated logarithm.  _The Annals of Mathematical Statistics_ , 41(5):1397–1409, 1970.
  * Schwarz et al. [1962] Gideon Schwarz et al.  Asymptotic shapes of bayes sequential testing regions.  _The Annals of mathematical statistics_ , 33(1):224–236, 1962.
  * Wald [1945] Abraham Wald.  Sequential tests of statistical hypotheses.  _The Annals of Mathematical Statistics_ , 16(2):117–186, 1945.


# Heavy QCD axion model in light of pulsar timing arrays

Moslem Ahmadvand<EMAIL_ADDRESS>School of Particles and Accelerators,
Institute Research
in Fundamental Sciences (IPM) P. O. Box 19395-5531, Tehran, Iran  School of
Physics, Damghan University, Damghan 3671645667, Iran  Ligong Bian
<EMAIL_ADDRESS>Department of Physics and Chongqing Key Laboratory for
Strongly Coupled Physics,
Chongqing University, Chongqing 401331, China  Center for High Energy Physics,
Peking University, Beijing 100871, China  Soroush Shakeri<EMAIL_ADDRESS>Department of Physics, Isfahan University of Technology, Isfahan 84156-83111,
Iran  ICRANet, Piazza della Repubblica 10, I-65122 Pescara, Italy

###### Abstract

Recently, pulsar timing array experiments reported the observation of a
stochastic gravitational wave background in the nanohertz range frequency
band. We show that such a signal can be originated from a cosmological first-
order phase transition (PT) within a well-motivated heavy (visible) QCD axion
model. Considering the Peccei-Quinn symmetry breaking at the TeV scale in the
scenario, we find a supercooled PT, in the parameter space of the model,
prolonging the PT with the reheating temperature at the GeV scale.

##  I I. Introduction

Gravitational wave (GW) experiments have provided new paths to explore the
Universe and examine models of physics. Recently, pulsar timing array (PTA)
experiments, including the North American Nanohertz Observatory for
Gravitational Waves (NANOGrav) Agazie  _et al._ (2023), and the European
Pulsar Timing Array (EPTA) Antoniadis  _et al._ (2023a) along with the Parkes
Pulsar Timing Array (PPTA) Reardon  _et al._ (2023) and Chinese Pulsar Timing
Array Xu  _et al._ (2023) released their latest results observing significant
evidence of a GW following the Hellings-Downs pattern Hellings and Downs
(1983) in the angular cross-correlation of pulsar timing residuals, supporting
a stochastic gravitational background (SGWB).

Various sources can be the origin of these signals, for instance, supermassive
black hole binary mergers, have been proposed as astrophysical candidates
Afzal  _et al._ (2023); Antoniadis  _et al._ (2023b); Rajagopal and Romani
(1995); Ellis  _et al._ (2023a); Bi  _et al._ (2023). However, such stochastic
signals may come from a relic cosmic background Madge  _et al._ (2023) similar
to the cosmic microwave background (CMB). In addition, cosmological sources
such as first-order phase transition (FOPT) around the QCD scale Ellis  _et
al._ (2023b); Ratzinger and Schwaller (2021); Addazi  _et al._ (2023); Ghosh
_et al._ (2023); Rezapour  _et al._ (2022a, b); Neronov  _et al._ (2021);
Ahmadvand and Bitaghsir Fadafan (2018, 2017); Athron  _et al._ (2023); Han
_et al._ (2023); Li and Xie (2023), cosmic strings Blasi  _et al._ (2021);
Ellis and Lewicki (2021); Buchmuller  _et al._ (2020); Lazarides  _et al._
(2023); Wang  _et al._ (2023); Bian  _et al._ (2022a), domain walls Bai  _et
al._ (2023); Kitajima  _et al._ (2023); Blasi  _et al._ (2023); Li  _et al._
(2023); Bian  _et al._ (2022b); Babichev  _et al._ (2023); Sakharov  _et al._
(2021); Zhang  _et al._ (2023), primordial black holes Depta  _et al._ (2023);
Gouttenoire  _et al._ (2023); Franciolini  _et al._ (2023); Guo  _et al._
(2023); Hosseini Mansoori  _et al._ (2023) and inflation Firouzjahi and Riotto
(2023); Vagnozzi (2021); Ashoorioon  _et al._ (2022); Niu and Rahat (2023);
Servant and Simakachorn (2023) may fit the recent data better Agazie  _et al._
(2023). There have also been some recent discussions about the lake of strong
preferences for any specific SGWB sources with the current data based on
Bayesian analysis Bian  _et al._ (2023). However, future PTA datasets have the
potential to distinguish between the cosmological and astrophysical scenarios,
for example, by anisotropies, circular polarization or individual sources
Ellis  _et al._ (2023c).

An intriguing candidate that we propose in this paper is a supercooling phase
transition (PT) which ends around the $\mathrm{GeV}$ scale; such a PT can
naturally be accommodated in the context of heavy QCD axion models Rubakov
(1997); Di Luzio  _et al._ (2020). A cosmic FOPT may be accomplished due to a
spontaneous broken symmetry and the process is accompanied by the nucleation
of bubbles separating true and false vacua. In a supercooled PT, the vacuum
energy is dominated and the Universe remains in the false vacuum for a long
period and supercools, increasing the PT duration Ellis  _et al._ (2019);
Brdar  _et al._ (2019); Kobakhidze  _et al._ (2017); Von Harling  _et al._
(2020); Athron  _et al._ (2023); Yang  _et al._ (2023).

In this work, we consider a scenario within heavy axion models which are very
appealing in that not only can the strong CP problem Peccei and Quinn (1977)
be addressed but also the models have richer phenomenology Berezhiani  _et
al._ (2001); Dimopoulos  _et al._ (2016) relative to the invisible QCD axion
models Zhitnitsky (1980); Dine  _et al._ (1981); Kim (1979); Shifman  _et al._
(1980). Moreover, the relation of the axion mass, $m_{a}$, and its decay
constant, $f_{a}$ is modified so that larger $m_{a}$ and lower $f_{a}$ are
allowed. These models are also motivated in connection with the Peccei-Quinn
(PQ) quality problem Georgi  _et al._ (1981); Kamionkowski and March-Russell
(1992); Holman  _et al._ (1992) (see Appendix A).

We explore a supercooled PT with a U(1) PQ symmetry breaking at the TeV scale.
Based on the nearly conformal dynamics of the PQ scalar field, we analytically
find important parameters encoding the GW spectrum and show that the
corresponding signals can be well fitted to the PTA data.

##  II II. PQ phase transition in the model

The CP-violating $\bar{\theta}$ parameter in the strong interactions is
experimentally bounded $\bar{\theta}\lesssim 10^{-10}$ Baker  _et al._ (2006)
and the nature of this smallness, known as the strong CP problem, is a big
puzzle in particle physics. One of the interesting solutions to this problem
is based on a $U(1)$ global symmetry, $U(1)_{\mathrm{PQ}}$, first proposed by
Peccei and Quinn Peccei and Quinn (1977). At energies higher than the
electroweak (EW) scale, the symmetry is spontaneously broken and the axion as
the pseudo-Nambu-Goldstone boson is generated. The interaction of the axion
with gluons at the QCD scale, $(a/f_{a}+\bar{\theta})G\widetilde{G}$, and as a
result the vacuum expectation value of the axion can cancel the $\bar{\theta}$
term Peccei (2008). Considering the astrophysical constraints on the axion
decay constant, $10^{8}\,\mathrm{GeV}\lesssim f_{a}\lesssim
10^{17}\,\mathrm{GeV}$ Raffelt (2008); Arvanitaki  _et al._ (2010); Shakeri
and Hajkarim (2023), there are two classes of models, known as invisible axion
models Zhitnitsky (1980); Dine  _et al._ (1981); Kim (1979); Shifman  _et al._
(1980). However, as mentioned before, $f_{a}$ can be lowered,
$f_{a}\sim(1-100)\,\mathrm{TeV}$ in the so-called visible axion models, still
addressing the strong CP problem. In this sense, there are diverse scenarios,
e.g. enlarging the QCD color gauge group $SU(3+\mathcal{N})$ Gherghetta  _et
al._ (2016) or assuming some hidden sector with a confining scale larger than
the QCD cutoff $\Lambda_{\mathrm{QCD}}$ such that additional terms are
supplemented to the axion potential and change the axion mass Berezhiani  _et
al._ (2001); Fukuda  _et al._ (2015). In the case of adding a hidden sector,
one can consider a whole copy of the standard model (SM), with the $SM\times
SM^{\prime}$ gauge group, and particles of the two sectors can interact
gravitationally or by some very feebly couplings. Furthermore, a $Z_{2}$
mirror symmetry between particles of the two sectors can be imposed and soft
terms breaking the symmetry can induce
$\Lambda^{\prime}_{\mathrm{QCD}}\gg\Lambda_{\mathrm{QCD}}$ Berezhiani  _et
al._ (2001); Fukuda  _et al._ (2015).

Here, we consider a DFSZ axion model Zhitnitsky (1980); Dine  _et al._ (1981),
containing a gauge-singlet PQ scalar field, $\Phi$, and two Higgs doublets
under $SU(2)_{L}$, $H_{u}$ and $H_{d}$, which is supplemented with its hidden
copy. Indeed, we consider a theory with the $SM\times SM^{\prime}$ gauge group
and each sector interacting with the PQ scalar (see Appendix B) and has the
same PQ symmetry.111An analogous procedure can also be applied to a KSVZ axion
model Kim (1979); Shifman  _et al._ (1980). Therefore, the axion couples to
both QCD sectors

$\frac{a}{f_{a}}\left(G\widetilde{G}+G^{\prime}\widetilde{G}^{\prime}\right),$
(1)

and its mass would be $m_{a}\sim\Lambda_{\mathrm{QCD}}^{\prime 2}/f_{a}$
Berezhiani  _et al._ (2001); Fukuda  _et al._ (2015) with
$f_{a}\sim\mathrm{TeV}$. In this work, we do not go through axion interactions
and focus on the UV theory and study the PQ PT associated with the PQ symmetry
breaking. We explore a supercooled PQ PT along the direction of PQ scalar
dynamic, $\langle\Phi\rangle=v_{\phi}/\sqrt{2}$. We consider the case where
mass parameters are small, $\mu_{d}^{2},\mu_{u}^{2},\lambda_{\phi}f_{a}^{2}\ll
f_{a}^{2}$ (see Appendix B) with the similar condition in the hidden sector,
and main contributions to the potential arises from the one-loop Coleman-
Weinberg (CW) quantum correction Coleman and Weinberg (1973) and the leading
contribution from thermal corrections Quiros (1999); thereby masses and the
symmetry breaking scale are generated radiatively Gildener and Weinberg
(1976).

At the zero-temperature limit, quantum corrections (see Appendix C) contribute
to the potential that is approximately scale invariant as

$V=\left(A+B\ln\frac{\phi^{2}}{f_{a}^{2}}\right)\phi^{4}\equiv\lambda(\phi)\phi^{4}.$
(2)

where

$A=\frac{1}{64\pi^{2}}\left(\kappa_{u}^{2}\left(\ln(\kappa_{u})-\frac{3}{2}\right)+\kappa_{u}^{\prime
2}\left(\ln(\kappa^{\prime}_{u})-\frac{3}{2}\right)\right),$ (3)
$B=\frac{1}{64\pi^{2}}\left(\kappa_{u}^{2}+\kappa_{u}^{\prime 2}\right).$ (4)

For the sake of simplicity, we only considered $\kappa_{u}$ and
$\kappa_{u}^{\prime}$ couplings which are determined at the vacuum and are
then evolved to a scale $\phi$ through the beta function (see Appendix C).
Thus, at the nonzero vacuum, $v_{\phi}(T=0)\equiv f_{a}$, from
$(dV/d\phi)|_{f_{a}}=0$ we find

$-4\lambda=\frac{d\lambda}{d\ln\phi},~{}~{}~{}~{}~{}\Delta V\equiv
V(0)-V(f_{a})=\frac{Bf_{a}^{4}}{2}=\frac{(\bar{\kappa}_{u}^{2}+\bar{\kappa}_{u}^{\prime
2})f_{a}^{4}}{128\pi^{2}}.$ (5)

Considering thermal corrections (see Appendix C) and high temperatures, the
origin point would be the minimal, and the potential would be

$V=DT^{2}\phi^{2}+\lambda\phi^{4}+\cdots,$ (6)

where

$D=\frac{1}{24}\left(\kappa_{u}+\kappa_{u}^{\prime}\right).$ (7)

As temperature goes down, negative values of $\lambda$ can change the
potential slope and induce a bump so that at the critical temperature,
$T_{c}$, two degenerate states, $\Delta V(T)=V(0,T)-V(v_{\phi}(T),T)=0$,
indicating a FOPT, can be generated due to these values of $\lambda$, Fig.
(1). At temperatures below $T_{c}$, the vacuum $v_{\phi}(T)$ is the favorable
one and at some temperature $T_{n}$ bubbles of the new phase are nucleated. In
fact, bubble nucleation occurs when the bubble formation probability per unit
Hubble space-time volume, which is proportional to
$\exp\left(-S_{3}(T)/T\right)$ where $S_{3}(T)$ is the bounce action
quantifying the tunneling process, would be of the order of 1. As a result, we
can find the nucleation temperature $T_{n}$ via the relation Linde (1983);
Arnold and Vokos (1991)

$\frac{S_{3}(T_{n})}{T_{n}}\sim 4\ln\Big{(}\frac{T_{n}}{H(T_{n})}\Big{)},$ (8)

where the Hubble parameter $H$ is given by

$H^{2}=\frac{1}{3M_{\mathrm{Pl}}^{2}}\left(\xi T^{4}+\Delta V(T)\right),$ (9)

where the reduced Planck mass is $M_{\mathrm{pl}}\simeq 2.43\times
10^{18}\,\mathrm{GeV}$, $\xi=\pi^{2}g_{*}/30$, and $g_{*}\simeq 107$ is the
effective number of relativistic degrees of freedom, assuming that it does not
change during the PT. At low temperatures, $T\ll T_{c}$, it is expected
$v_{\phi}(T)=f_{a}$, and $\Delta V(T)$ is well approximated by $\Delta V$.
Furthermore, in this limit, the bounce action can be well approximated by
Witten (1981)

$\frac{S_{3}}{T}\simeq-18.9\frac{\sqrt{2D}}{\lambda(T)}\simeq-5.4\frac{\sqrt{\kappa_{u}+\kappa_{u}^{\prime}}}{\lambda(T)}$
(10)

and hence we can find the nucleation temperature as

$4\ln\Big{(}\frac{\sqrt{3}M_{\mathrm{Pl}}T_{n}}{\sqrt{\xi T_{n}^{4}+\Delta
V}}\Big{)}\simeq-5.4\frac{\sqrt{\kappa_{u}+\kappa_{u}^{\prime}}}{\lambda(T)}.$
(11)

At the time when the transition completes, from the conservation of energy,
one obtains $\rho_{R}(T_{*})=\rho_{R}(T_{p})+\Delta V$ Ellis  _et al._ (2019)
where $\rho_{R}(T)=\xi T^{4}$ is the radiation energy and $T_{p}$ is the
percolation temperature. Assuming $T_{p}\simeq T_{n}$, the reheating
temperature can be computed by

$T_{*}^{4}=\frac{\Delta V}{\xi}+T_{n}^{4}.$ (12)

For the case of fast reheating $H_{*}\equiv H(T_{*})\simeq H(T_{n})$. Another
important quantity in characterizing the generated GWs is the inverse of PT
duration calculated by the following relation

$\frac{\beta}{H_{*}}=T_{n}\frac{d}{dT}\left(\frac{S_{3}(T)}{T}\right)\Bigg{|}_{T_{n}}\simeq-\frac{S_{3}(T_{n})}{T_{n}}\frac{\beta_{\lambda}(T_{n})}{\lambda(T_{n})},$
(13)

where $(d\lambda(T)/d\ln T)|_{T_{n}}=\beta_{\lambda}(T_{n})$. From Eq. (6) and
$V(\phi_{0},T)=0$, for $T\ll T_{c}$, and $\phi_{0}\sim T$ we obtain
$\lambda(T_{n})\sim A+B\ln(T_{n}/f_{a})$. Therefore, expressing
$-\lambda(T_{n})/\beta_{\lambda}(T_{n})\sim\gamma\ln(f_{a}/T_{n})$, we find

$\frac{\beta}{H_{*}}\simeq\frac{4}{\gamma\ln(f_{a}/T_{n})}\ln\Big{(}\frac{\sqrt{3}M_{\mathrm{Pl}}T_{n}}{\sqrt{\xi
T_{n}^{4}+\Delta V}}\Big{)}.$ (14)

Expecting $\beta/H_{*}\sim\mathcal{O}(1-10)$ for the case of supercooling, we
fix the parameter $\gamma$ with the data.

The strength of the supercooled PT is also obtained by

$\alpha=\frac{\Delta V}{\rho_{R}(T_{n})}.$ (15)

In the next section, calculating the aforementioned quantities, we obtain the
GW energy density spectrum.

Figure 1: For three different values of temperature as labeled and
representative values of $A=-0.0053$, $B=0038$, and $D=0.09$, the radiatively
induced potential is illustrated. Above $T_{c}$, the origin is the only
minimum. At $T_{c}$, two degenerate states appear and at lower temperatures
the nonzero minimum is favored. The Universe remains in the false vacuum until
the tunneling process is well developed at low temperatures due to
$\lambda(T)$.

##  III III. SGWB signal

In this section, based on the obtained PT quantities, we calculate the GW
spectrum of the supercooled PT. In this sense, sources contributing to the GW
production are bubble wall collisions and fluid motions. Relying on the recent
numerical computation Lewicki and Vaskonen (2023), the present-day GW energy
density spectrum is obtained by

$h^{2}\Omega_{\mathrm{GW}}(f)=1.67\times
10^{-5}\Big{(}\frac{H_{*}}{\beta}\Big{)}^{2}\Big{(}\frac{\kappa\alpha}{1+\alpha}\Big{)}^{2}\Big{(}\frac{100}{g_{*}}\Big{)}^{\frac{1}{3}}S(f),$
(16)

where the spectral shape of the GW is obtained by

$S(f)=\frac{\bar{A}(a+b)^{c}}{\left[b\left(\frac{f}{f_{p,0}}\right)^{-\frac{a}{c}}+a\left(\frac{f}{f_{p,0}}\right)^{\frac{b}{c}}\right]^{c}},$
(17)

and $\kappa$ is the fraction of the vacuum energy converted to the kinetic
energy of bubble walls, $\kappa_{\mathrm{vac}}\simeq 1/(1+5/(\beta
R_{\mathrm{eq}}))$, and fluid motions
$\kappa_{\mathrm{f}}=1-\kappa_{\mathrm{vac}}$ Lewicki and Vaskonen (2023). As
for the spectral shape, Eq. (17), the spectrum’s low-and high-frequency power-
law tails are controlled by $a$ and $b$, its transition width from one power
law to another is determined by $c$, and its amplitude depends on $\bar{A}$.

In our case with $\alpha\gg 1$, bubbles collide in the vacuum and the bubble
wall velocity can reach the speed of light; thus we consider
$\kappa_{\mathrm{vac}}\simeq 1$ corresponding to $\beta R_{\mathrm{eq}}\gg 1$.
Therefore, according to Table I of Ref. Lewicki and Vaskonen (2023),
corresponding fit parameters of the spectral function appropriate for a global
broken symmetry, $T_{rr}\propto R^{-2}$ (where $T_{rr}$ is the maximum of the
stress-energy tensor and $R$ is the bubble radius), would be as $\bar{A}\simeq
5.93\times 10^{-2}$, $a\simeq 1.03$, $b\simeq 1.84$, and , $c\simeq 1.45$,
whereas for the envelope approximation, $\bar{A}\simeq 3.78\times 10^{-2}$,
$a\simeq 3.08$, $b\simeq 0.98$, and , $c\simeq 1.91$.

The present redshifted peak frequency is given by

$f_{p,0}=16.5\times
10^{-6}[\mathrm{Hz}]\Big{(}\frac{f_{p}}{\beta}\Big{)}\Big{(}\frac{\beta}{H_{*}}\Big{)}\Big{(}\frac{T_{*}}{100~{}\mathrm{GeV}}\Big{)}\Big{(}\frac{g_{*}}{100}\Big{)}^{\frac{1}{6}},$
(18)

where for the case of $T_{rr}\propto R^{-2}$, $f_{p}/\beta\simeq 0.64/(2\pi)$
and for the envelope approximation $f_{p}/\beta\simeq 1.33/(2\pi)$ Lewicki and
Vaskonen (2023).

As can be seen from Fig. 2, using datasets of NANOGrav 15 yr, EPTA release 2,
and PPTA data release 3 (DR3), for $f_{a}=1\,\mathrm{TeV}$, and some
representative values of parameters,
$\bar{\kappa}_{u}=\bar{\kappa}^{\prime}_{u}=0.0001$,
$\kappa_{u}=\kappa^{\prime}_{u}=0.001$ and $\lambda(T)=-0.0015$, which
correspond to $T_{n}=100\,\mathrm{MeV}$, $T_{*}=1\,\mathrm{GeV}$,
$\alpha=8995$, and $\beta/H_{*}=17,3.4,1.7$ with the corresponding parameters
$\gamma=1,5,10$, respectively, the GW signals can be very consistent with the
PTA data from the bubble collision spectra for both the spectral shape of
$T_{rr}\propto R^{-2}$ (solid line) and for the envelope (dashed line).
Remarkably, we see that the high-frequency range of the GW signals with the
envelope approximation falls within the sensitivity range of future space-
based GW experiments such as Laser Interferometer Space Antenna (LISA), Deci-
hertz Interferometer Gravitational Wave Observatory (DECIGO), and Big Bang
Observer (BBO) and can be probed by these detectors.

Furthermore, using the PTA data, we perform a best-fit analysis over a range
of $\gamma$ parameter ($\beta/H_{*}$) through the $\chi^{2}$ test based on the
relation

$\chi^{2}=\sum_{i=1}^{N}\frac{\left(\log_{10}\Omega_{\mathrm{th}}h^{2}-\log_{10}\Omega_{\mathrm{exp}}h^{2}\right)^{2}}{2\bar{S_{i}}^{2}},$
(19)

where $\Omega_{\mathrm{th}}h^{2}$ denotes the GW predicted by the model,
$\Omega_{\mathrm{exp}}h^{2}$ represents the observed GW signal by PTA
experiments, and $\bar{S_{i}}$ is the deviation from the midpoint value of
each data point in $\log_{10}\Omega_{\mathrm{exp}}h^{2}$ within the
uncertainty range. As shown in Figs. 3, 4 the best-fit point of the GW energy
spectrum taking into account all datasets derived from the $T_{rr}\propto
R^{-2}$ case for the aggregated dataset is $\gamma=3.28$ ($\beta/H_{*}=5.1$)
and for the envelope case, the best-fit point would be $\gamma=13.3$
($\beta/H_{*}=1.3$).

Figure 2:  GW spectra with the spectral shape for the case $T_{rr}\propto
R^{-2}$ (solid line) and for envelope approximation (dashed line) are
displayed. Three different values $\gamma=1,5,10$ correspond to the inverse
duration of the PT $\beta/H_{*}=17,3.4,1.7$, respectively. The recent results
of PTA experiments are also shown. The shaded regions represent the
sensitivity range of ongoing GW experiments: LISA Amaro-Seoane (2017), DECIGO
et al. (2011), BBO Crowder and Cornish (2005), Einstein Telescope (ET) Punturo
_et al._ (2010) and Cosmic Explorer (CE) Reitze  _et al._ (2019). Figure 3:
Based on the GWs with the spectral shape of $T_{rr}\propto R^{-2}$, we the
find best-fit point in terms of $\gamma$ parameter, which is 2.52 for EPTA
(solid purple), 2.17 for NANOGrav (dashed green), and 1.81 for PPTA (dotted
orange). The best-fit value for the aggregated dataset is given by 3.28 (solid
blue) corresponding to $\beta/H_{*}=5.1$. The insert shows GW spectra for the
best-fit values of $\gamma$ in addition to observed signals.  Figure 4: Based
on the GWs with the envelope approximation.we find the best-fit point in terms
of $\gamma$ parameter, which is 13.55 for EPTA (solid purple), 9.49 for
NANOGrav (dashed green), and 7.93 for PPTA (dotted orange). The best-fit value
for the aggregated dataset is given by 13.30 (solid blue) corresponding to
$\beta/H_{*}=1.3$. The insert shows GW spectra for the best-fit values of
$\gamma$ in addition to observed signals.

##  IV IV. Discussion

Another PQ PT consequence of the model is that the symmetry breaking may
induce a similar or different EW symmetry breaking scale in the ordinary and
hidden sectors. For instance, as for $v_{\mathrm{EW}}^{\prime}$, this can be
realized by these terms of the Lagrangian
$\kappa^{\prime}_{u}h_{u}^{2}f_{a}^{2}+\lambda^{\prime}_{u}h_{u}^{4}$,
expressed in terms of the electrically neutral component of the fields, and
hence $v_{\mathrm{EW}}^{\prime
2}\simeq\kappa^{\prime}_{u}f_{a}^{2}/\lambda^{\prime}_{u}$. Setting
$\kappa^{\prime}_{u}=0.0001$ and $\lambda^{\prime}_{u}=0.1$, we obtain
$v_{\mathrm{EW}}^{\prime}=1\,\mathrm{GeV}$. Even for these values of
$v_{\mathrm{EW}}^{\prime}$ its QCD scale could be
$\Lambda^{\prime}_{\mathrm{QCD}}\gtrsim v^{\prime}_{\mathrm{EW}}$ Fukuda  _et
al._ (2015). Consequently, in this case, we obtain the axion mass
$m_{a}\gtrsim 1\,\mathrm{MeV}$, while considering $v_{\mathrm{EW}}\lesssim
v^{\prime}_{\mathrm{EW}}<f_{a}$ and assuming
$\Lambda^{\prime}_{\mathrm{QCD}}\lesssim v_{\mathrm{EW}}$, one finds
$m_{a}\lesssim 10\,\mathrm{GeV}$.

Heavy axions may interact with visible and hidden (dark) photons

$\frac{a}{f_{a}}\left(F\widetilde{F}+F^{\prime}\widetilde{F}^{\prime}\right),$
(20)

so that the axion decay rate would be

$\Gamma_{a\gamma\gamma}=\frac{g_{a\gamma}^{2}m_{a}^{3}}{64\pi},\quad\Gamma_{a\gamma^{\prime}\gamma^{\prime}}=\frac{g_{a\gamma}^{\prime
2}m_{a}^{3}}{64\pi}$ (21)

where

$g_{a\gamma}\simeq\frac{\bar{\alpha}}{2\pi}\frac{E}{Nf_{a}}\frac{2z}{1+z},~{}~{}~{}~{}~{}~{}g_{a\gamma}^{\prime}\simeq\frac{\bar{\alpha^{\prime}}}{2\pi}\frac{E}{Nf_{a}}\frac{2z^{\prime}}{1+z^{\prime}},$
(22)

$\bar{\alpha}=1/137$ is the fine structure constant; $z=m_{u}/m_{d}$; and $E$
and $N$ are the electromagnetic and QCD anomaly coefficients, respectively, Di
Luzio  _et al._ (2020). The heavy QCD axions can also couple to SM fermions

$\frac{c_{f}}{2f_{a}}\partial_{\mu}a\bar{f}\gamma^{\mu}\gamma_{5}f+\frac{c^{\prime}_{f}}{2f_{a}}\partial_{\mu}a\bar{f^{\prime}}\gamma^{\mu}\gamma_{5}f^{\prime}.$
(23)

Decaying the heavy axions to these photons contributes to the effective number
of relativistic species which is strongly constrained by CMB measurements,
$N_{\mathrm{eff}}=2.96^{+0.34}_{-0.33}$ Aghanim  _et al._ (2020), putting the
most stringent constraints on the heavy QCD axion mass. According to the
analysis of Ref. Dunsky  _et al._ (2022), for $E/N=8/3$ which may be the case
of a DFSZ model Ahmadvand (2021), $m_{a}<100\,\mathrm{MeV}$ would be excluded.
However, lower masses, $1\,\mathrm{MeV}<m_{a}<100\,\mathrm{MeV}$ where the
dark photon contribution to $N_{\mathrm{eff}}$ compensates neutrino dilution
from axion decays to photons, are allowed for $E/N=1/3$.

We assume a baryon asymmetric dark sector. Thus, considering
$v_{\mathrm{EW}}\lesssim v^{\prime}_{\mathrm{EW}}<f_{a}$, the mass range and
relic abundance of stable particles in the hidden sector are approximately the
same as those of the ordinary sector. As a result, the relic abundance of
these particles cannot give rise to cosmological problems. Moreover, we assume
that in the hidden sector there is no seesaw mechanism, which allows the Dirac
neutrino mass in this sector as Fukuda  _et al._ (2015)

$m_{v^{\prime}}\sim\frac{v^{\prime}_{\mathrm{EW}}}{v_{\mathrm{EW}}}\sqrt{M_{R}m_{v}},$
(24)

where $M_{R}$ is the right-handed neutrino mass in the ordinary sector. For
$M_{R}\gtrsim 10^{7}\,\mathrm{GeV}$ and $m_{a}\gtrsim 1\,\mathrm{MeV}$
($\Lambda^{\prime}_{\mathrm{QCD}}\gtrsim 1\,\mathrm{GeV}$), the hidden
neutrinos decay to the hidden electron and pion, $v^{\prime}\rightarrow
e^{\prime}+\pi^{\prime}$, and thereby these neutrinos do not contribute to
$N_{eff}$.

Concerning astrophysical bounds, the model parameters can be constrained from
SN 1987A Raffelt (1990), with the supernova core cooling rate due to axion
emission. According to Eq. (21), one can obtain the axion mean free path as

$\lambda_{mfp}=\frac{vE_{a}}{m_{a}\Gamma_{a\gamma\gamma}}\approx 10\
\text{m}\left(\frac{f_{a}}{\text{TeV}}\right)^{2}\left(\frac{E_{a}}{200\
\text{MeV}}\right)\left(\frac{100\text{MeV}}{m_{a}}\right)^{4}$ (25)

where $E_{a}$ stands for typical energies of the axions. Hence, for
$E_{a}\gtrsim 200\,\mathrm{MeV}$, $f_{a}=1\,\mathrm{TeV}$ and $m_{a}\gtrsim
100\,\mathrm{MeV}$, we find that $\lambda_{mfp}$ is much smaller than the
approximate radius of the supernova core $R\sim 10\,\mathrm{km}$. Thus,
because of the scattering and reabsorption, the axions are trapped within the
supernova, and cannot take away any energy Essig  _et al._ (2010); Jaeckel
_et al._ (2018).

In addition, the emission of dark photons can lead to a very fast core cooling
rate of the supernova and the associated energy loss rate, which is bounded
$L_{\gamma^{\prime}}<10^{52}\,\mathrm{erg}/\mathrm{s}$, and can be obtained as
$L_{\gamma^{\prime}}\simeq\Gamma_{a\gamma^{\prime}\gamma^{\prime}}m_{a}(4\pi
R^{3}/3)(1.2T^{3}/\pi^{2})$ Berezhiani  _et al._ (2001). Thus, taking into
account a core temperature $T\simeq 10\,\mathrm{MeV}$,
$f_{a}=1\,\mathrm{TeV}$, $m_{a}=1\,\mathrm{MeV}$, and $z^{\prime}=z\simeq
0.5$, the bound can be satisfied, while for $m_{a}\gtrsim 100\,\mathrm{MeV}$,
it may require $z^{\prime}\lesssim 10^{-4}$.

Considering the axion-electron interactions, the astrophysical bounds on
helium burning lifetime of the horizontal branch stars imply $g_{ae}\sim
m_{e}/f_{a}<(2.5-6)\times 10^{-13}$ Berezhiani  _et al._ (2001); Di Luzio  _et
al._ (2020), and hence restrict $f_{a}$ to its invisible values. However, this
bound does not apply to the heavy axions with $m_{a}>300\,\mathrm{keV}$. In
fact, in the stellar cores whose typical temperatures reach
$10\,\mathrm{keV}$, the axion production rate is suppressed by the exponential
factor $\exp(-m_{a}/T)$. Therefore, astrophysical bounds allow heavy axions
with masses above $100\,\mathrm{MeV}$ for $f_{a}=1\,\mathrm{TeV}$.

The viable parameter space of the model can also be probed via other
experiments such as proton beam dump Bergsma  _et al._ (1985); Aloni  _et al._
(2019), kaon decays Alves and Weiner (2018); Ahn  _et al._ (2019), and
colliders Knapen  _et al._ (2017). For $m_{a}\gtrsim 3m_{\pi}$, the heavy
axions may be produced via $B\rightarrow K\,a$, where $B$ and $K$ denote $B$
mesons and kaons, respectively, and then mainly decay into three pions. For
$m_{a}\lesssim 3m_{\pi}$, the dominant decay channels are the axion decay into
two photons, electrons, and muons, for $m_{a}>2m_{e}$ and $m_{a}>2m_{\mu}$,
respectively. In this sense, $m_{a}\lesssim 400\,\mathrm{MeV}$ for $f_{a}$ at
the TeV scale is excluded by experiments such as kaon decays and the proton
beam dump, while for $m_{a}\gtrsim 400\,\mathrm{MeV}$ the search is
challenging due to the dominance of the hadronic decay mode Bertholet  _et
al._ (2022).

It is noteworthy that the heavy axion cannot be a dark matter candidate in
this context. However, some stable particles in the hidden sector may
contribute to the dark matter relic abundance.222Moreover, according to Eq.
(VIII), the interaction of the Higgs boson with the heavy axions can
contribute to its invisible decays Aad  _et al._ (2022); Tumasyan  _et al._
(2022), and since the heavy axions are unstable, they can decay to stable
particles in the hidden sector. We leave detailed calculations in connection
with axion interactions and possible dark matter candidates in the presented
scenario to a future work.

##  V V. Conclusion

With regard to the recently detected GW signals by PTA experiments, we have
proposed a heavy QCD axion scenario, with a PQ symmetry breaking scale around
the TeV scale, based on a (DFSZ) axion model which is supplemented with its
hidden copy. We investigated a supercooled PQ FOPT which is derived by CW
quantum corrections and analytically found important PT quantities, including
the reheating temperature, the inverse duration of the PT, and the strength of
the PT. We have shown that within the parameter space of the model the
generated GWs from the bubble wall collisions can be consistently fitted with
the recently observed PTA data. Furthermore, it is shown that the high-
frequency range of such GW signals can be probed by future space-based GW
detectors.

In addition, Heavy axion models are well-motivated scenarios allowing the
$(m_{a},f_{a})$ relation to be relaxed, still addressing the strong CP
problem. It is shown that considering the EW symmetry breaking scale in the
ordinary and hidden sectors induced after the PQ symmetry breaking, a range of
axion masses, $1\,\mathrm{MeV}\lesssim m_{a}\lesssim 10\,\mathrm{GeV}$ can be
obtained. Nevertheless, viable regions of the mass space are limited via CMB
observation due to the dark photon contribution to $N_{\mathrm{eff}}$, and
also other experimental setups such as the proton beam dump.

As a result, interpreting the recent GW data, the model provides a setup which
can be probed by future GW detectors, CMB telescopes and collider experiments.

##  VI ACKNOWLEDGMENTS

SS thanks Fazlollah Hajkarim and Seyed Mohammad Mahdi Sanagostar for helpful
discussions about data analysis. This work is supported in part by the
National Key Research and Development Program of China under Grant No.
2021YFC2203004, and the National Natural Science Foundation of China (NSFC)
under Grants No. 12075041 and No. 12147102.

##  VII Appendix A: The PQ quality problem

Invisible axion models with $f_{a}\gtrsim 10^{8}\,\mathrm{GeV}$ suffer from
the axion quality problem. Indeed, it is expected that the PQ symmetry is
explicitly broken by higher-dimensional Planck-suppressed operators which can
spoil the axion solution. In our model, the PQ symmetry may be explicitly
broken by the following five-dimensional operator

$\frac{\lambda_{a}}{5!M_{P}}\phi_{a}^{5}$ (26)

where $M_{P}$ denotes the Planck mass and $\phi_{a}$ is the radial component
of the PQ scalar. Such an explicit PQ symmetry breaking term induces a mass
term in the axion field, $m_{*}^{2}\sim 10^{-2}\lambda_{a}(f_{a}^{3}/M_{P})$,
and moves the axion potential minimum away from the CP conserving minimum,
$\Delta\bar{\theta}\sim m_{*}^{2}/m_{a}^{2}$ Di Luzio  _et al._ (2020).
Therefore, assuming $\lambda_{a}=\mathcal{O}(1)$, for $f_{a}=1\,\mathrm{TeV}$
and $m_{a}\gtrsim 100\,\mathrm{MeV}$, the model satisfies the current bound
$\Delta\bar{\theta}\lesssim 10^{-10}$.

##  VIII Appendix B: The scalar potential

The tree-level scalar potential with the PQ symmetry is given by Ahmadvand
(2021)

$\displaystyle
V(\Phi,H_{u},H_{d})=\lambda_{\phi}\left(|\Phi|^{2}-f^{2}\right)^{2}+\left|H_{d}\right|^{2}\left(\kappa_{d}\,|\Phi|^{2}-\mu_{d}^{2}\right)$
$\displaystyle+\left|H_{u}\right|^{2}\left(\kappa_{u}\,|\Phi|^{2}-\mu_{u}^{2}\right)-\left(\kappa_{ud}\,\Phi
H_{u}^{\dagger}\,H_{d}+\text{H.c.}\right)$
$\displaystyle+\lambda_{d}\left|H_{d}\right|^{4}+\lambda_{u}\left|H_{u}\right|^{4}+\lambda_{ud}\left(\left|H_{u}\right|^{2}\left|H_{d}\right|^{2}-\left|H_{u}^{\dagger}H_{d}\right|^{2}\right).$
(27)

We also assume an analogous potential with $H^{\prime}_{d}$ and
$H^{\prime}_{u}$ from the hidden sector.

##  IX Appendix C: Quantum and thermal corrections

The one-loop Coleman-Weinberg quantum correction is given by

$\displaystyle V_{\mathrm{CW}}\left(\phi\right)$
$\displaystyle=\sum_{i}(-1)^{F_{b/f}}g_{i}\frac{m_{i}^{4}\left(\phi\right)}{64\pi^{2}}\left[\ln\left(\frac{m_{i}^{2}\left(\phi\right)}{\Lambda^{2}}\right)-c_{i}\right]$
(28)

where $F_{b/f}=1(0)$ for fermions (bosons), $g_{i}$ is the number of degrees
of freedom for a given field, and $c_{i}=3/2(5/2)$ for scalars and fermions
(vectors).

For our model, considering only $\kappa_{u}$ and $\kappa^{\prime}_{u}$
couplings, $m^{2}(\phi)=\kappa_{u}\phi^{2}$ (and the same with
$\kappa^{\prime}_{u}$), one can write

$V_{\mathrm{CW}}=\left(A+B\ln\frac{\phi^{2}}{\Lambda^{2}}\right)\phi^{4}\equiv\lambda(\phi)\phi^{4},$
(29)

where

$A=\frac{1}{64\pi^{2}}\left(\kappa_{u}^{2}\left(\ln(\kappa_{u})-\frac{3}{2}\right)+\kappa_{u}^{\prime
2}\left(\ln(\kappa^{\prime}_{u})-\frac{3}{2}\right)\right),$ (30)
$B=\frac{1}{64\pi^{2}}\left(\kappa_{u}^{2}+\kappa_{u}^{\prime 2}\right).$ (31)

The running couplings at the scale $\phi$ are also obtained by the leading-
order beta function

$\frac{d\kappa_{u}}{d\ln\phi}\simeq\frac{\kappa_{u}^{2}}{16\pi^{2}}+\frac{\kappa_{u}\kappa_{u}^{\prime}}{8\pi^{2}}$
(32) $\frac{d\kappa_{u}^{\prime}}{d\ln\phi}\simeq\frac{\kappa_{u}^{\prime
2}}{16\pi^{2}}+\frac{\kappa_{u}\kappa_{u}^{\prime}}{8\pi^{2}}$ (33)

The one-loop thermal correction is expressed as

$V_{T}\left(\phi,T\right)=\sum_{i}(-1)^{F_{b/f}}g_{i}\frac{T^{4}}{2\pi^{2}}J_{b/f}\left[\frac{m_{i}^{2}\left(\phi\right)}{T^{2}}\right]$
(34)

where

$J_{b/f}\left(y^{2}\right)=\int_{0}^{\infty}dx\,x^{2}\ln\left[1\mp
e^{-\sqrt{x^{2}+y^{2}}}\right]$ (35)

and these thermal functions in the high-temperature limit are given by

$\displaystyle J_{b}\left(\frac{m^{2}}{T^{2}}\right)$
$\displaystyle=-\frac{\pi^{4}}{45}+\frac{\pi^{2}}{12}\left(\frac{m}{T}\right)^{2}-\frac{\pi}{6}\left(\frac{m^{2}}{T^{2}}\right)^{3/2}$
$\displaystyle-\frac{1}{32}\left(\frac{m}{T}\right)^{4}\ln\left(\frac{m^{2}}{a_{b}T^{2}}\right)+\cdots$
(36)
$J_{f}\left(\frac{m^{2}}{T^{2}}\right)=\frac{7\pi^{4}}{360}-\frac{\pi^{2}}{24}\left(\frac{m}{T}\right)^{2}-\frac{1}{32}\left(\frac{m}{T}\right)^{4}\ln\left(\frac{m^{2}}{a_{f}T^{2}}\right)+\cdots$
(37)

where $\ln(a_{b})=5.4076$ and $\ln(a_{f})=2.6351$. Thus, in the high-
temperature limit one has

$V_{T}\left(\phi,T\right)=\frac{1}{24}\left(\kappa_{u}+\kappa_{u}^{\prime}\right)T^{2}\phi^{2}+\cdots.$
(38)

## References

  * Agazie  _et al._ (2023) G. Agazie _et al._ (NANOGrav), Astrophys. J. Lett. 951, L8 (2023), arXiv:2306.16213 [astro-ph.HE] .
  * Antoniadis  _et al._ (2023a) J. Antoniadis _et al._ , (2023a), arXiv:2306.16214 [astro-ph.HE] .
  * Reardon  _et al._ (2023) D. J. Reardon _et al._ , Astrophys. J. Lett. 951, L6 (2023), arXiv:2306.16215 [astro-ph.HE] .
  * Xu  _et al._ (2023) H. Xu _et al._ , Res. Astron. Astrophys. 23, 075024 (2023), arXiv:2306.16216 [astro-ph.HE] .
  * Hellings and Downs (1983) R. w. Hellings and G. s. Downs, Astrophys. J. Lett. 265, L39 (1983).
  * Afzal  _et al._ (2023) A. Afzal _et al._ (NANOGrav), Astrophys. J. Lett. 951, L11 (2023), arXiv:2306.16219 [astro-ph.HE] .
  * Antoniadis  _et al._ (2023b) J. Antoniadis _et al._ , (2023b), arXiv:2306.16227 [astro-ph.CO] .
  * Rajagopal and Romani (1995) M. Rajagopal and R. W. Romani, Astrophys. J. 446, 543 (1995), arXiv:astro-ph/9412038 .
  * Ellis  _et al._ (2023a) J. Ellis, M. Fairbairn, G. Hütsi, J. Raidal, J. Urrutia, V. Vaskonen, and H. Veermäe, (2023a), arXiv:2306.17021 [astro-ph.CO] .
  * Bi  _et al._ (2023) Y.-C. Bi, Y.-M. Wu, Z.-C. Chen, and Q.-G. Huang, (2023), arXiv:2307.00722 [astro-ph.CO] .
  * Madge  _et al._ (2023) E. Madge, E. Morgante, C. Puchades-Ibáñez, N. Ramberg, W. Ratzinger, S. Schenk, and P. Schwaller, (2023), arXiv:2306.14856 [hep-ph] .
  * Ellis  _et al._ (2023b) J. Ellis, M. Lewicki, C. Lin, and V. Vaskonen, (2023b), arXiv:2306.17147 [astro-ph.CO] .
  * Ratzinger and Schwaller (2021) W. Ratzinger and P. Schwaller, SciPost Phys. 10, 047 (2021), arXiv:2009.11875 [astro-ph.CO] .
  * Addazi  _et al._ (2023) A. Addazi, Y.-F. Cai, A. Marciano, and L. Visinelli, (2023), arXiv:2306.17205 [astro-ph.CO] .
  * Ghosh  _et al._ (2023) T. Ghosh, A. Ghoshal, H.-K. Guo, F. Hajkarim, S. F. King, K. Sinha, X. Wang, and G. White, (2023), arXiv:2307.02259 [astro-ph.HE] .
  * Rezapour  _et al._ (2022a) S. Rezapour, K. Bitaghsir Fadafan, and M. Ahmadvand, Annals Phys. 437, 168731 (2022a), arXiv:2006.04265 [hep-th] .
  * Rezapour  _et al._ (2022b) S. Rezapour, K. Bitaghsir Fadafan, and M. Ahmadvand, Phys. Scripta 97, 035301 (2022b).
  * Neronov  _et al._ (2021) A. Neronov, A. Roper Pol, C. Caprini, and D. Semikoz, Phys. Rev. D 103, 041302 (2021), arXiv:2009.14174 [astro-ph.CO] .
  * Ahmadvand and Bitaghsir Fadafan (2018) M. Ahmadvand and K. Bitaghsir Fadafan, Phys. Lett. B 779, 1 (2018), arXiv:1707.05068 [hep-th] .
  * Ahmadvand and Bitaghsir Fadafan (2017) M. Ahmadvand and K. Bitaghsir Fadafan, Phys. Lett. B 772, 747 (2017), arXiv:1703.02801 [hep-th] .
  * Athron  _et al._ (2023) P. Athron, A. Fowlie, C.-T. Lu, L. Morris, L. Wu, Y. Wu, and Z. Xu, (2023), arXiv:2306.17239 [hep-ph] .
  * Han  _et al._ (2023) C. Han, K.-P. Xie, J. M. Yang, and M. Zhang, (2023), arXiv:2306.16966 [hep-ph] .
  * Li and Xie (2023) S.-P. Li and K.-P. Xie, (2023), arXiv:2307.01086 [hep-ph] .
  * Blasi  _et al._ (2021) S. Blasi, V. Brdar, and K. Schmitz, Phys. Rev. Lett. 126, 041305 (2021), arXiv:2009.06607 [astro-ph.CO] .
  * Ellis and Lewicki (2021) J. Ellis and M. Lewicki, Phys. Rev. Lett. 126, 041304 (2021), arXiv:2009.06555 [astro-ph.CO] .
  * Buchmuller  _et al._ (2020) W. Buchmuller, V. Domcke, and K. Schmitz, Phys. Lett. B 811, 135914 (2020), arXiv:2009.10649 [astro-ph.CO] .
  * Lazarides  _et al._ (2023) G. Lazarides, R. Maji, and Q. Shafi, (2023), arXiv:2306.17788 [hep-ph] .
  * Wang  _et al._ (2023) Z. Wang, L. Lei, H. Jiao, L. Feng, and Y.-Z. Fan, (2023), arXiv:2306.17150 [astro-ph.HE] .
  * Bian  _et al._ (2022a) L. Bian, J. Shu, B. Wang, Q. Yuan, and J. Zong, Phys. Rev. D 106, L101301 (2022a), arXiv:2205.07293 [hep-ph] .
  * Bai  _et al._ (2023) Y. Bai, T.-K. Chen, and M. Korwar, (2023), arXiv:2306.17160 [hep-ph] .
  * Kitajima  _et al._ (2023) N. Kitajima, J. Lee, K. Murai, F. Takahashi, and W. Yin, (2023), arXiv:2306.17146 [hep-ph] .
  * Blasi  _et al._ (2023) S. Blasi, A. Mariotti, A. Rase, and A. Sevrin, (2023), arXiv:2306.17830 [hep-ph] .
  * Li  _et al._ (2023) Y. Li, L. Bian, and Y. Jia, (2023), arXiv:2304.05220 [hep-ph] .
  * Bian  _et al._ (2022b) L. Bian, S. Ge, C. Li, J. Shu, and J. Zong, (2022b), arXiv:2212.07871 [hep-ph] .
  * Babichev  _et al._ (2023) E. Babichev, D. Gorbunov, S. Ramazanov, R. Samanta, and A. Vikman, (2023), arXiv:2307.04582 [hep-ph] .
  * Sakharov  _et al._ (2021) A. S. Sakharov, Y. N. Eroshenko, and S. G. Rubin, Phys. Rev. D 104, 043005 (2021), arXiv:2104.08750 [hep-ph] .
  * Zhang  _et al._ (2023) Z. Zhang, C. Cai, Y.-H. Su, S. Wang, Z.-H. Yu, and H.-H. Zhang, (2023), arXiv:2307.11495 [hep-ph] .
  * Depta  _et al._ (2023) P. F. Depta, K. Schmidt-Hoberg, and C. Tasillo, (2023), arXiv:2306.17836 [astro-ph.CO] .
  * Gouttenoire  _et al._ (2023) Y. Gouttenoire, S. Trifinopoulos, G. Valogiannis, and M. Vanvlasselaer, (2023), arXiv:2307.01457 [astro-ph.CO] .
  * Franciolini  _et al._ (2023) G. Franciolini, A. Iovino, Junior., V. Vaskonen, and H. Veermae, (2023), arXiv:2306.17149 [astro-ph.CO] .
  * Guo  _et al._ (2023) S.-Y. Guo, M. Khlopov, X. Liu, L. Wu, Y. Wu, and B. Zhu, (2023), arXiv:2306.17022 [hep-ph] .
  * Hosseini Mansoori  _et al._ (2023) S. A. Hosseini Mansoori, F. Felegray, A. Talebian, and M. Sami, (2023), arXiv:2307.06757 [astro-ph.CO] .
  * Firouzjahi and Riotto (2023) H. Firouzjahi and A. Riotto, (2023), arXiv:2304.07801 [astro-ph.CO] .
  * Vagnozzi (2021) S. Vagnozzi, Mon. Not. Roy. Astron. Soc. 502, L11 (2021), arXiv:2009.13432 [astro-ph.CO] .
  * Ashoorioon  _et al._ (2022) A. Ashoorioon, K. Rezazadeh, and A. Rostami, Phys. Lett. B 835, 137542 (2022), arXiv:2202.01131 [astro-ph.CO] .
  * Niu and Rahat (2023) X. Niu and M. H. Rahat, (2023), arXiv:2307.01192 [hep-ph] .
  * Servant and Simakachorn (2023) G. Servant and P. Simakachorn, (2023), arXiv:2307.03121 [hep-ph] .
  * Bian  _et al._ (2023) L. Bian, S. Ge, J. Shu, B. Wang, X.-Y. Yang, and J. Zong, (2023), arXiv:2307.02376 [astro-ph.HE] .
  * Ellis  _et al._ (2023c) J. Ellis, M. Fairbairn, G. Franciolini, G. Hütsi, A. Iovino, M. Lewicki, M. Raidal, J. Urrutia, V. Vaskonen, and H. Veermäe, (2023c), arXiv:2308.08546 [astro-ph.CO] .
  * Rubakov (1997) V. A. Rubakov, JETP Lett. 65, 621 (1997), arXiv:hep-ph/9703409 .
  * Di Luzio  _et al._ (2020) L. Di Luzio, M. Giannotti, E. Nardi, and L. Visinelli, Phys. Rept. 870, 1 (2020), arXiv:2003.01100 [hep-ph] .
  * Ellis  _et al._ (2019) J. Ellis, M. Lewicki, J. M. No, and V. Vaskonen, JCAP 06, 024 (2019), arXiv:1903.09642 [hep-ph] .
  * Brdar  _et al._ (2019) V. Brdar, A. J. Helmboldt, and J. Kubo, JCAP 02, 021 (2019), arXiv:1810.12306 [hep-ph] .
  * Kobakhidze  _et al._ (2017) A. Kobakhidze, C. Lagger, A. Manning, and J. Yue, Eur. Phys. J. C 77, 570 (2017), arXiv:1703.06552 [hep-ph] .
  * Von Harling  _et al._ (2020) B. Von Harling, A. Pomarol, O. Pujolàs, and F. Rompineve, JHEP 04, 195 (2020), arXiv:1912.07587 [hep-ph] .
  * Yang  _et al._ (2023) A. Yang, J. Ma, S. Jiang, and F. P. Huang, (2023), arXiv:2306.17827 [hep-ph] .
  * Peccei and Quinn (1977) R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 (1977).
  * Berezhiani  _et al._ (2001) Z. Berezhiani, L. Gianfagna, and M. Giannotti, Phys. Lett. B 500, 286 (2001), arXiv:hep-ph/0009290 .
  * Dimopoulos  _et al._ (2016) S. Dimopoulos, A. Hook, J. Huang, and G. Marques-Tavares, JHEP 11, 052 (2016), arXiv:1606.03097 [hep-ph] .
  * Zhitnitsky (1980) A. R. Zhitnitsky, Sov. J. Nucl. Phys. 31, 260 (1980).
  * Dine  _et al._ (1981) M. Dine, W. Fischler, and M. Srednicki, Phys. Lett. B 104, 199 (1981).
  * Kim (1979) J. E. Kim, Phys. Rev. Lett. 43, 103 (1979).
  * Shifman  _et al._ (1980) M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Nucl. Phys. B 166, 493 (1980).
  * Georgi  _et al._ (1981) H. M. Georgi, L. J. Hall, and M. B. Wise, Nucl. Phys. B 192, 409 (1981).
  * Kamionkowski and March-Russell (1992) M. Kamionkowski and J. March-Russell, Phys. Lett. B 282, 137 (1992), arXiv:hep-th/9202003 .
  * Holman  _et al._ (1992) R. Holman, S. D. H. Hsu, T. W. Kephart, E. W. Kolb, R. Watkins, and L. M. Widrow, Phys. Lett. B 282, 132 (1992), arXiv:hep-ph/9203206 .
  * Baker  _et al._ (2006) C. A. Baker _et al._ , Phys. Rev. Lett. 97, 131801 (2006), arXiv:hep-ex/0602020 .
  * Peccei (2008) R. D. Peccei, Lect. Notes Phys. 741, 3 (2008), arXiv:hep-ph/0607268 .
  * Raffelt (2008) G. G. Raffelt, Lect. Notes Phys. 741, 51 (2008), arXiv:hep-ph/0611350 .
  * Arvanitaki  _et al._ (2010) A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, Phys. Rev. D 81, 123530 (2010), arXiv:0905.4720 [hep-th] .
  * Shakeri and Hajkarim (2023) S. Shakeri and F. Hajkarim, JCAP 04, 017 (2023), arXiv:2209.13572 [hep-ph] .
  * Gherghetta  _et al._ (2016) T. Gherghetta, N. Nagata, and M. Shifman, Phys. Rev. D 93, 115010 (2016), arXiv:1604.01127 [hep-ph] .
  * Fukuda  _et al._ (2015) H. Fukuda, K. Harigaya, M. Ibe, and T. T. Yanagida, Phys. Rev. D 92, 015021 (2015), arXiv:1504.06084 [hep-ph] .
  * Coleman and Weinberg (1973) S. R. Coleman and E. J. Weinberg, Phys. Rev. D 7, 1888 (1973).
  * Quiros (1999) M. Quiros, in  _ICTP Summer School in High-Energy Physics and Cosmology_ (1999) pp. 187–259, arXiv:hep-ph/9901312 .
  * Gildener and Weinberg (1976) E. Gildener and S. Weinberg, Phys. Rev. D 13, 3333 (1976).
  * Linde (1983) A. D. Linde, Nucl. Phys. B 216, 421 (1983), [Erratum: Nucl.Phys.B 223, 544 (1983)].
  * Arnold and Vokos (1991) P. B. Arnold and S. Vokos, Phys. Rev. D 44, 3620 (1991).
  * Witten (1981) E. Witten, Nucl. Phys. B 177, 477 (1981).
  * Lewicki and Vaskonen (2023) M. Lewicki and V. Vaskonen, Eur. Phys. J. C 83, 109 (2023), arXiv:2208.11697 [astro-ph.CO] .
  * Amaro-Seoane (2017) e. a. Amaro-Seoane, arXiv e-prints , arXiv:1702.00786 (2017), arXiv:1702.00786 [astro-ph.IM] .
  * et al. (2011) K. S. et al., Classical and Quantum Gravity 28, 094011 (2011).
  * Crowder and Cornish (2005) J. Crowder and N. J. Cornish, Phys. Rev. D 72, 083005 (2005), arXiv:gr-qc/0506015 .
  * Punturo  _et al._ (2010) M. Punturo _et al._ , Class. Quant. Grav. 27, 194002 (2010).
  * Reitze  _et al._ (2019) D. Reitze _et al._ , Bull. Am. Astron. Soc. 51, 035 (2019), arXiv:1907.04833 [astro-ph.IM] .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Dunsky  _et al._ (2022) D. I. Dunsky, L. J. Hall, and K. Harigaya, (2022), arXiv:2205.11540 [hep-ph] .
  * Ahmadvand (2021) M. Ahmadvand, JHEP 10, 109 (2021), arXiv:2108.00958 [hep-ph] .
  * Raffelt (1990) G. G. Raffelt, Phys. Rept. 198, 1 (1990).
  * Essig  _et al._ (2010) R. Essig, R. Harnik, J. Kaplan, and N. Toro, Phys. Rev. D 82, 113008 (2010), arXiv:1008.0636 [hep-ph] .
  * Jaeckel  _et al._ (2018) J. Jaeckel, P. C. Malta, and J. Redondo, Phys. Rev. D 98, 055032 (2018), arXiv:1702.02964 [hep-ph] .
  * Bergsma  _et al._ (1985) F. Bergsma _et al._ (CHARM), Phys. Lett. B 157, 458 (1985).
  * Aloni  _et al._ (2019) D. Aloni, Y. Soreq, and M. Williams, Phys. Rev. Lett. 123, 031803 (2019), arXiv:1811.03474 [hep-ph] .
  * Alves and Weiner (2018) D. S. M. Alves and N. Weiner, JHEP 07, 092 (2018), arXiv:1710.03764 [hep-ph] .
  * Ahn  _et al._ (2019) J. K. Ahn _et al._ (KOTO), Phys. Rev. Lett. 122, 021802 (2019), arXiv:1810.09655 [hep-ex] .
  * Knapen  _et al._ (2017) S. Knapen, T. Lin, H. K. Lou, and T. Melia, Phys. Rev. Lett. 118, 171801 (2017), arXiv:1607.06083 [hep-ph] .
  * Bertholet  _et al._ (2022) E. Bertholet, S. Chakraborty, V. Loladze, T. Okui, A. Soffer, and K. Tobioka, Phys. Rev. D 105, L071701 (2022), arXiv:2108.10331 [hep-ph] .
  * Aad  _et al._ (2022) G. Aad _et al._ (ATLAS), JHEP 08, 104 (2022), arXiv:2202.07953 [hep-ex] .
  * Tumasyan  _et al._ (2022) A. Tumasyan _et al._ (CMS), Phys. Rev. D 105, 092007 (2022), arXiv:2201.11585 [hep-ex] .