# LIGO as a probe of Dark Sectors

Fei Huanga,b<EMAIL_ADDRESS>Veronica Sanzc,d<EMAIL_ADDRESS>Jing
Shua,e,f,g,h,i<EMAIL_ADDRESS>Xiao Xuea,e<EMAIL_ADDRESS>a CAS Key
Laboratory of Theoretical Physics, Institute of Theoretical Physics, Chinese
Academy of Sciences, Beijing 100190, China  b Department of Physics and
Astronomy, University of California, Irvine, CA 92697 USA  c Instituto de
Física Corpuscular (IFIC), Universidad de Valencia-CSIC, E-46980 Valencia,
Spain  d Department of Physics and Astronomy, University of Sussex, Brighton
BN1 9QH, UK  e School of Physical Sciences, University of Chinese Academy of
Sciences, Beijing 100049, P. R. China  f CAS Center for Excellence in Particle
Physics, Beijing 100049, China  g Center for High Energy Physics, Peking
University, Beijing 100871, China  h School of Fundamental Physics and
Mathematical Sciences, Hangzhou Institute for Advanced Study, UCAS, Hangzhou
310024, China  i International Centre for Theoretical Physics Asia-Pacific,
Beijing/Hangzhou, China

###### Abstract

We show how current LIGO data is able to probe interesting theories beyond the
Standard Model, particularly Dark Sectors where a Dark Higgs triggers symmetry
breaking via a first-order phase transition. We use publicly available LIGO O2
data to illustrate how these sectors, even if disconnected from the Standard
Model, can be probed by Gravitational Wave detectors. We link the LIGO
measurements with the model content and mass scale of the Dark Sector, finding
that current O2 data is testing a broad set of scenarios where the breaking of
$SU(N)$ theories with $N_{f}$ fermions is triggered by a Dark Higgs at scales
$\Lambda\simeq 10^{8}-10^{9}$ GeV with reasonable parameters for the scalar
potential.

††preprint: UCI-HEP-TR-2021-06

##  I Introduction

Much of the Universe is dark, and many theories have been built trying to
explain it. Our hopes for probing these theories often rely on their possible
connection to regular matter via some form of non-gravitational interaction
For example, direct searches for Dark Matter hinges on some sort of coupling
to nucleons or electrons, and constraints on those couplings usually assume a
mechanism of communication between the Dark Sector and the rest of the
Universe which establishes some form of tracking between these two sectors.

The first observation of Gravitational Waves (GW) by the LIGO and Virgo
collaborations Abbott  _et al._ (2016a); Aasi  _et al._ (2015); Acernese  _et
al._ (2014) in 2015 initiated a new way to see the Universe, and since then
exciting new observations have provided information about astrophysical
objects like Black Holes Abbott  _et al._ (2019a, 2020). However, the physics
reach for LIGO is not circumscribed to detection of mergers Allen and Romano
(1999); Maggiore (2000); Regimbau (2011); Caprini and Figueroa (2018). The
detection, or the lack of, a stochastic GW background allows us to explore
interesting, non-standard sectors. We will explain how, with the current
public data from LIGO, one can probe plausible Dark Sector scenarios,
regardless of their non-gravitational interaction with visible matter.

These Dark Sectors could resemble Standard-Model dynamics, with new forces,
Dark Higgses and states charged under them. Influenced by thermal
contributions from the degrees of freedom in the Dark Sector, the thermal
history of the Dark Higgs could then lead to first-order phase transitions.
Many studies have been devoted to the prospects that future interferometers
could offer to explore Dark Sectors, e.g., Ref Croon  _et al._ (2018). In this
paper we explore the possibilities that LIGO and its current public dataset
present, and bridge the gap between generic studies of thermal parameters,
e.g., Ref Romero  _et al._ (2021), and specific particle-physics models.

The paper is structured as follows. In Sec. II we first describe the analysis
of GWs from first-order phase transitions, then discuss in Sec. III the
connection between the phase transition thermal parameters with classes of
particle-physics models, especially of $SU(N)/SU(N-1)$ Dark Sectors. We
finally link these models to current exclusions set by LIGO, and in Sec. IV we
conclude the discussion.

##  II Gravitational Waves from Phase Transitions

The stochastic gravitational wave background (SGWB) is often considered as an
isotropic, unpolarized, stationary and Gaussian background generated by a
large number of unresolved gravitational-wave sources Allen and Romano (1999);
Maggiore (2000); Regimbau (2011); Caprini and Figueroa (2018). Its power
spectrum is characterized by the dimensionless quantity

$\Omega_{\rm GW}(f)\equiv\frac{1}{\rho_{\rm crit}}\frac{d\rho_{\rm GW}}{d\ln
f},$ (1)

where $\rho_{\rm GW}$ is the energy density of the stochastic gravitational
wave background, $f$ is the frequency of the GW and

$\rho_{\rm crit}\equiv\frac{3c^{2}H_{0}^{2}}{8\pi G}$ (2)

is the critical energy density of the universe today.

In principle, the total SGWB is a superposition of all possible astrophysical
and cosmological sources. However, we can obtain a conservative upper limit
for the SGWB from phase transitions that occur in the early universe by
assuming that phase transition dynamics is the main source of SGWB.

Figure 1: The coloured curves show the spectral shape $S_{\rm sw}(f)$ by
varying the peak frequency $f_{\rm sw}$, whereas the dashed and dotted black
lines indicate its asymptotic behaviour. The two solid vertical lines at $20$
and $1726~{}\rm Hz$ show the minimum and maximum frequency considered when
obtaining the GW upper limits using LIGO data.

The SGWB generated from phase transitions in the early universe consists of
three parts Caprini  _et al._ (2008); Huber and Konstandin (2008); Caprini
_et al._ (2009); Espinosa  _et al._ (2010); Hindmarsh  _et al._ (2014);
Hindmarsh (2018):

$\Omega_{\rm GW}=\Omega_{\rm col}+\Omega_{\rm sw}+\Omega_{\rm turb}\,,$ (3)

in which the three terms on the right hand side correspond to the contribution
from bubble collisions, sound waves in the fluid and the turbulence,
respectively. For simplicity, we shall assume in this work that contributions
from sound waves are always dominant. We emphasize that this is typically the
case for models in which gauge bosons acquire masses during the phase
transition Bodeker and Moore (2017). However, the analysis we present can be
easily generalized to cases in which other types of contribution become more
important.

The phase transition is in general characterized by just a few parameters: the
velocity of the bubble wall $v_{w}$, the ratio of the free energy density
difference between the true and false vacuum and the total energy density,
$\xi$, the speed of the phase transition $\beta/H$, and the nucleation
temperature $T_{N}$. With these parameters, the GW power spectrum can be
expressed as Weir (2018)

$\Omega_{\rm GW}h^{2}\simeq 8.5\times
10^{-6}\left(\frac{g_{*}}{100}\right)^{1/3}\Gamma^{2}\bar{U}_{f}^{4}\left(\frac{\beta}{H}\right)^{-1}v_{w}S_{\rm
sw}(f),$ (4)

where $g_{*}$ is the effective number of relativistic degrees of freedom at
the time of the transition, $\Gamma\sim 4/3$ is the adiabatic index,
$\bar{U}_{f}^{4}\sim(3/4)\kappa_{f}\xi$ is the root-mean-square fluid velocity
with the efficiency parameter given by the approximate expressions Espinosa
_et al._ (2010)

$\kappa_{f}\sim\begin{cases}\displaystyle\frac{\xi}{0.73+0.83\sqrt{\xi}+\xi}&v_{w}\rightarrow
1\\\ \displaystyle\frac{\xi^{2/5}}{0.017+(0.997+\xi)^{2/5}}&v_{w}\approx
0.5\end{cases}\,,$ (5)

The spectral shape $S_{\rm sw}$ is given by

$S_{\rm sw}(f)=\left(\frac{f}{f_{\rm
pk}}\right)^{3}\left(\frac{7}{4+3\left(f/f_{\rm pk}\right)^{2}}\right)^{7/2}$
(6)

with the peak frequency

$f_{\rm pk}=8.9\times
10^{-8}\mathrm{Hz}\left(\frac{1}{v_{w}}\right)\left(\frac{\beta}{H}\right)\left(\frac{T_{N}}{\rm
GeV}\right)\left(\frac{g_{*}}{100}\right)^{1/6}.$ (7)

The shape of $S_{\rm sw}$ is shown in FIG. 1, noting that $S_{\rm sw}$ is
equal to $\Omega_{\rm GW}(f)/\Omega_{\rm GW}(f_{pk})$. In this figure one
observes that varying the peak frequency $f_{\rm sw}$ amounts to simply
shifting the spectrum horizontally. Also note that the asymptotic behavior of
$S_{\rm sw}$ goes as $\sim f^{3}$ for $f\ll f_{\rm sw}$ (dotted line), whereas
one expects a behaviour $\sim f^{-4}$ for $f\gg f_{\rm sw}$ (dashed line).

The behavior of $S_{\rm sw}$ in the LIGO frequency range (indicated by the two
vertical lines) can therefore transition from a simple descending power law,
to one with a peak in between, and eventually to a ascending power law as we
increase $f_{\rm pk}$.

Figure 2: The upper limit of $\Omega_{\rm ref}$ for each different $f_{\rm
pk}$ using the $2\sigma$ criterion. The curve with triangles corresponds to
the upper limit at 25 Hz, while the curve with circles uses the upper limit at
$f_{\rm pk}$. The solid horizontal line corresponds to the BBN bound at
$f_{\rm pk}$.

With this spectrum, we follow the procedure laid out in Refs Abbott  _et al._
(2017, 2019b) to compute the upper limit of $\Omega_{\rm GW}$ for different
values of $f_{\rm pk}$ using data from LIGO O2 Abbott  _et al._ (2019b). More
details can be found in Appendix A. The only difference respect to these
references is related to the choice of the optimized estimator, i.e., the
estimator $\hat{\Omega}_{f_{\rm ref}}$ is optimized by putting $\Omega_{\rm
GW}(f)/\Omega_{\rm GW}(f_{\rm ref})$ instead of some power of $f/f_{\rm ref}$.
Note that this means the upper limit essentially depends only on the shape of
$S_{\rm sw}$ within the LIGO frequency band, since all the other factors drop
out when taking the ratio. The 95% confidence level upper limit is then
obtained by setting

$\Omega_{\rm GW}^{\rm up.~{}lim.}(f_{\rm ref})=2\sigma(f_{\rm ref})\,.$ (8)

In FIG. 2 we show the upper limit at $f_{\rm ref}=25$ Hz by the curve with
triangles. For $f_{\rm pk}\lesssim 10$ Hz or $f_{\rm pk}\gtrsim 10^{3}$ Hz,
the bound approaches a constant since the spectrum is essentially a power law
with fixed exponent within the LIGO frequency band. This behaviour allows us
to easily extrapolate the constraint to even smaller or bigger values of
$f_{\rm pk}$. However, for intermediate $f_{\rm pk}$ there is a smooth
transition on the upper limit between the ascending and descending asymptotic
behavior. We see that LIGO provides a stronger constraint on ascending spectra
than on the descending spectra, and the difference can be as large as an order
of magnitude.

In practice, LIGO could only provide reliable constraints on SGWB in the band
$20-1726$ Hz and we have chosen $f_{\rm ref}=25$ Hz as it is the frequency
where LIGO is most sensitive Abbott  _et al._ (2017, 2019b). However, the
thermal parameters of the phase transition are connected to the amplitude at
the $f_{\rm pk}$. In order to place constraint on the thermal parameters, we
notice that, for a given GW spectrum as in Eq. (6), the constraint at a
particular frequency in the LIGO band can be mapped to the peak frequency
$f_{\rm pk}$ using

$\Omega_{\rm GW}^{\rm up.~{}lim.}(f_{\rm pk})=\Omega_{\rm GW}^{\rm
up.~{}lim}(f_{\rm ref})\frac{S_{\rm sw}(f_{\rm pk})}{S_{\rm sw}(f_{\rm
ref})}\,.$ (9)

The curve made by circles in FIG. 2 shows the result of this mapping. For
$f_{\rm pk}$ very small or very large, the upper bound from LIGO becomes
substantially weaker, even when comparing with the constraint from BBN, which
is obtained by integrating out the spectrum and requiring $h^{2}\int d\ln
f~{}\Omega_{\rm GW}(f)<1\times 10^{-6}$ Caprini and Figueroa (2018); Tanin and
Tenkanen (2020).

In what follows, we shall discuss the LIGO constraints on the thermal
parameters and how such constraint can be utilized to constrain particular
models of phase transition.

##  III Scenarios for phase transitions and their thermal parameters

Typically, one represents phase transitions as driven by the dynamics of a
scalar field which transitions from one vacuum to another under the influence
of the evolving thermal potential. In that context, the thermal parameters
$\beta/H$ and $\xi$ are defined by

$\displaystyle\frac{\beta}{H}$ $\displaystyle=$
$\displaystyle\left.T\frac{(S_{E}/T)}{dT}\right|_{T=T_{N}}$ (10)
$\displaystyle\xi$ $\displaystyle=$
$\displaystyle\frac{1}{\rho_{N}}\left.\left(\Delta
V-T\Delta\frac{dV}{dT}\right)\right|_{T=T_{N}}$ (11)

in which $S_{E}$ is the Euclidean action, $V$ is the thermal potential of the
scalar field, and the bubble nucleation temperature can be obtained by solving

$\frac{S_{E}}{T_{N}}\approx 177-4\ln\left(\frac{T_{N}}{\rm GeV}\right)-2\ln
g_{\star}(T_{N})\,,\\\ $ (12)

in which $\rho_{N}=g_{\star}(T_{N})\pi^{2}T_{N}^{4}/30$. On the other hand,
the calculation of the bubble-wall velocity $v_{w}$ for a particular model is
highly non-trivial. Therefore, instead of directly calculating it, we shall
follow the customary convention of considering a few reference values,
$v_{w}=0.5$ and $v_{w}=1$.

To connect the thermal parameters, $\xi$ and $\beta/H$ to specific models we
will follow the approach described in Ref. Croon  _et al._ (2018). We will
consider classes of potentials which consists of competing terms with
alternating signs. Specifically, we will look into two types of finite-
temperature potentials

$\displaystyle V(H,T)$ $\displaystyle=$
$\displaystyle\frac{1}{2}m^{2}(T)h_{D}^{2}-c_{3}(T)h_{D}^{3}+\frac{1}{4}\lambda(T)h_{D}^{4}$
(13) $\displaystyle V(H,T)$ $\displaystyle=$
$\displaystyle\frac{1}{2}m^{2}(T)h_{D}^{2}-\frac{1}{4}\lambda(T)h_{D}^{4}+c_{6}(T)h_{D}^{6}$
(14)

in which the coefficients of the scalar field $h_{D}$ are all positive at the
time of transition.

Indeed, most dark phase transitions can be mapped onto these effective
scenarios. For example, in a Dark Sector where its particles acquire mass from
a Dark Higgs as its gauge group $SU(N)$ breaks into $SU(N-1)$, the potential
in Eq. (13) can be realized with renormalisable operators, whereas the
potential in Eq. (14) can be obtained from non-renormalisable sextet
interaction Croon  _et al._ (2018).

In what follows, we shall discuss the Dark Higgs - $SU(N)/SU(N-1)$ models with
renormalisable and non-renormalisable operators specifically.

###  III.1 Exploring Dark Sectors with LIGO: $SU(N)/SU(N-1)$ models

####  III.1.1 Models with renormalisable operators

For the type of potential in Eq. (13), we can parametrise zero-temperature
parameters as

$m^{2}(0)=-\frac{\Lambda^{4}}{v_{D}^{2}},\,\,\,\lambda(0)=\frac{\Lambda^{4}}{v_{D}^{4}}$
(15)

in which $v_{D}$ is the zero temperature vacuum expectation value, and
$\Lambda$ is the scale of the potential. With this parametrisation, the finite
temperature potential can be expressed as

$\displaystyle V(H,T)$ $\displaystyle=$
$\displaystyle\Lambda^{4}\left[-\frac{1}{2}\left(\frac{h_{D}}{v_{D}}\right)^{2}+\frac{1}{4}\left(\frac{h_{D}}{v_{D}}\right)^{4}\right]+\frac{T^{4}}{2\pi^{2}}\left[\sum_{i\in\rm
bosons}n_{i}J_{B}(m_{i}^{2}/T^{2})-\sum_{i\in\rm
fermions}n_{i}J_{F}(m_{i}^{2}/T^{2})\right]$ (16) $\displaystyle=$
$\displaystyle\Lambda^{4}\Bigg{\\{}\left[-\frac{1}{2}+\left(\frac{1}{8}+\frac{N_{G}}{24}\right)\frac{T^{2}}{v_{D}^{2}}+\frac{3}{24}N_{\rm
GB}\frac{g^{2}}{4}\frac{T^{2}v_{D}^{2}}{\Lambda^{4}}+y^{2}N_{f}\frac{T^{2}}{48}\frac{v_{D}^{2}}{\Lambda^{4}}\right]\left(\frac{h_{D}}{v_{D}}\right)^{2}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}-\left[N_{\rm
GB}\left(\frac{g^{2}}{4}\right)^{3/2}\frac{1}{4\pi}\frac{v_{D}^{3}T}{\Lambda^{4}}\right]\left(\frac{h_{D}}{v_{D}}\right)^{3}+\frac{1}{4}\left(\frac{h_{D}}{v_{D}}\right)^{4}\Bigg{\\}}\,,$

where $N_{\rm GB}=2N-1$ is the number of gauge bosons that couple to the Dark
Higgs with coupling constant $g$, and which get a mass from the Dark Higgs
interactions. $N_{G}=2N-1$ is the number of Goldstone bosons, and
$N_{f}=N\times N_{\rm FL}$ is the number of self-adjoint fermions with Yukawa
coupling $y$ where $N_{\rm FL}$ is the number of flavors. Note that for Dirac
fermions, one would need to double the number of degrees of freedom. For
simplicity, we assume the Yukawa coupling is universal for those fermions. See
Ref. Croon  _et al._ (2018) for more details.

In the second equality, the following high temperature expansions are used

$\displaystyle J_{B}(m^{2}/T^{2})\sim
2\pi^{2}\left(\frac{m^{2}}{24T^{2}}-\frac{m^{3}}{12\pi T^{3}}\right),$ (17)
$\displaystyle~{}~{}\text{and}~{}~{}J_{F}(m^{2}/T^{2})\sim-2\pi^{2}\left(\frac{m^{2}}{48T^{2}}\right),$
(18)

in which the field dependent masses can be read as

$\displaystyle m_{H}^{2}$ $\displaystyle=$
$\displaystyle\partial_{h_{D}}^{2}V=\Lambda^{4}\left(3\frac{h_{D}^{2}}{v^{4}}-\frac{1}{v^{2}}\right)\,,$
(19) $\displaystyle m_{G}^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{h_{D}}\partial_{h_{D}}V=\Lambda^{4}\left(\frac{h_{D}^{2}}{v_{D}^{4}}-\frac{1}{v_{D}^{2}}\right)\,,$
(20) $\displaystyle m_{\rm GB}$ $\displaystyle=$
$\displaystyle\frac{gh_{D}}{2}\,,$ (21) $\displaystyle m_{f}$ $\displaystyle=$
$\displaystyle\frac{yh_{D}}{\sqrt{2}}.$ (22)

Note that in the second line of Eq. (16), only the massive gauge bosons are
taken into account, i.e., the higher order term ($\sim m^{3}/T^{3}$) from the
Goldstone boson is neglected.Besides, the part of the expansion which gives
rise to terms independent of $h_{D}$ is neglected, since it only amounts to a
constant shift in the potential $V$. Finally, the mapping to the temperature
dependent parameters in Eq. (13) is straightforward by matching the terms with
the same powers of $h_{D}$.

Eq. (16) is also often written in the following form

$\displaystyle V(H,T)=\Lambda^{4}(T)\Bigg{[}$
$\displaystyle\left(\frac{3-4\alpha(T)}{2}\right)\left(\frac{h_{D}}{v_{D}(T)}\right)^{2}$
$\displaystyle-\left(\frac{h_{D}}{v_{D}(T)}\right)^{3}+\alpha(T)\left(\frac{h_{D}}{v_{D}(T)}\right)^{4}\Bigg{]}\,,$

in which the minimum of the potential can be easily obtained by minimizing the
potential

$v_{D}(T)=\frac{3c_{3}(T)+\sqrt{9c_{3}^{2}(T)-4m^{2}(T)\lambda^{2}(T)}}{2\lambda(T)}.$
(24)

By identifying

$\displaystyle
c_{3}(T)=\frac{\Lambda^{4}(T)}{v_{D}^{3}(T)},~{}~{}~{}~{}~{}\frac{\lambda(T)}{4}=\alpha(T)\frac{\Lambda^{4}(T)}{v_{D}^{4}(T)}\,,$
(25)

one finds

$\displaystyle\alpha(T)=\frac{\lambda(T)v_{D}(T)}{4c_{3}(T)},~{}\Lambda(T)=\left(c_{3}(T)v_{D}^{3}(T)\right)^{1/4}\,.$
(26)

With these, for $\alpha(T)\in[0.51,0.65]$, the effective action can be fitted
by Croon and White (2018)

$\frac{S_{E}}{T}=\frac{v_{D}^{3}(T)}{T\Lambda^{2}(T)}10^{a+b\frac{\left|\alpha(T)-0.75\right|^{c}}{\left|\alpha(T)-0.5\right|^{d}}}$
(27)

with the fitting parameters $a=-71.06,~{}b=71.62$, $c=0.008805$, and
$d=0.009263$.

####  III.1.2 Models with non-renormalisable operators

For the type of potential in Eq. (14), similar to the previous case, one can
perform the following parametrisation

$~{}m^{2}(0)=(2-3\alpha)\frac{\Lambda^{4}}{v_{D}^{2}},\,\,\,\lambda(0)=4\frac{\Lambda^{4}}{v_{D}^{4}},\,\,\,c_{6}(0)=\alpha\frac{\Lambda^{4}}{v_{D}^{6}}.$
(28)

The finite temperature potential then becomes

$\displaystyle V(H,T)$ $\displaystyle=$
$\displaystyle\Lambda^{4}\left[(2-3\alpha)\left(\frac{h_{D}}{v_{D}}\right)^{2}-\left(\frac{h_{D}}{v_{D}}\right)^{4}+\alpha\left(\frac{h_{D}}{v_{D}}\right)^{6}\right]$
(29) $\displaystyle+\frac{T^{4}}{2\pi^{2}}\left[\sum_{i\in\rm
bosons}n_{i}J_{B}(m_{i}^{2}/T^{2})-\sum_{i\in\rm
fermions}n_{i}J_{F}(m_{i}^{2}/T^{2})\right]$ $\displaystyle=$
$\displaystyle\Lambda^{4}\Bigg{\\{}\bigg{[}2-3\alpha-\left(\frac{1}{2}+\frac{N_{G}}{6}\right)\frac{T^{2}}{v_{D}^{2}}+\frac{3}{24}N_{\rm
GB}\frac{g^{2}}{4}\frac{T^{2}v_{D}^{2}}{\Lambda^{4}}+y^{2}N_{f}\frac{T^{2}}{48}\frac{v_{D}^{2}}{\Lambda^{4}}\bigg{]}\left(\frac{h_{D}}{v_{D}}\right)^{2}$
$\displaystyle~{}~{}~{}~{}~{}-\left[1-\frac{(30+6N_{G})\alpha
T^{2}}{24v_{D}^{2}}\right]\left(\frac{h_{D}}{v_{D}}\right)^{4}+\alpha\left(\frac{h_{D}}{v_{D}}\right)^{6}\Bigg{\\}}\,.$

Note that, in the high-temperature expansion, the cubic term is assumed to be
subdominant, i.e., we have only kept the part proportional to $m^{2}/T^{2}$
Croon  _et al._ (2018). Moreover, the field-dependent masses of the Goldstone
bosons and the Dark Higgs which goes into the thermal correction are

$\displaystyle m_{H}^{2}$ $\displaystyle=$
$\displaystyle\Lambda^{4}\left[\frac{2(2-3\alpha)}{v_{D}^{2}}-\frac{12h_{D}^{2}}{v_{D}^{4}}+\frac{30\alpha
h_{D}^{4}}{v_{D}^{6}}\right]\,,$ (30) $\displaystyle m_{G}^{2}$
$\displaystyle=$
$\displaystyle\Lambda^{4}\left[\frac{2(2-3\alpha)}{v_{D}^{2}}-\frac{4h_{D}^{2}}{v_{D}^{4}}+\frac{6\alpha
h_{D}^{4}}{v_{D}^{6}}\right]\,.$ (31)

The terms proportional to $\alpha h_{D}^{4}/v_{D}^{4}$ would give rise to the
thermal correction of the quartic term.

Just as we have done in the previous section, one can write Eq. (29) in terms
of temperature-dependent parameters:

$\displaystyle V(H,T)=\Lambda^{4}(T)\Bigg{[}$
$\displaystyle(2-3\alpha(T))\left(\frac{h}{v_{D}(T)}\right)^{2}-\left(\frac{h}{v_{D}(T)}\right)^{4}$
(32) $\displaystyle+\alpha(T)\left(\frac{h}{v_{D}(T)}\right)^{6}\Bigg{]}\,.$

The non-vanishing VEV

$v_{D}(T)=\left(\frac{\lambda(T)+\sqrt{\Lambda^{2}(T)-24c_{6}(T)m^{2}(T)}}{12c_{6}(T)}\right)^{1/2}$
(33)

is obtained by minimizing the potential. Suppose the non-vanishing VEV does
exists ($\Lambda^{2}(T)-24c_{6}(T)m^{2}(T)>0$), then $\alpha(T)$ can be
obtained by solving

$\lambda(T)=4\Lambda(T)^{4}/v(T)^{4}\,\,\,c_{6}(T)=\alpha(T)\Lambda(T)^{4}/v_{D}(T)^{6}.$
(34)

Therefore,

$\alpha(T)=4c_{6}(T)v_{D}^{2}(T)/\lambda(T).$ (35)

Following Croon  _et al._ (2018), the Euclidean action can be fitted by

$S_{E}=\frac{v_{D}^{3}(T)}{\Lambda^{2}(T)}10^{\sum_{i=1}^{3}a_{i}(\alpha(T)-2/3)^{i}}$
(36)

with $a_{i}=(-17.446,-132.404,-763.744)$ for $\alpha(T)\in[0.51,0.65]$.

Finally, for $g_{\star}(T_{N})$, at high temperatures, i.e.,
$T\gg\mathcal{O}(100)~{}\text{GeV}\gg m_{i},~{}m_{h_{D}}$, we shall assume
that the particles in the Dark Sector are the only degrees of freedom in
addition to the SM. Therefore, the effective number of relativistic degrees of
freedom will be given by

$g_{\star}(T_{N})\approx 106.75+3N_{\rm GB}+2(N^{2}-1-N_{\rm
GB})+1+\frac{7}{8}\times 2\times N_{f}\,,$ (37)

where we have included all degrees of freedom from the SM, as well as the
massless and massive gauge bosons and the dark fermions charged under $SU(N)$.

####  III.1.3 Results

Figure 3: 95% confidence level exclusion region from LIGO (curved contours)
and the constraint from BBN (vertical contours), projected on the
$T_{N}-\beta/H$ plane. The values of $\xi$ corresponding to this exclusion are
represented by the colour variation. The constraints are calculated for two
reference bubble-wall velocity $v_{w}=0.5$ and $v_{w}=1$. Markers with
different shapes are obtained from different models of phase transitions. The
top two panels correspond to the cases with renormalisable operators, whereas
the bottom two panels show the cases with non-renormalisable operators.

In FIG. 3, we present the LIGO constraints on the thermal parameters at 95%
confidence level together with the constraint from BBN. In the left panels,
the constraints are evaluated with $v_{w}=0.5$, whereas the bounds in the
right panels are evaluated at $v_{w}=1$.

Note that, with $v_{w}$ fixed, the thermal parameters $(\beta/H,T_{N},\xi)$
constitutes a 3-dimensional space. In our analysis, we find the the excluded
region can be conveniently projected on the 2-dimensional plane of $\beta/H$
and $T_{N}$ as the contours of $\xi$ do not intersect each other. This thus
enables us to use colour variation to represent the values of $\xi$ on the
exclusion contours. Indeed, within a particular contour, any value of $\xi$
larger than the value associated with the contour is excluded. As a result,
the regions in the 2D plots where the LIGO ellipsoidal contours have a colour
lighter than the BBN vertical contours suggest that LIGO can set a better
bound than BBN, and viceversa.

On the same figure, we also plot the values of the thermal parameters obtained
from different models defined by $\Lambda,~{}N$ and $N_{\rm FL}$. We show
examples with $N=2$ and 5, and a fixed number of flavours $N_{\rm FL}=3$, as
representative of the behaviour expected from the choices on matter content
and gauge groups. We also fix $g=1$ and $y=1$, although similar behaviour is
found for similar ${\cal O}(1)$ values.

If a set of thermal parameters obtained from a particular model lies within a
contour, and has a larger $\xi$ (darker colour) than the corresponding
contour, then this set of thermal parameter is excluded at 95% CL.

From the top figures 3, one can deduce that the constraint from LIGO is still
unable to probe the renormalisable classes of models, as the markers
representing the thermal parameters barely touch the excluded region.

However, in the non-renormalisable case, we indeed see how LIGO is able to
probe these Dark Sectors. In particular we observe that the dark scales
$\Lambda\sim 10^{8}$ \- $10^{9}$ GeV are best constrained by LIGO. Note that
models with larger group rank $N$ and number of flavours $N_{\rm FL}$ are more
constrained as they tend to produce a larger $\xi$.

The regions excluded by LIGO O2 correspond to $\Lambda$ in the region around
$10^{8}$-$10^{9}$ GeV. In this region we have varied $v_{w}$ from 0.5 and 1
and explored different options for $N$ and $N_{\rm FL}$ as illustrated by FIG.
3. The range of parameters $\alpha$ and $x\equiv v_{D}/\Lambda$ where LIGO
currently has sensitivity is given by

$\alpha\in[0.6,0.65]\,,\,x=\frac{v_{D}}{\Lambda}\in[1.4,2]\,.$ (38)

These regions can be translated into parameters in the scalar potential Eq.
(14) by inspecting Eq. (28):

$\displaystyle\frac{m}{\Lambda}\in[0.1,0.35],\,\lambda\in[0.2-1.1],\,c_{6}\Lambda^{2}\in[0.01,0.1].\,$
(39)

Clearly, these are reasonable choices of parameter space and indicate that
LIGO is testing interesting Dark-Sector theories.

##  IV Conclusions

When discussing Dark Sectors and their GW signatures, we usually think on
future probes like LISA, many years from now. Here we have shown that LIGO is
already probing interesting scenarios for Dark Sectors.

In the context of first-order phase transitions and LIGO data, the emphasis
has been placed in performing effective analyses, such as Ref. Romero  _et
al._ (2021) where the authors explore the bounds from LIGO O3 using
parametrisations of the power spectrum. In this paper, we take a complementary
step and focus on examining whether concrete particle-physics models could be
related to the tested regions.

In this work we answer the question whether the type of first-order phase
transitions LIGO is currently probing could be represented by concrete,
reasonable particle-physics models. For that reason, we have focused on
classes of models which capture a broad set of features of Dark Sectors, and,
at the same time, capable of producing interesting GW signatures. The
renormalisable and non-renormalisable benchmarks we used had been identified
in Ref. Croon  _et al._ (2018) as more promising for strong first-order phase
transitions from Dark Sectors.

We choose to set up those models with a breaking $SU(N)/SU(N-1)$ which should
be understood as an example in which some bosonic and fermionic degrees of
freedom influence the thermal history of the Dark Higgs. Such a choice also
allows a simple parametrisation in terms of the group rank $N$ and the number
of flavours $N_{\rm FL}$. We find that scales around $10^{8}-10^{9}$ GeV are
better probed by LIGO O2. We also find that the sensitive regions correspond
to moderate values for $N$ and $N_{\rm FL}$, evidencing that LIGO is not
testing extreme regions in the UV parameter space. Of course, various other
types of models for phase transition could also generate GW whose spectrum
lies in the frequency range relevant for LIGO, e.g., models motivated by grand
unification theories Croon  _et al._ (2019a, b) and models for confinement-
deconfinement phase transition Huang  _et al._ (2020). It is straightforward
to see whether LIGO might be able to constrain those models once the thermal
parameters are computed.

Note that, when obtaining the GW peak amplitude in Eq. (4), we have used the
standard formula from Ref. Weir (2018). Recent discussions in Ref. Guo  _et
al._ (2021) suggest the existence of an additional suppression factor due to
the finite lifetime of the sound waves. Moreover, it has been shown recently
that theoretical uncertainties in GW production could lead to changes in the
GW spectrum as large as several orders of magnitude Croon  _et al._ (2020). In
addition, we have assumed in our analysis that sound waves are the dominant
source for GW production. In other scenarios, for example, strongly
supercooled phase transitions Ellis  _et al._ (2019, 2020); Lewicki and
Vaskonen (2020a, b), or scenarios in which new heavy particles can provide
sufficiently large friction Azatov and Vanvlasselaer (2021), contribution from
turbulence or bubble collision could also be important and one needs to care
about their effects in the GW power spectrum. Although subject to those
uncertainties, our analysis nevertheless continues to provide a concrete
method to constrain Dark-Sector models with LIGO data.

Our analysis of the SGWB is based on publicly available O2 data. In Appendix
A, we have provided explanations on how to reproduce our analysis. Recent
papers from authors in the LIGO/Virgo collaboration, e.g., Romero  _et al._
(2021); Abbott  _et al._ (2021), make use of the O3 data. GW constraints on
particle-physics models considered in this paper are expected to improve as
more data becomes available.

We believe our results motivate a more systematic study of the particle-
physics scenarios that the LIGO experiment is able to test. We emphasize again
that traditional direct or indirect searches for dark particles assume that
the Dark Sector interacts non-gravitationally with the Standard Model. On the
other hand, since gravity is universal, methods for probing the Dark-Sector
via its gravitational effects such as structure formation (see Ref. Dienes
_et al._ (2020, 2021) and references in it for recent progress) and
gravitational waves do not rely on those assumptions. These gravitational
effects therefore offer unique opportunities to access Dark Sectors, which
would otherwise be hidden from us if they lack a connection to the Standard
Model.

## Acknowledgements

We would like to thank Djuna Croon for conversations at the beginning of this
project. V.S. acknowledges support from the UK Science and Technology
Facilities Council ST/L000504/1. J.S. and F.H. are supported by the National
Natural Science Foundation of China under Grants No. 12025507, No. 11690022,
No.11947302; and is supported by the Strategic Priority Research Program and
Key Research Program of Frontier Science of the Chinese Academy of Sciences
under Grants No. XDB21010200, No. XDB23010000, and No. ZDBS-LY-7003. F.H. is
also supported by the National Science Foundation of China under Grants No.
12022514 and No. 11875003.

##  Appendix A On the use of LIGO data

In this Appendix, we describe what we actually do with data downloaded from
LIGO. This consists of two parts: 1) the data selection in which we select
data that satisfy certain criterion, and 2) the analysis in which we estimate
the GW upper limit using data selected in the first part.

###  A.1 Data Selection

To perform this analysis, one needs to download the data of the LIGO detectors
at both Hanford and Livingston from O2 data release (https://www.gw-
openscience.org/data/) Abbott  _et al._ (2019b) with a 16 kHz sampling
frequency. Each file covers a 4096 s period of measurement, and the file name
contains the start time of the measurement, which is referred to as the “GPS
start time”. In each file, there are in general two types of data – the strain
time series $h(t_{i})$ and some auxiliary data such as the data quality (DQ)
mask label associated with each strain measurement. The DQ mask is a 7-bit
binary number each of which indicates whether a certain type of check is
passed (value=1) or not (value=0). We convert this binary number into a
decimal digit. For example, there is no data at time $t_{i}$ if
DQ$(t_{i})=(0000000)_{10}=0$. On the contrary, data is present if this value
is nonzero.

Following the LIGO stochastic gravitational wave analysis Abbott  _et al._
(2017, 2019b), we first downsample the 16 kHz strain data to 4 kHz. We then
select out the timestamps at which both detectors are taking data properly,
i.e., the times $t_{i}$ at which both $\text{DQ}_{H}(t_{i})\neq 0$ and
$\text{DQ}_{L}(t_{i})\neq 0$. After doing that we get, for each data file, a
list of GPS times at which data in both detectors is available. We then
combine the lists of GPS times within a file and across neighbouring files
into continuous segments. The segments whose duration $\geq 600$ s are further
picked out to perform a stationarity cut following Ref. Abbott  _et al._
(2009, 2017). When we perform the stationarity cut, we also notch out
frequencies at which the data exhibits narrowband coherent lines that are
known to be instrumental or environmental artifacts Abbott  _et al._ (2016b);
Thrane  _et al._ (2013, 2014); Covas  _et al._ (2018). The list of the notched
frequency bands we used can be found on the public data release page:
https://dcc.ligo.org/LIGO-T1900058/public.

###  A.2 Analysis

After selecting a clean list of strain data and GPS time which satisfies the
requirements on data quality and stationarity, we then use the cross-
correlation method to estimate the SGWB signal. The spectrum of the GW
background is estimated with the cross-correlation statistic $\hat{C}(f)$
Abbott  _et al._ (2019b), defined as

$\displaystyle\hat{C}(f)$ $\displaystyle\equiv$
$\displaystyle\frac{2}{T}\frac{\textrm{Re}[\tilde{s}^{*}_{1}(f)\tilde{s}_{2}(f)]}{\gamma_{T}(f)S_{0}(f)},$
(40) $\displaystyle\langle\hat{C}(f)\rangle$ $\displaystyle=$
$\displaystyle\Omega_{\rm GW}(f),$ (41)

where $S_{0}=3H_{0}^{2}/(10\pi^{2}f^{3})$, $H_{0}$ is the Hubble parameter,
$\tilde{s}_{1,2}(f)$ are the Fourier transfroms of the strain data of both
detectors, $T=192$ s is the segment duration of the Fourier transforms, and
$\langle...\rangle$ indicates the average over all such 192s-segments. In the
limit that the GW signal is negligible comparing to the instrumental noise,
the variance of $\hat{C}(f)$ is given by

$\displaystyle\sigma^{2}(f)=\frac{1}{2T\Delta
f}\frac{\bar{P}_{1}(f)\bar{P}_{2}(f)}{\gamma_{T}^{2}(f)S_{0}^{2}(f)},$ (42)

where $\Delta f=1/32$ Hz, and $\bar{P}_{1,2}(f)$ are the one-sided power
spectrum of each detector, which are obtained as an average over two
neighbouring segments 111This means we need at least 3 continuous
192s-segments to obtain a power spectrum. Those 192s-segments without a
neighbor on both sides do not have an associated power spectrum.,

$\bar{P}_{i,I}\equiv\frac{P_{i,I-1}+P_{i,I+1}}{2},$ (43)

in which the subscript $i$ and $I$ labels the detector and the 192s-segments,
respectively. For each 192s-segment, we use the broadband estimator for any
spectral shape of the gravitational wave background,

$\displaystyle\hat{\Omega}_{\rm ref}$ $\displaystyle\equiv$
$\displaystyle\frac{\sum_{k}\omega(f_{k})^{-2}\hat{C}(f_{k})\sigma^{-2}(f_{k})}{\sum_{k}\omega(f_{k})^{-1}\sigma^{-2}(f_{k})}$
(44) $\displaystyle\langle\hat{\Omega}_{\rm ref}\rangle$ $\displaystyle=$
$\displaystyle\Omega_{\rm GW}(f_{\rm ref})$ (45)

where the weight function $\omega(f)\equiv\Omega_{\rm GW}(f_{\rm
ref})/\Omega_{\rm GW}(f)$, $f_{k}$ are discrete frequencies between 20 and
1726 Hz with the interval of $1/32$ Hz. The uncertainty of the optimal
estimator is,

$\sigma^{-2}_{\Omega}=\sum_{k}\omega(f_{k})^{-1}\sigma^{-2}(f_{k}).\\\ $ (46)

After calculating $\hat{\Omega}_{\rm ref}$ and $\sigma^{2}_{\Omega}$ for all
192s-segments, the ensemble average over all the segments is obtained from

$\displaystyle\sigma^{2}_{\Omega,{\rm tot}}$ $\displaystyle=$
$\displaystyle\frac{1}{\sum_{I}\sigma_{\Omega,I}^{-2}}\,,$ (47)
$\displaystyle\mu$ $\displaystyle=$
$\displaystyle\frac{\sum_{I}{\hat{\Omega}_{{\rm
ref},I}}{\sigma_{\Omega,I}^{-2}}}{\sum_{I}\sigma_{\Omega,I}^{-2}}\,.$ (48)

The signal-to-noise ratio (SNR) can be calculated by
$\mathrm{SNR}=\mu/\sigma$, where $\sigma\equiv\sqrt{\sigma_{\Omega,{\rm
tot}}^{2}}$. In the absence of detection signal, we set the 95% confidence
level upper limit by

$\Omega_{\rm ref}^{\rm up.~{}lim.}=\Omega_{\rm GW}^{\rm up.~{}lim}(f_{\rm
ref})=2\times\sigma_{\Omega,{\rm tot}}(f_{\rm ref})\,.$ (49)

Note that this upper limit depends on the choice of the reference frequency
$f_{\rm ref}$. In the calculation of Fourier transforms $\tilde{s_{i}}(f)$ and
power spectral density $P_{i}(f)$, we use the 50% overlapping Hann windows to
avoid spectral leakage Abbott  _et al._ (2017, 2019b, 2021).

## References

  * Abbott  _et al._ (2016a) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016a), arXiv:1602.03837 [gr-qc] .
  * Aasi  _et al._ (2015) J. Aasi, B. Abbott, R. Abbott, T. Abbott, M. Abernathy, K. Ackley, C. Adams, T. Adams, P. Addesso, R. Adhikari, _et al._ , Classical and quantum gravity 32, 074001 (2015).
  * Acernese  _et al._ (2014) F. a. Acernese, M. Agathos, K. Agatsuma, D. Aisa, N. Allemandou, A. Allocca, J. Amarni, P. Astone, G. Balestri, G. Ballardin, _et al._ , Classical and Quantum Gravity 32, 024001 (2014).
  * Abbott  _et al._ (2019a) B. Abbott, R. Abbott, T. Abbott, S. Abraham, F. Acernese, K. Ackley, C. Adams, R. Adhikari, V. Adya, C. Affeldt, _et al._ , Physical Review X 9, 031040 (2019a).
  * Abbott  _et al._ (2020) R. Abbott, T. Abbott, S. Abraham, F. Acernese, K. Ackley, A. Adams, C. Adams, R. Adhikari, V. Adya, C. Affeldt, _et al._ , arXiv preprint arXiv:2010.14527 (2020).
  * Allen and Romano (1999) B. Allen and J. D. Romano, Phys. Rev. D 59, 102001 (1999), arXiv:gr-qc/9710117 .
  * Maggiore (2000) M. Maggiore, Physics Reports 331, 283 (2000).
  * Regimbau (2011) T. Regimbau, Research in Astronomy and Astrophysics 11, 369 (2011).
  * Caprini and Figueroa (2018) C. Caprini and D. G. Figueroa, Class. Quant. Grav. 35, 163001 (2018), arXiv:1801.04268 [astro-ph.CO] .
  * Croon  _et al._ (2018) D. Croon, V. Sanz, and G. White, JHEP 08, 203 (2018), arXiv:1806.02332 [hep-ph] .
  * Romero  _et al._ (2021) A. Romero, K. Martinovic, T. A. Callister, H.-K. Guo, M. Martínez, M. Sakellariadou, F.-W. Yang, and Y. Zhao, (2021), arXiv:2102.01714 [hep-ph] .
  * Caprini  _et al._ (2008) C. Caprini, R. Durrer, and G. Servant, Phys. Rev. D 77, 124015 (2008), arXiv:0711.2593 [astro-ph] .
  * Huber and Konstandin (2008) S. J. Huber and T. Konstandin, JCAP 09, 022 (2008), arXiv:0806.1828 [hep-ph] .
  * Caprini  _et al._ (2009) C. Caprini, R. Durrer, and G. Servant, JCAP 12, 024 (2009), arXiv:0909.0622 [astro-ph.CO] .
  * Espinosa  _et al._ (2010) J. R. Espinosa, T. Konstandin, J. M. No, and G. Servant, JCAP 06, 028 (2010), arXiv:1004.4187 [hep-ph] .
  * Hindmarsh  _et al._ (2014) M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, Phys. Rev. Lett. 112, 041301 (2014), arXiv:1304.2433 [hep-ph] .
  * Hindmarsh (2018) M. Hindmarsh, Phys. Rev. Lett. 120, 071301 (2018), arXiv:1608.04735 [astro-ph.CO] .
  * Bodeker and Moore (2017) D. Bodeker and G. D. Moore, JCAP 05, 025 (2017), arXiv:1703.08215 [hep-ph] .
  * Weir (2018) D. J. Weir, Phil. Trans. Roy. Soc. Lond. A 376, 20170126 (2018), arXiv:1705.01783 [hep-ph] .
  * Abbott  _et al._ (2017) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 118, 121101 (2017), [Erratum: Phys.Rev.Lett. 119, 029901 (2017)], arXiv:1612.02029 [gr-qc] .
  * Abbott  _et al._ (2019b) B. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. D 100, 061101 (2019b), arXiv:1903.02886 [gr-qc] .
  * Tanin and Tenkanen (2020) E. H. Tanin and T. Tenkanen, (2020), arXiv:2004.10702 [astro-ph.CO] .
  * Croon and White (2018) D. Croon and G. White, JHEP 05, 210 (2018), arXiv:1803.05438 [hep-ph] .
  * Croon  _et al._ (2019a) D. Croon, T. E. Gonzalo, and G. White, JHEP 02, 083 (2019a), arXiv:1812.02747 [hep-ph] .
  * Croon  _et al._ (2019b) D. Croon, T. E. Gonzalo, L. Graf, N. Košnik, and G. White, Front. in Phys. 7, 76 (2019b), arXiv:1903.04977 [hep-ph] .
  * Huang  _et al._ (2020) W.-C. Huang, M. Reichert, F. Sannino, and Z.-W. Wang, (2020), arXiv:2012.11614 [hep-ph] .
  * Guo  _et al._ (2021) H.-K. Guo, K. Sinha, D. Vagie, and G. White, JCAP 01, 001 (2021), arXiv:2007.08537 [hep-ph] .
  * Croon  _et al._ (2020) D. Croon, O. Gould, P. Schicho, T. V. I. Tenkanen, and G. White, (2020), arXiv:2009.10080 [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] .
  * Ellis  _et al._ (2020) J. Ellis, M. Lewicki, and V. Vaskonen, JCAP 11, 020 (2020), arXiv:2007.15586 [astro-ph.CO] .
  * Lewicki and Vaskonen (2020a) M. Lewicki and V. Vaskonen, Eur. Phys. J. C 80, 1003 (2020a), arXiv:2007.04967 [astro-ph.CO] .
  * Lewicki and Vaskonen (2020b) M. Lewicki and V. Vaskonen, (2020b), arXiv:2012.07826 [astro-ph.CO] .
  * Azatov and Vanvlasselaer (2021) A. Azatov and M. Vanvlasselaer, JCAP 01, 058 (2021), arXiv:2010.02590 [hep-ph] .
  * Abbott  _et al._ (2021) R. Abbott _et al._ (LIGO Scientific, Virgo, KAGRA), (2021), arXiv:2101.12130 [gr-qc] .
  * Dienes  _et al._ (2020) K. R. Dienes, F. Huang, J. Kost, S. Su, and B. Thomas, Phys. Rev. D 101, 123511 (2020), arXiv:2001.02193 [astro-ph.CO] .
  * Dienes  _et al._ (2021) K. R. Dienes, F. Huang, J. Kost, K. Manogue, and B. Thomas, (2021), arXiv:2101.10337 [astro-ph.CO] .
  * Abbott  _et al._ (2009) B. P. Abbott _et al._ (LIGO Scientific, VIRGO), Nature 460, 990 (2009), arXiv:0910.5772 [astro-ph.CO] .
  * Abbott  _et al._ (2016b) B. P. Abbott, R. Abbott, T. Abbott, M. Abernathy, F. Acernese, K. Ackley, M. Adamo, C. Adams, T. Adams, P. Addesso, _et al._ , Classical and Quantum Gravity 33, 134001 (2016b).
  * Thrane  _et al._ (2013) E. Thrane, N. Christensen, and R. M. Schofield, Physical Review D 87, 123009 (2013).
  * Thrane  _et al._ (2014) E. Thrane, N. Christensen, R. M. Schofield, and A. Effler, Physical Review D 90, 023013 (2014).
  * Covas  _et al._ (2018) P. Covas, A. Effler, E. Goetz, P. Meyers, A. Neunzert, M. Oliver, B. Pearlstone, V. Roma, R. Schofield, V. Adya, _et al._ , Physical Review D 97, 082002 (2018).


# InProC: Industry and Product/Service Code Classification

Simerjot Kaur  JPMorgan AI ResearchPalo AltoCaliforniaUSA
<EMAIL_ADDRESS>,  Andrea Stefanucci  JPMorgan AI ResearchNew
YorkNYUSA<EMAIL_ADDRESS>and  Sameena Shah  JPMorgan AI
ResearchNew YorkNYUSA<EMAIL_ADDRESS>
(2022)

###### Abstract.

Determining industry and product/service codes for a company is an important
real-world task and is typically very expensive as it involves manual curation
of data about the companies. Building an AI agent that can predict these codes
automatically can significantly help reduce costs, and eliminate human biases
and errors. However, unavailability of labeled datasets as well as the need
for high precision results within the financial domain makes this a
challenging problem. In this work, we propose a hierarchical multi-class
industry code classifier with a targeted multi-label product/service code
classifier leveraging advances in unsupervised representation learning
techniques. We demonstrate how a high quality industry and product/service
code classification system can be built using extremely limited labeled
dataset. We evaluate our approach on a dataset of more than 20,000 companies
and achieved a classification accuracy of more than 92%. Additionally, we also
compared our approach with a dataset of 350 manually labeled product/service
codes provided by Subject Matter Experts (SMEs) and obtained an accuracy of
more than 96% resulting in real-life adoption within the financial domain.

small datasets, representation learning, classification

††copyright: acmcopyright††journalyear: 2022††conference: ICAIF’22; Small
Data, Big Opportunities: Making the Most of AI; 2022††booktitle: ICAIF, Small
Data, Big Opportunities: Making the Most of AI††price:

##  1\. Introduction

Determining the industry sectors to which a company belongs to, and the class
of products and services it is involved in, is an important real-world task
with many applications in economic, financial, advertising, marketing and
other business firms. For instance, such categorization is commonly used for
prospective client-lead generation by sales and marketing teams, for trading
strategies like pairs-trading which rely upon highly correlated stocks within
the same industry and dealing with similar products, recommending companies to
investors and vice-versa based on historical industry and product patterns.
Moreover, the knowledge about company’s industry, products and services helps
in performing industry analysis about the company like general industry
economics, industry participants, competition, distribution and buying
patterns, etc.

However, correctly classifying which industries each company belongs to and
categorizing its products/services is typically very expensive as it involves
manual curation of data about the companies, analyzing all the sectors they
might be involved in by going through details of the various products they
have built, etc. This quickly becomes a tedious task involving significant
amounts of manual labour. Further, it brings about human errors and the
results are often dependent on the individual experience and expertise of
humans involved in this process. Additionally, with the interdisciplinary
nature of the new products and services that companies are creating, a
majority of the companies belong to more than one industry, for instance, a
company building an educational platform on a mobile app belongs to three
industries: education, software and mobile applications, and correctly
classifying each company to multiple classes manually is extremely tedious and
prone to errors. Furthermore, the lack of a common taxonomy of industry
sections and classes of products and services that is useful and agreed upon
across different institutions and firms further exacerbates the problem, and
necessities such internal development activities at each individual firm. It
may be noted that we are not interested in classifying only the available
public traded companies, but rather all companies, i.e. public or private,
small or large, etc.

Figure 1. Problem Statement

To overcome the above limitations and the need for a human expert with years
of experience, we propose to build a system for classifying the company’s
industries and categorizing its products and services in an automated manner
using a hybrid deep-learning based hierarchical classifier. However, most
current classification techniques either involve building a supervised multi-
class multi-label classifier that require thousands of labeled dataset or
unsupervised hierarchical clustering approach which are not easily scalable
and have low precision. One of the standard taxonomy datasets for industry
code classification is the North American Industry Classification System
(NAICS) (Association, 2017) codes but due to its lack of adaptation to
changing economy, NAICS codes have not gained real-life adoption by most
financial organizations. Hence, in case of financial applications,
unavailability of labeled dataset and direct financial consequences makes it a
hard problem to solve and poses a serious challenge to adoption within the
financial domain. In this work, we propose a novel algorithm that leverages
recent advances in representation learning to build a hierarchical multi-class
classifier along with a targeted multi-label classifier which predicts which
industries each company belongs to and which products and services the company
builds within each industry by only using company’s description as input.

Formally, Fig 1 above summarizes our problem. Given a limited sample-set of
companies A, B…N, together with company descriptions, can we correctly predict
which custom-defined industry codes each company belongs to and what products
and services within each industry the company is involved in?

Our proposed approach consists of first constructing a labeled dataset for
training a multi-class industry classifier by leveraging human subject matter
experts (SMEs) and unsupervised techniques. For this work, a pre-defined set
of custom named industries and products and services within each industry were
provided by the SMEs 111To protect proprietary data, the names of pre-defined
set of custom industries and products/services have not been disclosed.. The
SMEs also provided brief one-to-two line description for each of the industry
as well as the product and services taxonomy. We then build a high dimensional
vector representation for each company using its description as available
input data and predict top 3 industries in which the company might belong to.
Finally, we construct a targeted multi-class products and services classifier
by first building a high dimensional vector representation for each company as
well as predefined products and services taxonomy within the predicted
industries using their descriptions as inputs and generate similarity scores
between them which are then used to identify top 2 products and services that
the company might be involved in.

The rest of the paper is outlined as follows. Section 2 describes related work
and highlights the unique challenges that come about in building an industry
and product/service code classifier, section 3 describes the details of how we
constructed the dataset and section 4 provides details of our proposed
hierarchical architecture and implementation in more detail. In section 5 we
review the results of our approach and also demonstrate how the human targeted
SME verification helped in the adoption of our model in the real-world of
financial domain. Finally, section 6 concludes our work and lays out
foundation for future work in this area.

##  2\. Related Work

There has been numerous works that focus on solving hierarchical text
classification which try to leverage recent advances in deep learning and
solve the task using supervised techniques such as (Yao et al., 2019), (Song
et al., 2020), (Luo, 2021), and (Lavanya and Sasikala, 2021) has done a
comprehensive survey on these techniques. While many of these techniques have
achieved remarkable results, they require huge amounts of labeled hierarchical
data which is very hard to obtain in the financial domain. As mentioned in
Section 1, one of the available taxonomy datasets for industry code
classification is the NAICS (Association, 2017) codes, a six-digit code,
containing 1057 codes which comprises of 20 sectors, 99 sub-sectors, and 311
industry groups. However, this dataset is very focused on North America
economy and has not been adaptive enough to capture the changing economy and
hence has not been adopted by most financial organizations. This makes the
problem very hard to solve as every financial organization defined their own
custom industry and products and services codes.

Some groups have tried to solve the text classification problem by using
unsupervised techniques like hierarchical clustering and topic modeling such
as (Aljedani et al., 2021), (Zhang et al., 2021), (Rafea et al., 2021).
However, these models are hard to get adopted in real-world financial domain
as they do not perform well on large scale data and the low precision results
can have direct financial consequences as explained in Section 1.

Additionally, there also has been some work in the literature on automated
industry classification (Lim and Lee, 2005), (Phillips and Ormsby, 2016),
(Husmann et al., 2020), almost all these approaches rely on self- reporting,
manual entry, and possibly naive algorithms. (Wood et al., 2017) and (Oehlert
et al., 2022) explored the possibility of using NAICS code for building a
supervised learning algorithm and leverage advancements in deep learning.
(Papagiannidis et al., 2018) tried to use clustering algorithms to capture
similarities across companies and identify industry codes. These methodologies
are not scalable and are difficult to be adopted and deployed in financial
domain since the companies these days belong to multiple industries and the
products and services are interdisciplinary in nature. Moreover, there has
been no literature on hierarchical classification of the companies into
identifying what products and services the companies are involved in. In this
work, we have tried to leverage advances in deep learning and unsupervised
learning to build a hierarchical multi-class multi-label classifier which uses
the company’s descriptions as input and predicts top 3 industries that the
company might belong to and top 2 product/services within each industry that
the company might be involved in.

##  3\. Data

This section describes the details for our dataset. Our model has been
developed and tested using licensed data from Pitchbook (Pitchbook Dataset,
2021). The dataset contains details on companies such as a brief description
on what the company does as well as pitchbook-defined industry verticals.
Moreover, we were provided with $\sim$40 custom defined industry codes and
$\sim$400 products and services codes by the SMEs into which the companies had
to be classified in. The following sections detail how we constructed the
dataset to build an industry code classifier as well as the products and
services classifier.

###  3.1. Dataset Construction for Industry Code Classification

To classify the companies with high accuracy and precision into which
industries they belong to, we first construct a labeled dataset with input as
company descriptions and output as the industry codes. In order to construct
this dataset, we leveraged the pitchbook industry codes. The pitchbook
industry codes are defined in three levels: (a)Primary Industry Sector: A
broad category that contains industry groups and codes, (b)Primary Industry
Group: A sub-category that provides more specific classification, (c)Primary
Industry Code: The primary industry the company operates in.

In order to obtain the mapping between pitchbook descriptions and SME’s custom
defined industry codes, we leveraged our SMEs to create a mapping between the
pitchbook’s three level codes and SME’s industry codes. Through this mapping
we were able to generate a labeled dataset containing pitchbook descriptions
as inputs and SME’s industry codes as output. However, through this mapping we
were only able to generate the labeled dataset for 75% of the industry codes.

Figure 2. Construction of Labeled Data for the remaining 25% of the industry
codes for Industry Code Classification

For the remaining 25%, we used an unsupervised approach to obtain the labeled
dataset. Fig 2 outlines the methodology used to construct the dataset for the
remaining industry codes. Since the industry codes in itself do not contain
much contextual meaning, the SME’s provided us with brief one-to-two line
descriptions of all the custom defined industry codes. In order to enable us
to easily perform comparisons and estimate how close/far the companies and the
industry codes are, we then used a pre-trained large language model to encode
the company descriptions and industry code descriptions into a high
dimensional distributed vector representation. The key guiding principle
behind the high dimensional vector representation is that similar
sentences/words, or phrases within a similar contexts should map close to each
other in a high-dimensional space, and sentences/words or phrases that are
very different should be far away from each other. We then estimate the
similarity scores between the two vector representations using cosine
similarity, $cos(comp,ind\\_code)$. The industry code which has the highest
similarity score and is above a given threshold, $thresh$ (hyperparameter), is
then assigned as the label to the corresponding company description. Please
note that these comparisons are performed only for the remaining 25% of the
industry codes.

(1)
$cos(comp,ind\\_code)=\frac{comp\\_vec.ind\\_vec}{||comp\\_vec||.||ind\\_vec||}$

###  3.2. Dataset Construction for Products and Services Code Classification

Since we use an unsupervised approach for this classification, we only used
the company descriptions from pitchbook dataset. As the product and service
codes in itself do not contain much contextual meaning, the SME’s provided us
with one-to-two line descriptions of all the custom defined product and
service codes.

##  4\. Proposed Approach

This section describes our algorithm used to classify the companies into the
industry codes to which it belongs, in Subsection 4.1, and subsequently how it
predicts the product and service codes in which the companies are involved in,
in Subsection 4.2.

Figure 3. Industry Codes Classification Methodology

###  4.1. Industry Code Classifier

Figure 3 describes our proposed approach to solve the industry code
classification problem i.e. ’What are the industry codes to which Company $i$
belongs to?’ . Our proposed approach consists of two sequential steps. In the
first step, we pass the company descriptions through a pre-trained large
language model, specifically Roberta-base model (Liu et al., 2019), and obtain
high dimensional distributed vector representations for each company. In the
second step, we use these representations as input features to train a
multilayer perceptron with a classification layer using the labeled dataset
constructed in Section 3.1. This trained classification model is then used to
predict industry codes for new companies. Moreover, since the Roberta-base
model is pre-trained over generic Wikipedia data, this makes our trained
classifier generalizable to predict industry codes for company descriptions
obtained from any data source.

###  4.2. Products and Services Code Classifier

Figure 4. Product/Service Codes Classification Methodology

Figure 4 describes our proposed approach to solve the product and service code
classification problem i.e. ’What are the product and services codes within a
particular industry that the Company $i$ is involved in?’ . This is the second
targeted classification layer wherein we try to predict which product and
service codes the company is involved in within each predicted industry code
($\sim$8-15 product/services codes within each industry) in an unsupervised
manner, obtained from Subsection 4.1. In this approach, we obtained the
products and services code descriptions within each industry code from the
SMEs, as described in Section 3.2. In order to enable us to easily perform
comparisons and estimate how close/far the companies and the product and
service codes are, we then used a pre-trained large language model,
specifically Roberta-base model, to encode the company descriptions and
product and service code descriptions into a high dimensional distributed
vector representation. We then estimate the similarity scores between the two
vector representations using cosine similarity, Equation 1. The product and
service codes which have the top 2 similarity scores are then assigned as the
predicted product and service codes to the corresponding company.

The Algorithm 1 describes our overall industry and product/service codes
prediction algorithm in a step-by-step manner.

Algorithm 1 Industry & Product/Service Codes Prediction Algo

1: Input 1: Companies, say $comp$ where $comp\in[1,2,3,\ldots,m]$, its
descriptions, say $desc$ where $desc\in[1,2,3,\ldots,m]$ and its labels, say
$l$ where $l\in[1,2,\ldots,m]$, where $m$ is the number of companies in the
dataset

2: Input 2: Products and services taxonomy, say $ps\\_code$ where
$ps\\_code\in[A,B,\ldots,N]$ and its descriptions, say $ps\\_code\\_desc$
where $ps\\_code\\_desc\in[A,B,\ldots,N]$, where $N$ is number of products and
services codes belonging to each industry code ($N$ varies for each industry
code)

3: (1) Industry Codes Classifier:

4: Split the dataset into training set, say
$comp\\_tr,desc\\_tr,l\\_tr\in[1,2,3,\ldots,k]$, and test set, say
$comp\\_ts,desc\\_ts,l\\_ts\in[1,2,3,\ldots,k]$

5: for $comp\\_tr,desc\\_tr=[1,2,\ldots,k]$ do

6: Obtain $tr\\_v=$ vector rep. for $desc\\_tr$

7: end for

8: Train classification model, say $ind\\_code\\_mo$, using $tr\\_v$ as inputs
and $l\\_tr$ as actual labels

9: for $comp\\_ts,desc\\_ts=[1,2,\ldots,k]$ do

10: Obtain $ts\\_v=$ vector rep. for $desc\\_ts$

11: Use $ind\\_code\\_mo$ and predict top 3 industry codes for $comp\\_ts$,
say $comp\\_InC\in[X,Y,Z]$

12: end for

13: (2) Products and Services Codes Classifier:

14: for $comp\\_ts,desc\\_ts=[1,2,\ldots,k]$ do

15: for $comp\\_InC=[X,Y,Z]$ do

16: for $ps\\_code=[A,B,\ldots,N]$ do

17: Obtain $ps\\_v=$ vector rep. for $ps\\_code$ description

18: Obtain $ts\\_v=$ vector rep. for $desc\\_ts$

19: Calculate cosine
similarity,$sim(comp,ps\\_code)=\frac{ts\\_v.ps\\_v}{||ts\\_v||.||ps\\_v||}$

20: end for

21: Obtain top 2 $prod\\_codes$ for which the $sim(comp,ps\\_code)$ are the
highest

22: end for

23: end for

##  5\. Experiments and Results

Figure 5. Confusion Matrix for a sample of 10 industry codes

In order to first evaluate the performance of our industry code classifier, we
first split the labeled dataset constructed in Section 3.1 into 80% training
and 20% test set. Moreover, as discussed in Sections 1 and 2, since the
companies are interdisciplinary in nature we generated the confusion matrix
and observed that majority of the companies cover the span of at least 3
industry codes. Fig 5 shows the confusion matrix for a sample of 10 industry
codes222Due to proprietary reasons, the industry codes have been masked. with
red circles depicting the span size for $Industry$ $Code$ $5$ . Finally, Table
1 shows top-3-accuracy score of the industry code classifier on this labeled
test set.

Further, in order to test the accuracy of product and services classifier as
well as for real-life adoption, we were provided a small set of 350 companies
which were manually labeled by Subject Matter Experts (SMEs) containing the
ground truth for both industry and product and service. Table 1 shows
top-3-accuracy score of industry classifier and top-2-accuracy score of
product and services classifier on this SME-labeled dataset.

Industry Code Classifier Products and Services Classifier # of Samples
Constructed Test Set 92.38% - 20,000 SME-Labeled Dataset 94.87% 96.85% 350

Table 1. Results of Industry Code Classifier and Product and Service Code
Classifier

##  6\. Conclusion and Future Work

Our work proposes a novel hierarchical classification algorithm that leverages
recent advances in representation learning to predict high precision custom-
defined industry codes as well as products and services codes for a particular
company by only using company’s description as input. We also highlight how a
high quality labeled industry code dataset can be constructed through mapping
and unsupervised techniques, and demonstrate how a targeted unsupervised
multi-class products and services classifier can yield high precision
predictions resulting in real-life adoption and deployment.

The work opens numerous avenues to further build upon. For instance, we could
further enhance the algorithm to automatically generate custom industry codes
and product and services codes based on company descriptions and then using
public information extraction techniques to extract code definitions hence
eliminating the need for predefining these codes by a subject matter expert.

## References

  * (1)
  * Aljedani et al. (2021) Nawal Aljedani, Reem Alotaibi, and Mounira Taileb. 2021\.  Hmatc: Hierarchical multi-label arabic text classification model using machine learning.  _Egyptian Informatics Journal_ 22, 3 (2021), 225–237.
  * Association (2017) NAICS Association. 2017\.  North American Industry Classification System (NAICS) code.  _http://www.naics.com/frequently-asked-questions/NAICSfaq_ (2017).
  * Husmann et al. (2020) Sven Husmann, Antoniya Shivarova, and Rick Steinert. 2020\.  Company classification using machine learning.  https://doi.org/10.48550/ARXIV.2004.01496
  * Lavanya and Sasikala (2021) PM Lavanya and E Sasikala. 2021.  Deep learning techniques on text classification using Natural language processing (NLP) in social healthcare network: A comprehensive survey. In _2021 3rd International Conference on Signal Processing and Communication (ICPSC)_. IEEE, 603–609.
  * Lim and Lee (2005) Heui Seok Lim and Seong Hoon Lee. 2005.  An Application of Information Retrieval Technique to Automated Code Classification. In _Knowledge-Based Intelligent Information and Engineering Systems_ , Rajiv Khosla, Robert J. Howlett, and Lakhmi C. Jain (Eds.). Springer Berlin Heidelberg, Berlin, Heidelberg, 90–96.
  * Liu et al. (2019) Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov. 2019.  RoBERTa: A Robustly Optimized BERT Pretraining Approach.  https://doi.org/10.48550/ARXIV.1907.11692
  * Luo (2021) Xiaoyu Luo. 2021\.  Efficient english text classification using selected machine learning techniques.  _Alexandria Engineering Journal_ 60, 3 (2021), 3401–3409.
  * Oehlert et al. (2022) Christine Oehlert, Evan Schulz, and Anne Parker. 2022\.  NAICS Code Prediction Using Supervised Methods.  _Statistics and Public Policy_ 9, 1 (2022), 58–66.
  * Papagiannidis et al. (2018) Savvas Papagiannidis, Eric WK See-To, Dimitris G Assimakopoulos, and Yang Yang. 2018\.  Identifying industrial clusters with a novel big-data methodology: Are SIC codes (not) fit for purpose in the Internet age?  _Computers & Operations Research_ 98 (2018), 355–366.
  * Phillips and Ormsby (2016) Ryan L. Phillips and Rita Ormsby. 2016.  Industry classification schemes: An analysis and review.  _Journal of Business & Finance Librarianship_ 21, 1 (2016), 1–25.  https://doi.org/10.1080/08963568.2015.1110229 arXiv:https://doi.org/10.1080/08963568.2015.1110229
  * Pitchbook Dataset (2021) Pitchbook Dataset 2021.  Pitchbook Dataset: All Company, Investor, Deals Dataset. Retrieved April 15, 2021 from PitchBook database.  https://get.pitchbook.com/pitchbook-data/.
  * Rafea et al. (2021) Layth Rafea, Abdulrahman Ahmed, and Wisam D Abdullah. 2021\.  Classification of a COVID-19 dataset by using labels created from clustering algorithms.  _Indonesian Journal of Electrical Engineering and Computer Science_ 21, 1 (2021), 164–173.
  * Song et al. (2020) Min Song, Keun Young Kang, Tatsawan Timakum, and Xinyuan Zhang. 2020.  Examining influential factors for acknowledgements classification using supervised learning.  _PloS one_ 15, 2 (2020), e0228928.
  * Wood et al. (2017) Sam Wood, Rohit Muthyala, Yi Jin, Yixing Qin, Nilaj Rukadikar, Amit Rai, and Hua Gao. 2017.  Automated industry classification with deep learning. In _2017 IEEE International Conference on Big Data (Big Data)_. IEEE, 122–129.
  * Yao et al. (2019) Liang Yao, Chengsheng Mao, and Yuan Luo. 2019.  Clinical text classification with rule-based features and knowledge-guided convolutional neural networks.  _BMC medical informatics and decision making_ 19, 3 (2019), 31–39.
  * Zhang et al. (2021) Yu Zhang, Zhihong Shen, Yuxiao Dong, Kuansan Wang, and Jiawei Han. 2021.  MATCH: Metadata-aware text classification in a large hierarchy. In _Proceedings of the Web Conference 2021_. 3246–3257.