# A Multi-Sensor Interface to Improve the Learning Experience
in Arc Welding Training Tasks

Hoi-Yin Lee, , Peng Zhou, , Anqing Duan, ,
Jiangliu Wang, Victor Wu and David Navarro-Alarcon This work has been
submitted to the IEEE for possible publication. Copyright may be transferred
without notice, after which this version may no longer be accessible
https://journals.ieeeauthorcenter.ieee.org/become-an-ieee-journal-
author/publishing-ethics/guidelines-and-policies/post-publication-
policies/This work is supported by the Research Grants Council (RGC) under
grant 15212721. Corresponding author: D. Navarro-Alarcon.H.-Y. Lee, P. Zhou,
A. Duan, V. Wu and D. Navarro-Alarcon are with The Hong Kong Polytechnic
University, Department of Mechanical Engineering, Hung Hom, KLN, Hong Kong.
(contact e-mail: dna@ieee.org)J. Wang is with the CUHK T Stone Robotics
Institute, The Chinese University of Hong Kong, Shatin, NT, Hong Kong.

###### Abstract

This paper presents the development of a multi-sensor user interface to
facilitate the instruction of arc welding tasks. Traditional methods to
acquire hand-eye coordination skills are typically conducted through one-to-
one instruction where trainees must wear protective helmets and conduct
several tests. This approach is inefficient as the harmful light emitted from
the electric arc impedes the close monitoring of the process; Practitioners
can only observe a small bright spot. To tackle these problems, recent
training approaches have leveraged virtual reality to safely simulate the
process and visualize the geometry of the workpieces. However, the synthetic
nature of these types of simulation platforms reduces their effectiveness as
they fail to comprise actual welding interactions with the environment, which
hinders the trainees’ learning process. To provide users with a real welding
experience, we have developed a new multi-sensor extended reality platform for
arc welding training. Our system is composed of: (1) An HDR camera, monitoring
the real welding spot in real-time; (2) A depth sensor, capturing the 3D
geometry of the scene; and (3) A head-mounted VR display, visualizing the
process safely. Our innovative platform provides users with a “bot trainer”,
virtual cues of the seam geometry, automatic spot tracking, and performance
scores. To validate the platform’s feasibility, we conduct extensive
experiments with several welding training tasks. We show that compared with
the traditional training practice and recent virtual reality approaches, our
automated multi-sensor method achieves better performances in terms of
accuracy, learning curve, and effectiveness.

###### Index Terms:

Multisensory displays, human-computer interface, virtual reality, automation,
manufacturing.

##  I Introduction

Arc welding is one of the most common material fusing methods in modern
manufacturing [1]. In its most basic form, it uses a controllable electric
current to melt a joining metal that, once cooled, binds two metallic parts
together [2]. Due to its strong and enduring joining properties, welding is
used across numerous economically important fields, such as automotive and
aerospace industries, shipbuilding, steel construction, oil and gas pipelines,
to name a few instances [3]. Despite its widespread use, teaching and learning
a proper welding technique has historically presented many challenges to both
instructors and trainees [4].

Traditional approaches to acquire welding skills are typically conducted as
follows: The instructor first explains the working principle of the process
and demonstrates it with sample welds [5]; Trainees proceed to conduct hands-
on welding tasks under the guidance from the instructor, see Fig. 1. The main
limitation of this approach comes from the inability to clearly observe the
workpiece geometry and the process through the protective helmet, which must
be worn at all times [6, 7, 8] (beginners generally struggle to obtain a clear
spatial notion of what is occurring on the other side of the helmet). The
near-dark experience makes it difficult to learn how the different
configurations affect the quality of a weld; It is hard for trainees to
understand a process they can barely observe. Furthermore, as there is a
limited number of instructors that generally participate in a session,
trainees cannot receive real-time support during their practice (advice is
received only after the task has been completed and the helmet removed). All
these factors complicate the instruction and skill acquisition of arc welding
tasks.

Figure 1: (a) Arc welding; (b) Illustration of traditional training
approaches, where the trainer needs to demonstrate the welding skill one by
one.

To address the limitations of traditional training approaches, researchers
have developed various didactic platforms based on virtual/mixed reality
(VR/MR) systems [9, 10, 11, 12, 13, 14]. These platforms are typically
designed to help beginners to get familiar with the procedure and its torch
movements through the use of virtual environments that simulate the _high-
energy_ welding process. As no actual interaction with the environment occurs,
these simulation-based systems can only provide a synthetic perceptual
experience to the users, which may hinder the psychological adaptation that
practitioners acquire by conducting real tasks in the field (e.g., fear
management, thermal and noise sensations, etc. [15]).

To address the above mentioned issues, in this paper we propose an automated
extended reality (XR) training assistant platform that provides the user with
a multimodal experience. The proposed system consists of a high dynamic range
(HDR) camera [16, 17] to monitor the welding spot in real-time, an RGB-D
sensor [18, 19] to capture the 3D geometry of the scene, and a VR headset to
safely visualize the process; Images from these sensors [20] are registered
into a common frame. The system uses 3D vision to detect the welding region
and seam, and to automatically generate the desired welding path [21]. These
visual cues from the virtual assistant provide the user with valuable extended
reality guidance on the weld motion in real-time. The proposed automated “bot
trainer” provides the user with a scoring system that quantifies the
performance of a task upon its completion. As visual information is displayed
on a screen in XR mode, the process can be simultaneously observed by the
trainee, instructor, and other participants, see Fig. 2.

The goal of this XR system is to help trainees to learn the proper torch
movements and gain confidence with real welding tasks while reducing the
required human instruction to a minimum. We conducted a series of experiments
to evaluate our proposed approach in terms of usability, performance, learning
curve, and teaching/learning effectiveness.

Compared with existing approaches, the developed automated training system has
the following original features:

  * •

It detects and overlays a virtual welding path over the workpiece for trainees
to follow.

  * •

It provides the user with instantaneous motion recommendations to improve the
task performance.

  * •

It quantifies and visualizes the performance of the conducted welding task
based on sensory feedback.

  * •

It enables to display the process in real-time to the user and other
participants in the training session.

The rest of this paper is organized as follows: Sec. II introduces the
architecture of the system. Sec. III presents the experiments. Sec. IV gives
the final conclusion.

##  II Methods

Figure 2: Conceptual representation of the proposed training system.

###  II-A System Overview

The data flow among the different components of the XR bot trainer is shown in
Fig. 3. It includes an RGB-D sensor, an HDR camera, a VR headset with a
controller, a computer, and a welding torch. Cameras are cross-calibrated and
placed at around 10-degrees from the surface normal and facing towards the
welding region. Multimodal visual feedback is sent to the computer for
processing, and then fed into the VR headset to provide the user with 3D
models of the workpiece and a video stream of the welding process.

Figure 3: Workflow between the different components in the system.

The workflow depicted in Fig. 3 describes the following steps: Before the
process begins, (1) the RGB-D sensor captures a depth image of the welding
area to produce a point cloud of the workpiece; (2) A depth-based localization
algorithm segments the area of interest and computes the welding seam/path
[22, 23]. During welding, (3) the HDR camera captures greyscale images of the
process and traces the center of the electric arc in real-time, where (4) all
visual information is registered into a 2D image and streamed to a webpage;
(5) The VR headset access the 2D live streaming video and the 3D model via Wi-
Fi.

The electrode of the welding torch generates a strong electric arc (whose
emitted light is harmful to the human eyes) when approaching an electrified
workpiece. The position of this bright spot is tracked by the system, and its
trace is visualized to the user to provide valuable visual cues. Upon
completion of the task, an evaluation metric is calculated to assess the
overall performance; The welding trajectory, average error, and score are
displayed on the interface to quantify the performance of the user.

###  II-B Seam Localization

A groove is a channel between the edges of two metal workpieces [24]. The
RGB-D sensor captures an initial color-depth image that covers the whole
welding area of the workpiece; This sensor data forms a point cloud
$\mathbf{P}$ of the scene. The proposed detection algorithm uses the
difference in edge intensity to automatically find the path of the groove,
where local neighborhoods in $\mathbf{P}$ are created to segment the groove’s
approximated location [22, 23]. In Fig. 4, the groove points
$G=\\{g_{1},\ldots,g_{n}\\}$ are conceptually presented in green color over
the channel. Computing the bulk of data from this 3D point cloud results in a
long processing time; Several points are involved in finding the seam, which
(for beginners) is typically a simple straight line. Therefore, the following
efficient method for locating the seam is implemented: After the groove has
been segmented in the image, the seam is calculated by projecting all points
in ${G}$ into a 2D image $\mathbf{I}$, as represented with red color in Fig.
5. The groove’s edges are then calculated from 2D coordinates.

Noisy data from sensor measurements may affect the precision of the seam
localization, thus, a kernel convolution noise filter is implemented on
$\mathbf{I}$ to remove noisy/redundant points. Canny edge detection [25] is
then applied to $\mathbf{I}$ to sharpen the edges of the groove. Inspired by
the method in [26] (which uses an edge map to extract line segments), we
compute a point-line connection edge map to improve the precision of our
results. Points are interconnected with line segments, as depicted in Fig. 5.
Euclidean distance $\lVert{a-b}\rVert$ is then calculated to determine the
length of lines with endpoints $\mathnormal{a}$ and $\mathnormal{b}$. As the
welding seam is expected to be longer than the sides of the workpiece, our
method removes all short lines.

Figure 4: (a) Two different types of workpieces: Fillet welding and Butt
welding workpieces; (b) The groove detection algorithm is applied to locate
the possible welding area. Groove points $G$ are displayed in green color; (c)
Welding seams are indicated in orange. Figure 5: (a) 3D coordinates of the
groove ${G}$ are projected onto a 2D image $\mathbf{I}$ in red; (b) The
computed lines around the seam are the orange and green long lines; Blue lines
represent the side edges of the workpiece, which are generally shorter; (c)
Lines in blue and green in (b) are removed due to its length and slope.

To extract the characteristic line $l$ representing the seam, we compute the
average slope of $n$ (long) line segments as:

$\overline{s}=\frac{1}{n}\sum_{i=1}^{n}s_{i}$ (1)

for $\mathnormal{s}_{i}$ as the slope of the $i$th line segment. Lines with a
slope outside a range $[\overline{s}-c,\overline{s}+c]$, for
$\mathnormal{c}>0$ as clearance scalar, are discarded. This enables us to
narrow down the possible seam region to the orange region depicted in Fig. 5.

A KD-Tree algorithm [27] is introduced to classify the endpoints
$a,b=\mathnormal{E(x,y)}$ of these lines into groups according to their image
quadrants. The two quadrants that contain most of the endpoints
$\mathnormal{E}$ indicate the orientation of the seam. The mean value of
$\mathnormal{E}$ in these two quadrants is used to solve the seam
characteristic line $\mathnormal{l}$. By mapping the depth information in
$\mathbf{P}$ related to the points in $\mathnormal{l}$, the respective 3D
coordinates $\zeta=\\{\zeta_{1},\ldots,\zeta_{n}\\}$ can be obtained. Segments
are assigned to $\mathnormal{l}$ based on the distance between ${\zeta_{1}}$
and ${\zeta_{n}}$. A point $Q$ is added to visually indicate in the
interface’s screen the target to be followed by the user during the task.

###  II-C Electric Arc Localization

Figure 6: Overview of the electric arc and its center during the process.

Once the welding process starts, a strong light is generated by the electric
arc [6], which is perceived through a helmet as a small bright spot; An
effective welding task requires to move this spot in the correct direction,
i.e., along the seam. By computing the arc’s center (shown in Fig. 6), the
intersection point between the center of the electrode and the seam can be
located by our multi-sensor system. As the intensity of the HDR image is
affected by the strong light, thus, only the region around the seam is
processed and displayed to the user in real-time. Areas far from the seam are
statically displayed based on the scene’s initial observation (i.e., before
the welding process began). The region around the electric arc has high-
intensity levels in the captured images, thus, this area of interest is
extracted with a binarization approach that generates a feature image.
Nevertheless, some light reflections may still appear in the binary map as its
computation is susceptible to noises [26].

####  II-C1 Light Intensity-Based Confidence Map

To remove flashing noises and minimize the computational time, a confidence
map is constructed for each image [28, 29, 30]. A light intensity-based
confidence (LIC) map is applied to the HDR image, where the confidence is
determined by the past trajectory and the present observations, see Fig. 7.
The frame is first divided into multiple tiles, with an average light
intensity $\bar{\xi}$ computed for each of them. The light emitted from the
welding spot makes the light intensity of the surroundings increase radially
from the center of the electric arc. As we consider the slow motions of the
welding torch by the user, it is reasonable to assume that for a continuous
high intensity in a tile for the past few frames, there is an exponentially
high likelihood that an equal or greater light intensity value will be
obtained in that same tile in the next frame. Similarly, for a continuous low
intensity in the past few frames, it is unlikely there will be an electric arc
shown in the coming frames.

Our method aims to reduce the probability density distribution of the whole
image. Only high values are considered to create a high contrast image; Thus,
an exponential growth/decline is implemented to model the likelihood for the
welding spot to be located in a tile of a frame. The confidence level of each
grid is then calculated with the following function:

${p}[t]=\text{min}\bigg{(}\frac{(b\cdot\text{norm}(\bar{\xi}))^{p[{t-1}]}}{b}(\sigma+p[{t-1}]),1\bigg{)}$
(2)

where $\text{norm}(\bar{\xi})={\bar{\xi}}/{255}$ represents the normalized
average intensity in each tile; $b>0$ denotes a user-defined parameter to
specify the desired intensity range; $\sigma>0$ controls the sensitivity to
the current observation, where small values result in delayed responses and
large values make it susceptible to fast changes (a good compromise is
$\sigma=1$, which is used for the rest of our derivations). The current
probability (at time instance $t$) is denoted by $p[t]$, which describes the
confidence of a welding spot to be found in a tile; The method is initialized
as $p[t-1]=0$ for $t=0$. The minimum value in (2) ensures that $0\leq p[t]\leq
1$ is satisfied.

Different torch poses may affect the light intensity values, making $p[t]$ to
fluctuate. To achieve stable and consistent values, the confidence is then
normalized $p[t]/P_{max}$ based on the maximum $p[t]$ value of the image
$P_{max}$. Using this confidence map, we can estimate the location of the
welding spot from the high confidence areas (red regions in Fig. 7). We can
also predict the next location of the electric arc, labeled in green color in
Fig. 7.

Figure 7: Average light intensity $\bar{\xi}$ is calculated for each grid over
the image. An updated LIC ${p[t]}$ is solved for each grid based on the past
confidence level $p[t-1]$ and the current light intensity $\bar{\xi}$ value.
Three colors are used to visualize the probability of the potential welding
spot location. The red color tile indicates it has the highest confidence
level; Green indicates high confidence level; Blue indicates low confidence.

####  II-C2 Dimension Filter

To speed-up the computation time, regions extracted with the binarization
approach that overlap with a high probability area are only considered for
tracking the welding spot. To this end, we first represent these regions by a
series contours $K=\\{K_{1},\ldots,K_{n}\\}$ computed from the binary image.
Then, the areas $\alpha=\\{\alpha_{1},\ldots,\alpha_{n}\\}$ of the contours
$K$ are computed by using the Shoelace method
$\alpha_{i}=\frac{1}{2}\left|{u}-\hat{u}\right|$, for $u$ and $\hat{u}$
defined as follows [31]:

${u}=q_{nx}q_{iy}+\sum_{i=1}^{n-1}q_{ix}q_{(i+1)y},\quad\hat{u}=q_{ny}q_{ix}+\sum_{i=1}^{n-1}q_{iy}q_{(i+1)x}$
(3)

where $q=\\{q_{1},\ldots q_{m}\\}$ denotes the set of $m$ image points
$q_{j}=[q_{jx},q_{jy}]$ of the $i$th contour $\mathnormal{K}_{i}$. We then
evaluate the dimension of the areas $\alpha_{i}$ and apply the following
condition:

${\tilde{K}}_{i}=\begin{cases}K_{i},&\text{if
$A_{min}<\alpha_{i}<A_{max}$}\end{cases}$ (4)

where ${\tilde{K}}=\\{{\tilde{K}}_{1},\ldots,{\tilde{K}}_{f}\\}$ denote the
_valid_ contours that can be used in the tracking algorithm, for
$A_{min},A_{max}>0$ as scalars to determine this selection.

Figure 8: (a) The welding torch blocks the view; (b) Some light reflections
from the surroundings share a similar light intensity level; (c) A LIC map is
generated; (d) Noises in (b) are filtered out; The center of the contour
labeled in pink is away from the actual center of the electrode; (e) The
correct center can be found with the proposed method; (f) Comparison of the
center before and after the use of the method; (g) Demonstration of the
verification process. Figure 9: The conceptual interface of the XR bot trainer
where $Q$ is the recommended point and $C$ the current welding spot.

####  II-C3 Electric Arc Tracking

During the task, the user may partially block the view of the electric arc
with the torch. This may lead to a shift between the contour’s center and the
actual center of the arc. To deal with this issue, our method uses a circle
that includes all contour points $q$ in $\mathnormal{\alpha}$ with a minimum
radius $r_{i}$ available. To this end, we model the center of the bounding
circle by $C_{i}=[C_{ix},C_{iy}]$. As the region that indicates the electric
arc generally produces the largest bounding circle, our algorithm verifies
these regions starting from the largest circle. However, noises in the image
may give rise to the false detection of similar-sized contours, which creates
ambiguity for the system.

Algorithm 1 presents the proposed verification process, which is depicted in
Fig. 8. All valid centers $\mathnormal{C_{i}}$ and its radius
$\mathnormal{r_{i}}$ found from the contour $\mathnormal{{\tilde{K_{i}}}}$ are
sorted based on the size of $\mathnormal{r_{i}}$ in descending order, i.e.,
$\mathnormal{C}_{1}$ denotes the center of the largest circle. By comparing
the distance between the $\mathnormal{C_{i}}$ with the previous state location
$\mathnormal{C_{prev}}$, the circle can be verified. Through iterating each
$\mathnormal{C_{i}}$, the largest reasonable circle around
$\mathnormal{C_{prev}}$ can be located. Kalman filter [32] is applied to
further verify the current estimation. The point $\mathnormal{C(x,y)}$ denotes
the detected center of the bounding circle containing the welding spot. Fig. 9
conceptually depicts the usage of this point in the proposed interface.

Input: $C_{prev}$, ${\tilde{K}}$ $\triangleright$ Previous center, current
contours

Output: $C$, $r$ $\triangleright$ Center and radius of circle

$N\leftarrow\text{Compute\\_Length\\_of\\_Array}({\tilde{K}})$

for  _$i\leftarrow 1$ to $N-1$_ do

$C_{i}$, $r_{i}$ = Find_Minimum_Circle($\tilde{K}_{i}$)

Sort_Cicles_by_Maximum_Radius($C$, $r$)

if  _$C_{prev}$ is not empty_ then

$i\leftarrow 1$

$valid\leftarrow False$

while  _$i\leq N$ & $valid$ is $False$_ do

$valid$ = Compare_Center_Distance($C_{i}$, $C_{prev}$)

if  _$valid$_ then

$C_{prev}\leftarrow C_{i}$

return $C_{i}$, $r_{i}$

Algorithm 1 Electric Arc Tracking.

###  II-D Motion Direction Estimation

The direction of the welding movement depends on the handedness of a user and
the location of the seam. For example, a right-hand welder will usually
perform right to left motions in a horizontal workpiece; Vertical weld motions
are typically from bottom to top. To visually determine the direction of this
welding spot, our automated system records the movement of the center
$\mathnormal{C}$ for 20 frames. The starting point $\mathnormal{C_{start}}$ is
computed based on sensory data captured from these first few frames. The
moving direction is determined by comparing the signed coordinate difference
between the starting point $\mathnormal{C_{start}}$ and the current point
$\mathnormal{C}$. For instance, if $\mathnormal{C}$ is on the left side of
$\mathnormal{C_{start}}$ after 20-frames, it implies that the welder is going
from left to right. This way, guidance from the XR bot trainer system can be
adjusted accordingly to fit the users’ torch motions.

Figure 10: Overview of the welding training system. After we place a
workpiece, RGB-D and HDR images are captured. The groove detection and the
seam localization algorithms are executed to find the welding path. The
electric arc localization is applied to the HDR image to compute the center of
the arc. The yellow path indicates the welding trajectory with the seam line
shown in blue. Images and results are registered into a 2D frame and uploaded
to the server. Electronic devices, like a VR headset, are used to access the
real-time output via a webpage.

###  II-E Virtual XR Bot Trainer

The proposed XR bot trainer is responsible for providing instant welding
guidance to the user via a head-mounted VR display. The aim of the interface
is to provide trainees with the path to be followed by the torch as well as
its desired speed. For that, the target point $\mathnormal{Q=[Q_{x},Q_{y}]}$
is displayed into the VR headset as a moving circle that must be followed by
the user by closely placing the electric arc’s center $C$ over it. By
comparing the coordinates of the target and feedback points, we can compute
the _relative_ motion of the torch as $\vec{v}=Q-C$, and thus, generate visual
cues for the user. For example, if the position and velocity are close to the
bot’s suggestion, the point $\mathnormal{Q}$ is displayed with a green circle;
Deviations from the target location are displayed with either a red circle or
a blue circle, depending on whether $C$ lags or leads the point $Q$.

The instantaneous position error $\|Q-C\|$ and the average error
$\overline{\epsilon}$ are calculated by the XR bot trainer based on the user’s
motions; This information is displayed in real-time. After the completion of
the task, the system generates a final score and plots the executed welding
trajectory. The average error $\overline{\epsilon}$ is calculated as follows:

$\overline{\epsilon}=\frac{1}{N}{\sum_{i=1}^{N}\bigg{(}{|C_{xi}-Q_{xi}|}/{Q_{xi}}+{|C_{yi}-Q_{yi}|}/{Q_{yi}}}\bigg{)}$
(5)

where $N$ is the total number of target points $Q$ involved from
$\mathnormal{C_{start}}$ to $\mathnormal{C}$. With this average error, the
system then quantifies the task’s final score as
$\gamma=100(1-\overline{\epsilon})$.

As the proposed method is aimed at lowering the entry barrier for welding, it
supports cross-platform usage. The real-time results and the 3D model of the
workpiece are available on a webpage, which can be accessed through any
electronic device via a wireless connection. This valuable feature enables the
instructor and participants to observe the live welding process through a
display device.

##  III Results

Figure 11: Experiment setup of the proposed XR welding training system.

###  III-A Experimental Setup

To validate the proposed methodology, a prototype XR interface is built and
instrumented with vision sensors to capture the whole welding process. The
data processing pipeline is shown in Fig. 10, which contains an RGB-D sensor
(Intel Realsense SR305), an HDR camera (New Imaging Technology MC1003), a VR
headset with a controller (Oculus Quest 2), a PC (GPU RTX 3060), and a
standard TIG welding torch. Experiments are conducted to evaluate the system,
with two groups of trainees (who have no any prior welding experience) being
asked to perform several tasks. One group learns with the traditional (i.e.
current practice) approach while the other group learns with the proposed
multi-sensor interface. In the control group, trainees practice with the
instructor and receive one-to-one guidance when necessary; In the experimental
group, after the instructor gives a sample demonstration, the trainees
practice with the interface and receive no human guidance. Success and failure
are evaluated by the human instructor at the end of each task. Experiments
with the interface further evaluate: (1) Accuracy of the welding path, (2)
accuracy of the arc location, and (3) effectiveness of the use of this system
in teaching and learning.

###  III-B Experimental Design

A horizontal fillet TIG welding is chosen for all experiments, see Fig. 11. A
non-consumable electrode is used to generate the electric arc. The weld and
the molten pool are shielded from environmental contamination by inert
shielding gas [6]. Mild steel workpieces of $10$ mm thickness and a straight
welding path are used. Although a fillet weld workpiece is prepared, trainees
are not required to add any filler at the beginning of their training until
they are capable to do so. To obtain a fair comparison between the two methods
in terms of learning effectiveness, the task is judged as success/failure by
the human instructor. If success is obtained before filler inclusion, then,
the trainee can start adding filler in the next trial. If success is reached
with filler inclusion, the trainee is considered to have mastered the welding
skill and the practice is completed for the trainee.

In the experiments, the system first records the traditional welding process
that is conducted by the instructor. Then, the technique is explained to the
trainees from these video recordings. Afterwards, trainees perform the welding
task by using the interface; The instructor and other participants observe the
process from a display in _real-time_. After completion, the multimodal
recordings can be replayed to analyze various details of the trainees’
execution. As proficiency in the task is directly related to the number of
repetitions a user has performed [33, 34, 35], we quantify the system’s
effectiveness by computing a learning curve that relates the trial number with
the obtained score.

###  III-C Performance Analysis of the Multi-Sensor Interface

Figure 12: Various workpieces, where blue lines indicate the identified seam
location. (a)-(b) Fillet workpieces; (c) Butt workpieces.

As different groove geometries may lead to different seam localization
results, we conduct experiments to validate the performance of the automatic
seam detection algorithm, see Fig. 12. Fillet and butt workpieces with various
orientations are used for these validation tests, which show that the
algorithm can correctly locate different types of seams with various poses. We
also conducted experiments to validate the performance of the spot
localization and denoising algorithms; We did it in a room with multiple
_unknown_ light sources to test its robustness. The LIC map is applied to the
raw HDR image for noise removal, with $b=4.5$ in (2). We define tiles with
$p\geq 0.65$ as areas with high probability, and $p\geq 0.95$ as the estimated
welding spot region. These two types of regions are depicted with green and
red tiles in Fig. 13, where we can see that the LIC map predicts the location
of the welding spot from past and current intensity measurements.

In our performance experiments, the welding torch is arranged at various
angles and configurations, which sometimes occlude the electric arc. Fig. 14
depicts configurations where the welding spot is not entirely captured by the
HDR camera. The minimum radius circle $\mathnormal{r}$ is formed, and thus,
the center $\mathnormal{C}$ can still be computed in these situations. Based
on the measured trajectories of the electric arc, the bot assistant provides
_real-time_ suggestions for user, as depicted in Fig. 15. When the trainee has
a slow pace, the suggested target point is displayed in red to alert lagging
motion. When the speed of welding point is too high (which results in
insufficient time to melt the material) the circle is displayed in blue. When
the welding point is close to the target, the circle is shown in green. These
features of our system provide users with valuable guidance to help them
improve their technique.

Figure 13: Visualization of the welding spot prediction that is generated by
the LIC map. (a) The orange arrow indicates the welding direction; (b) Red and
Green tiles indicate the potential spot location, from high to medium. Blue
indicates low probability areas. Figure 14: Visualization of the contours $K$,
displayed in blue color. Red circles indicate the radius $\mathnormal{r}$ and
the center $\mathnormal{C}$.  Figure 15: (a)-(c) Examples of the welding path
and the guidance support provided by the XR bot trainer; The orange and white
circles indicate the start and end points of the seam. (d) Welding trajectory,
errors, and score of the conducted task; (e) User’s view in the VR headset.
Figure 16: Comparison between the proposed LIC Map method and the traditional
Softmax method in estimating the electric arc location. The black cross
indicates the center of the electric arc. The green tile represents a medium-
high belief ($\geq 0.65$). The red tile represents a high probability value
($\geq 0.95$). (a)–(b) The confidence value of the surrounding in the Softmax
method is higher. (c) A larger error is obtained in estimating the center with
the Softmax method. Figure 17: Evolution of the probability computed with the
LIC Map method and the Softmax method for one tile. The green and red
horizontal lines indicate the probability level at $0.65$ and $0.95$. The pink
area is the period when noise occurs. A steeper slope and a narrower density
are obtained by our proposed method. Noise is also greatly suppressed with our
exponential decline function. Figure 18: Comparison of the histograms obtained
with the proposed LIC Map method and the Softmax method. More tiles with $\geq
0.65$ probability are found in the Softmax approach, which indicates a higher
uncertainty in the estimation. Figure 19: Comparison between the proposed
method with the state-of-the-art practice in electric arc localization. The
red dot is the center determined by our electric arc localization algorithm.
The dark blue dot is the center found with a contour-based approach. The green
dot is the center identified by the intensity-based approach. Figure 20:
Quantitative comparison between the proposed method with the state-of-the-art
practice in electric arc localization. (a) The welding process with results
overlaid; (b) The position of the electric arc is located by three approaches.
Multiple results are obtained from the light intensity-based method.

###  III-D System Comparison

####  III-D1 LIC Map

There are other prediction methods for estimating features similar to the
welding spot. The Softmax function estimates the belief based on the current
information and the weighted past events [36, 37]. In contrast with this
approach, our method (2) provides exponential changes in the estimated belief.
We compare our method with the classical Softmax [36, 37]:

$p[t]=\text{norm}(\bar{\xi})+w\cdot p[t-1]$ (6)

where $\text{norm}(\bar{\xi})$ is the normalized light intensity at the
current time instance, and $w>0$ is a scalar weight (which we set to
$w=0.01$). The prediction performance of these methods can be visualized in
Fig. 16, which shows that the tiles around the welding spot (i.e. the red
tile) present a larger (undesired) belief in the Softmax approach than in our
method.

To get an insight into the effect of our exponential approach, in Fig. 17 we
compare the belief changes of a grid when the welding spot passes by it. The
figure shows that the confidence level increases and decreases rapidly when
the arc approaches the tile; When noises appear, our method can suppress their
effect and maintain a continuous low confidence level. The developed LIC map
produces a narrow probability distribution (see Fig. 18) that results in fewer
potential regions for the spot localization (green and red tiles); This leads
to a faster denoise process. With our exponential approach, we can sharpen the
differences between the useful welding light area and the rest of the image.

####  III-D2 Welding Spot Determination

There are various state-of-the-art approaches to locate the center of the
electric arc. Contour-based and light intensity-based methods from the OpenCV
library can be used to determine this spot. Fig. 19 compares the performance
of these methods with our approach. For fairness sake, the same input data is
used for all methods, as a binary image generated by thresholding. In the
contour-based method, the center is determined based on the largest contour’s
center in an image. In the light intensity-based approach, the center is found
from all the present high light intensity regions without considering the
dimension.

The main difference between these three methods is the accuracy in determining
the actual welding spot with unpredictable noise inclusion. In Fig. 19 and
Fig. 20, the contour-based approach misidentifies the starting point of the
electric arc as the welding spot. In the light intensity-based approach, more
than one “welding spot” are found. With this method, extra image processing
for arc localization is normally required. A quantitative comparison between
the three method is presented in Fig. 20. A significant error can be observed
in the state-of-the-art approaches. This proves that the proposed algorithm
can effective determine the the welding spot that is manipulated by the user
in noisy environments.

###  III-E Effectiveness in Teaching and Learning

As the trainees in the _control_ group have no prior experience with the XR
bot trainer, thus, poor welding performances is observed at first. As they
familiarize with the system, it becomes easier to perform the welding task;
Filler is successfully added to the weld (a sign of skill mastery) within a
short time of hands-on testing, see Fig. 21. The performance of the executed
welding trajectory is shown in Fig. 22.

Figure 21: Examples of the trainees’s performance by independently conducting
welding tests. (a) Standard learning approach without filler after 5 trials;
(b) Proposed method without filler after 5 trials; (c) Proposed method after
11 trials; (a) Standard learning approach led to consistent failures to
maintain a proper movement speed; (b)–(c) Proposed method helped trainees to
properly learn the technique, and thus, to master the complete TIG welding
skill. Figure 22: Sample welding trajectories performed by the trainees with
the aid of our XR system. (a) Spatial x-y motions of the torch, captured
during 3 different trials; (b) Instant errors of the trajectories
corresponding to the motions in (a). Figure 23: Average scores and errors of
different welding trials conducted by the group of trainees with our system.

Compared with the standard instructional practice, it is easier for
instructors to explain the welding process and techniques to beginners with
the proposed interface. In terms of teaching efficiency, the instructor has to
keep monitoring and providing guidance to each student one by one in the
traditional approach, which is a time-consuming task for them. In the
_experimental_ group that uses the interface, the instructor is only required
to evaluate the task as success or failure once it is completed. This greatly
reduces the instructor’s workload since minimal human support is needed. It
can be concluded that our XR system can lead to a boost in teaching
efficiency. This new approach also enables to increase the number of trainees
that can participate in one session.

With respect to the learning effectiveness, it took around 5 hours for
trainees in the control group to master the welding skill while only 3 hours
in the experimental group. On average, around 11 trials were needed for a
trainee to master the TIG welding skill with our new approach; The multi-
sensor interface enabled these participants to _independently_ conduct the
training tasks. Around 6 more trials were needed when using the traditional
approach, with most participants showing a high dependency on the instructor.
The horizontal axis in Fig. 23 indicates the number of trials conducted by the
trainees in the experimental group. The left vertical axis represents the
score obtained during the corresponding trial; The right vertical axis shows
the computed average error. As some received success at the eighth trial, no
more practices are conducted. These results demonstrate that trainees’
proficiency in the task consistently improves with the number of trials.
Effectiveness in learning grows with the use of the developed XR bot trainer.

Compared with other VR training studies, our automated XR bot assistant
provides immediate _real-time_ support to the user, during real-world welding
practice. The XR interface provides motion guidance and performance metrics
based on sensor feedback. Despite this instant support, our system cannot
update the seam’s geometric information in real-time since the welding path is
computed based on the initial (static) point cloud of the workpiece. However,
the instant performance of trainees can be accessed via various types of
displays that the instructor can remotely supervise. As virtual welding is
only capable of providing a synthetic perceptual experience, trainees cannot
address the typical psychological adaptation that they gain through practice
in the field.

In the accompanying multimedia file, we demonstrate the performance of the
system with multiple experimental videos. https://github.com/romi-
lab/VR_Welding/raw/main/video.mp4

TABLE I: Comparison with the features and functions in state-of-the-art VR methods for welding training. | Virtual Welding | Real Welding
---|---|---
Methods |  Seam & Arc Localization |  Guidance |  Seam & Arc Localization |  Guidance
[9] |  ✓ |  ✓ |  - |  -
[11] |  ✓ |  ✓ |  - |  -
[12] |  ✓ |  ✓ |  - |  -
[13] |  ✓ |  ✓ |  - |  -
Ours |  - |  - |  ✓ |  ✓
Figure 24: Some faulty examples. (a)–(b) Seam between sheet metal workpieces
cannot be located as the depth difference is small; (c) Due to the limited
angle of view from the system, the welding torch is hard to control with the
system.

##  IV Conclusion

In this paper, we presented an automated XR training bot assistant for
teaching and learning arc welding tasks. It involves the use of a welding
torch, an RGBD camera, an HDR camera, a VR headset, and image processing
algorithms. By using this multi-sensor interface, the seam can be
automatically located with 3D vision. The instant welding spot from the
electric arc is recognized and immediate XR advice is provided to the user to
improve the technique. Task scores, errors, and paths are displayed onto the
interface to provide the user with valuable feedback information. The
effectiveness of the proposed method is experimentally validated with a group
of beginners. Compared to the current practice, our method allows the
instructor to have a clearer and more convenient way to demonstrate the
process, and to quantify the trainees performance. With the multi-sensor
interface, users can independently practice the skill with minimal human
supervision, hence, the proposed method has the potential to increase the
number of trainees participating in a single session.

There are various limitations of the developed system. For example, only thick
materials can be considered at present. A thin layer such as sheet metal is
too small to be recognized, as the depth difference in the groove is
insignificant for the algorithm to locate the seam and it may be misjudged as
one large workpiece, see e.g. Fig. 24; The current system works for grooves of
least 5-mm thickness. The interface may also pose difficulties in controlling
the angle of the welding torch due to the limited camera view, see e.g., in
Fig. 24(c). This may lead to failures in joining two workpieces without a
filler.

Future work includes the integration of more sensing devices (e.g., 3D
trackers) to robustly detect the pose of the torch. A hierarchical ranking
system may be introduced to classify the users’ proficiency in the task based
on quantitative metrics [38]. Data collected from the practice can be used to
customize the bot trainer and to provide tailor-made learning exercises. This
can help trainees to acquire the skill faster through optimizing the learning
stages. The valuable feedback information of the interface can also be used to
guide the motion of robotic welder. A motorized robot arm is currently being
developed by our team to automatically perform the task based on a simple
multi-sensor system.

## References

  * [1] P. K. Ghosh, _Pulse current gas metal arc welding_. Springer, 2017.
  * [2] M. I. Khan, _Welding science and technology_. New Age International, 2007.
  * [3] C. Latella, Y. Tirupachuri, L. Tagliapietra, L. Rapetti, B. Schirrmeister, J. Bornmann, D. Gorjan, J. Čamernik, P. Maurice, L. Fritzsche _et al._ , “Analysis of human whole-body joint torques during overhead work with a passive exoskeleton,” _IEEE Transactions on Human-Machine Systems_ , 2021.
  * [4] Y. Liu, W. Zhang, and Y. M. Zhang, “A tutorial on learning human welder’s behavior: Sensing, modeling, and control,” _Journal of Manufacturing processes_ , vol. 16, no. 1, pp. 123–136, 2014.
  * [5] S.-B. Asplund and N. Kilbrink, “Lessons from the welding booth: theories in practice in vocational education,” _Empirical Research in Vocational Education and Training_ , vol. 12, no. 1, pp. 1–23, 2020.
  * [6] K. Weman, _Welding processes handbook_. Elsevier, 2011.
  * [7] J. M. Antonini, “Health effects of welding,” _Critical reviews in toxicology_ , vol. 33, no. 1, pp. 61–103, 2003.
  * [8] A. D. Althouse, C. H. Turnquist, W. A. Bowditch, K. E. Bowditch, and M. A. Bowditch, _Modern welding_. Goodheart-Wilcox Publisher, 2004.
  * [9] C.-c. Chung, C.-C. Tung, and S.-J. Lou, “Research on optimization of vr welding course development with anp and satisfaction evaluation,” _Electronics_ , vol. 9, no. 10, p. 1673, 2020.
  * [10] C.-Y. Huang, S.-J. Lou, Y.-M. Cheng, and C.-C. Chung, “Research on teaching a welding implementation course assisted by sustainable virtual reality technology,” _Sustainability_ , vol. 12, no. 23, p. 10044, 2020.
  * [11] M. I. M. Isham, H. N. H. Haron, F. bin Mohamed, C. V. Siang, M. K. Mokhtar, and A. S. binti Azizo, “Mobile vr and marker tracking method applied in virtual welding simulation kit for welding training,” in _2020 6th International Conference on Interactive Digital Media (ICIDM)_. IEEE, 2020, pp. 1–5.
  * [12] U. Yang, G. A. Lee, Y. Kim, D. Jo, J. Choi, and K.-H. Kim, “Virtual reality based welding training simulator with 3d multimodal interaction,” in _2010 International Conference on Cyberworlds_. IEEE, 2010, pp. 150–154.
  * [13] K. Kobayashi, S. Ishigame, and H. Kato, “Skill training system of manual arc welding,” in _Entertainment Computing_. Springer, 2003, pp. 389–396.
  * [14] R. S. Vergel, P. M. Tena, S. C. Yrurzum, and C. Cruz-Neira, “A comparative evaluation of a virtual reality table and a hololens-based augmented reality system for anatomy training,” _IEEE Transactions on Human-Machine Systems_ , vol. 50, no. 4, pp. 337–348, 2020.
  * [15] R. Agrawal and J. S. Pillai, “Augmented reality application in vocational education: A case of welding training,” in _Companion Proceedings of the 2020 Conference on Interactive Surfaces and Spaces_ , 2020, pp. 23–27.
  * [16] K. Wu, J. Tan, H. L. Xia, and C. B. Liu, “An exposure fusion-based structured light approach for the 3d measurement of a specular surface,” _IEEE Sensors Journal_ , vol. 21, no. 5, pp. 6314–6324, 2020.
  * [17] M. Purohit, M. Singh, A. Kumar, and B. K. Kaushik, “Enhancing the surveillance detection range of image sensors using hdr techniques,” _IEEE Sensors Journal_ , vol. 21, no. 17, pp. 19 516–19 528, 2021.
  * [18] H. Yu, Q. Fu, Z. Yang, L. Tan, W. Sun, and M. Sun, “Robust robot pose estimation for challenging scenes with an rgb-d camera,” _IEEE Sensors Journal_ , vol. 19, no. 6, pp. 2217–2229, 2018.
  * [19] S. Huo, B. Zhang, M. Muddassir, D. T. Chik, and D. Navarro-Alarcon, “A sensor-based robotic line scan system with adaptive roi for inspection of defects over convex free-form specular surfaces,” _IEEE Sensors Journal_ , 2021.
  * [20] L. Yang, E. Li, T. Long, J. Fan, and Z. Liang, “A novel 3-d path extraction method for arc welding robot based on stereo structured light sensor,” _IEEE Sensors Journal_ , vol. 19, no. 2, pp. 763–773, 2018.
  * [21] Y. Tian, H. Liu, L. Li, G. Yuan, J. Feng, Y. Chen, and W. Wang, “Automatic identification of multi-type weld seam based on vision sensor with silhouette-mapping,” _IEEE Sensors Journal_ , vol. 21, no. 4, pp. 5402–5412, 2020.
  * [22] P. Zhou, R. Peng, M. Xu, V. Wu, and D. Navarro-Alarcon, “Path planning with automatic seam extraction over point cloud models for robotic arc welding,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 3, pp. 5002–5009, 2021\.
  * [23] R. Peng, D. Navarro-Alarcon, V. Wu, and W. Yang, “A point cloud-based method for automatic groove detection and trajectory generation of robotic arc welding tasks,” _2020 17th International Conference on Ubiquitous Robots (UR)_ , pp. 380–386, 2020.
  * [24] A. W. Society, “Aws a3. 0m/a3. 0: 2010: standard welding terms and definition,” 2010.
  * [25] L. Ding and A. Goshtasby, “On the canny edge detector,” _Pattern Recognition_ , vol. 34, no. 3, pp. 721–725, 2001.
  * [26] X. Xiong, G. Jin, Q. Xu, H. Zhang, and J. Xu, “Robust line detection of synthetic aperture radar images based on vector radon transformation,” _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , vol. 12, no. 12, pp. 5310–5320, 2019.
  * [27] J. L. Bentley, “Multidimensional binary search trees used for associative searching,” _Communications of the ACM_ , vol. 18, no. 9, pp. 509–517, 1975\.
  * [28] Y. Fang, S. Ko, and G.-S. Jo, “Robust visual tracking based on global-and-local search with confidence reliability estimation,” _Neurocomputing_ , vol. 367, pp. 273–286, 2019.
  * [29] R. Yasarla and V. M. Patel, “Confidence measure guided single image de-raining,” _IEEE Transactions on Image Processing_ , vol. 29, pp. 4544–4555, 2020.
  * [30] H. Idrees, I. Saleemi, C. Seibert, and M. Shah, “Multi-source multi-scale counting in extremely dense crowd images,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2013, pp. 2547–2554.
  * [31] B. Polster, _The shoelace book: a mathematical guide to the best (and worst) ways to lace your shoes_. American Mathematical Soc., 2006, no. 24.
  * [32] B. Millidge, A. Tschantz, A. Seth, and C. Buckley, “Neural kalman filtering,” _arXiv preprint arXiv:2102.10021_ , 2021.
  * [33] M. J. Anzanello and F. S. Fogliatto, “Learning curve models and applications: Literature review and research directions,” _International Journal of Industrial Ergonomics_ , vol. 41, no. 5, pp. 573–583, 2011.
  * [34] G. Fioretti, “The organizational learning curve,” _European Journal of Operational Research_ , vol. 177, no. 3, pp. 1375–1384, 2007.
  * [35] A. Newell and P. S. Rosenbloom, “Mechanisms of skill acquisition and the law of practice,” _Cognitive skills and their acquisition_ , vol. 1, no. 1981, pp. 1–55, 1981.
  * [36] E. Jang, S. Gu, and B. Poole, “Categorical reparameterization with gumbel-softmax,” _arXiv preprint arXiv:1611.01144_ , 2016.
  * [37] F. Wang, J. Cheng, W. Liu, and H. Liu, “Additive margin softmax for face verification,” _IEEE Signal Processing Letters_ , vol. 25, no. 7, pp. 926–930, 2018.
  * [38] M. Lorenzini, W. Kim, and A. Ajoudani, “An online multi-index approach to human ergonomics assessment in the workplace,” _IEEE Transactions on Human-Machine Systems_ , 2022.


# Single-step first order phase transition and gravitational waves in a SIMP
dark matter scenario

Nabarun Chakrabarty<EMAIL_ADDRESS>Department of Physics, Indian
Institute of Technology Kanpur, Kanpur, Uttar Pradesh-208016, India  Himadri
Roy<EMAIL_ADDRESS>Department of Physics, Indian Institute of
Technology Kanpur, Kanpur, Uttar Pradesh-208016, India  Tripurari Srivastava
<EMAIL_ADDRESS>Department of Physics and Astrophysics, University
of Delhi, Delhi 110007, India

###### Abstract

We investigate the non-zero temperature dynamics of a sub-GeV dark matter
scenario freezing-out via self-interactions. As a prototype, we take up the
case of a scalar dark matter species undergoing $3\to 2$ number changing
annihilations catalysed by another scalar. We study the shape of the thermal
potential of this scenario in a parameter region accounting for the observed
relic abundance. An analysis reveals the possibility of a first order phase
transition with bubble nucleation occurring at sub-GeV temperatures. This
finding can be correlated with the typical sub-GeV masses in the framework.
The gravitational wave spectra associated with such a phase transition is
subsequently computed.

##  I Introduction

It has been concurred by numerous observations Rubin:1970zza ; Bertone:2004pz
; Hu:2001bc ; WMAP:2012nax that around a fourth of the total energy budget of
the universe is some non-luminous matter known as dark matter (DM). The amount
of such matter in the universe is also well measured from various
observations, e.g., cosmic microwave background (CMB) Planck:2018vyg and large
scale structure surveys Blumenthal:1984bp . If DM is thought to be an
elementary particle, then there is no information available from the current
experiments on its mass and quantum numbers. More importantly, that there is
no such particle candidate for DM within the Standard Model (SM), advocates
the presence of additional dynamics. While DM can be hypothesized to stem from
some beyond-the-Standard Model (BSM) framework, a possibility in this
direction is that DM is a weakly interacting massive particle (WIMP)
Kolb:1990vq ; Jungman:1995df ; Bertone:2004pz ; Feng:2010gw ; Arcadi:2017kky ;
Roszkowski:2017nbc . That is, the DM-SM interactions are in the _weak_
ballpark. Moreover, the DM relic abundance in the WIMP paradigm is generated
typically through $2_{\text{DM}}\to 2_{\text{SM}}$ annihilation processes
followed by a thermal _freeze out_.

Search for DM through its scattering with nuclei at terrestrial detectors is
known as direct detection XENON:2018voc ; LUX:2017ree . However, non-
observation of such scattering at have put stringent upper limits on the
corresponding cross sections, and ultimately, on the DM-SM interactions. Such
small interactions are also consistent with the null results obtained in DM
searches at the energy frontier, the Large Hadron Collider (LHC) being the
prime example here Abercrombie:2015wmb . These results have pushed a plethora
of WIMP scenarios to a corner Roszkowski:2017nbc ; Arcadi:2017kky . An
interesting alternative therefore could a strongly interacting massive
particle (SIMP) Bernal:2015xba ; Choi:2016hid ; Choi:2016tkj ; Bernal:2017mqb
; Chauhan:2017eck ; Mohanty:2019drv ; Bhattacharya:2019mmy ; Lee:2015gsa ;
Yamanaka:2015tba ; Hochberg:2015vrg ; Hochberg:2014kqa ; Choi:2017zww ;
Choi:2018jsb ; Herms:2018ajr ; Ma:2015mjd . In this case, the freeze-out
dynamics is driven by DM matter _number changing_ processes, i.e., processes
that feature an imbalance of DM particles in the initial and final states. The
main advantage of the SIMP paradigm is that it allows such number changing
processes to be predominantly driven by DM self interactions. The DM-SM
interaction stengths are thus permitted to be appropriately small so as to
conform with the limits from direct detection and colliders. A SIMP therefore
typically has sub-GeV mass and a large self-scattering cross section.

With the advent of gravitational wave (GW) interferometry Thrane:2013oya ;
Sathyaprakash:2012jk ; LISA:2017pwj ; Crowder:2005nr ; Seto:2001qf ;
Sato:2017dkf ; Lentati:2015qwp ; Janssen:2014dka ; Fairbairn:2019xog , it has
become possible to probe a dark sector through its GW imprints
Breitbach:2018ddu ; Schwaller:2015tja ; Jaeckel:2016jlh ; Addazi:2016fbj ;
Baldes:2017rcu ; Tsumura:2017knk ; Baldes:2018emh ; Croon:2018erz . In this
study, we aim to study the possibility of a strong first order phase
transition (SFOPT) Kamionkowski:1993fg ; Witten:1984rs ; Hogan:1986qda ;
Grojean:2006bp and the GW spectrum from associated with a sub-GeV SIMP
framework. The model we take up is where the SM is extended by the scalars
$\phi$ and $\delta$ that are singlets under the SM gauge symmetry
Mohanty:2019drv . The scalar $\phi$, stabilised by a governing
$\mathbb{Z}_{2}$ symmetry, becomes the DM candidate here. The number changing
process for this case that is consistent with the $\mathbb{Z}_{2}$ is
$\phi\phi\phi\to\phi h_{2}$, where $h_{2}$ is an admixture of $\delta$ and the
Higgs doublet $H$. The scalar $\delta$ and its interaction with the DM $\phi$
thus play crucial roles in generating the observed relic density through the
aforesaid $3\to 2$ dynamics. Given the importance of $\delta$ in this setup, a
pertinent question to ask is whether the same can trigger a SFOPT and how
strong is the resulting GW amplitude. We plan to explore this possibility here
by incorporating finite temperature corrections to the scalar potential along
the direction of $\delta$.

The study is organised as follows. We introduce the setup and outline the
$3\to 2$ dynamics in section II. Thermal corrections to the scalar potential
are computed in section III. The same section presents the GW amplitude for
representative benchmark points. We finally conclude in section IV. Some
important formulae can be found in the Appendix V.

##  II The sub-GeV dark matter framework and $3\to 2$ dynamics

We extend the SM with the real gauge singlet scalars $\phi$ and $\delta$
Mohanty:2019drv . A $\mathbb{Z}_{2}$ symmetry is imposed under which
$\phi\to-\phi$ while all other fields transform trivially. This negative
$\mathbb{Z}_{2}$ charge for $\phi$ stabilises the same and makes it a
potential DM candidate. The scalar potential of the setup reads

$\displaystyle V$ $\displaystyle=$
$\displaystyle\frac{\mu_{\phi}^{2}}{2}\phi^{2}+\frac{\mu_{\delta}^{2}}{2}\delta^{2}+\mu_{H}^{2}H^{\dagger}H+\frac{\mu_{21}}{2}\phi^{2}\delta+\frac{\mu_{22}}{3!}\delta^{3}$
(1)
$\displaystyle+\mu_{23}H^{\dagger}H\delta+\frac{\lambda_{11}}{4!}\phi^{4}+\frac{\lambda_{12}}{4}\phi^{2}\delta^{2}+\frac{\lambda_{22}}{4!}\delta^{4}$
$\displaystyle+\frac{\lambda_{13}}{2}\phi^{2}H^{\dagger}H+\frac{\lambda_{23}}{2}\delta^{2}H^{\dagger}H+\lambda_{33}\left[H^{\dagger}H\right]^{2}.$

In the above, all parameteres are taken real and $H$ denotes the SM Higgs
doublet. The doublet $H$ and the singlet $\delta$ receive vacuum expectation
values (VEVs) $v$ and $v_{\delta}$ respectively. One then writes

$\displaystyle H=\begin{pmatrix}G^{+}\\\
\frac{1}{\sqrt{2}}(v+h+iG^{0})\end{pmatrix};~{}~{}~{}\delta=v_{\delta}+\delta^{\prime}.$
(2)

Here, $v$ = 246 GeV. It is noted that $\phi$ does not receive a VEV and thus
the $\mathbb{Z}_{2}$ remains intact. In all, such a configuration of the VEVs
introduces a mixing between $h$ and $\delta^{\prime}$. The mixing of $\phi$ is
however forbidden by the $\mathbb{Z}_{2}$ symmetry. The mass matrix in the
$(h~{}~{}\delta^{\prime})$ basis is

$\displaystyle\mathcal{M}^{2}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}2\lambda_{33}v^{2}&\lambda_{23}\hskip
1.42262ptvv_{\delta}+\mu_{23}v\\\ \lambda_{23}\hskip
1.42262ptvv_{\delta}\hskip
1.42262pt+\mu_{23}v&\frac{\mu_{22}v_{\delta}}{2}+\frac{\lambda_{22}}{3}v_{\delta}^{2}-\frac{\mu_{23}}{2}\frac{v^{2}}{v_{\delta}}\end{pmatrix}.$
(3)

We have made use of the tadpole conditions $\frac{\partial V}{\partial
h}=0=\frac{\partial V}{\partial\delta^{\prime}}$ while deriving Eq.(3). One
derives the tadpole conditions to be

$\displaystyle\mu_{h}^{2}$ $\displaystyle=$
$\displaystyle-\big{(}\mu_{23}\hskip
2.84526ptv_{\delta}+\frac{\lambda_{23}}{2}\hskip
2.84526ptv_{\delta}^{2}+\lambda_{33}\hskip 2.84526ptv^{2}\big{)},$ (4a)
$\displaystyle\mu_{\delta}^{2}$ $\displaystyle=$
$\displaystyle-\big{(}\frac{\mu_{22}}{2}v_{\delta}+\frac{\lambda_{22}}{6}v_{\delta}^{2}+\frac{\mu_{23}}{2}\frac{v^{2}}{v_{\delta}}+\frac{\lambda_{23}}{2}v^{2}\big{)}.$
(4b)

Eq.(3) is diagonalised through a $2\times 2$ rotation parameterised by a
mixing angle $\theta$. That is,

$\displaystyle\begin{pmatrix}h\\\
\delta^{\prime}\end{pmatrix}=\begin{pmatrix}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{pmatrix}\begin{pmatrix}h_{1}\\\
h_{2}\end{pmatrix},$ (5)

which leads to the eigenstates $(h_{1},h_{2})$. We identify $h_{1}$ with the
discovered Higgs boson with $m_{h_{1}}$ = 125 GeV. The DM mass is given by

$\displaystyle
m_{\phi}^{2}=\mu^{2}_{\phi}+\mu_{21}v_{\delta}+\frac{1}{2}\lambda_{12}v^{2}_{\delta}+\frac{1}{2}\lambda_{13}v^{2}.$
(6)

Some of the parameters in the scalar potential can be expressed in terms of
the physical quantities of the theory such as masses and mixing angles using
the relations

$\displaystyle\lambda_{33}$ $\displaystyle=$
$\displaystyle\frac{(m_{h_{1}}^{2}\cos^{2}\theta+m_{h_{2}}^{2}\sin^{2}\theta)}{2\hskip
1.42262ptv^{2}},$ (7a) $\displaystyle\mu_{23}$ $\displaystyle=$
$\displaystyle\frac{(m_{h_{2}}^{2}-m_{h_{1}}^{2})\cos\theta\sin\theta}{v}-\lambda_{23}\hskip
1.42262ptv_{\delta},$ (7b) $\displaystyle\mu_{22}$ $\displaystyle=$
$\displaystyle\frac{2}{v_{\delta}}(m_{h_{1}}^{2}\sin^{2}\theta+m_{h_{2}}^{2}\cos^{2}\theta)-\frac{\lambda_{22}}{3}\hskip
1.42262ptv_{\delta}$ (7c)
$\displaystyle+\mu_{23}\frac{v^{2}}{v^{2}_{\delta}},$
$\displaystyle\mu_{\phi}^{2}$ $\displaystyle=$ $\displaystyle
m^{2}_{\phi}-\mu_{21}v_{\delta}-\frac{1}{2}\lambda_{12}v^{2}_{\delta}-\frac{1}{2}\lambda_{13}v^{2}.$
(7d)

The independent variables of the present framework thus are
$\\{m_{h_{2}},m_{\phi},\mu_{21},v_{h},v_{\delta},\text{sin}\theta,\lambda_{11},\lambda_{12},\lambda_{13},\lambda_{22},\lambda_{23}\\}$.

The relic abundance in a strongly interacting DM setup is generated through
number changing processes. In this setup, the leading number changing
processes consistent with the $\mathbb{Z}_{2}$ is $\phi\phi\phi\to\phi S$
where $S=h_{1},h_{2}$. We list below the Feynman diagrams for illustration.

Figure 1: $\phi$-mediated Feynman diagram responsible for annihilation of
$\phi$.

Figure 2:  Feynman diagrams show annihilation of $\phi$ mediated by $h_{1}$,
$h_{2}$.

Figure 3:  Annihilation diagrams mediated by $h_{1}$, $h_{2}$ and $\phi$.

Figs.1, 2 and 3 show the $3\to 2$ annihilation diagrams. We take the limit
where sin$\theta$ and $\lambda_{13}$ are tiny. The tiny sin$\theta$ conforms
with the constraints on $\delta-H$ mixing from direct search and Higgs signal
strengths. The $\phi-\phi-h_{1}$ trilinear coupling becomes negligibly small
in this limit while the $\phi-\phi-h_{2}$ coupling becomes
$\simeq(\mu_{21}+\lambda_{12}v_{\delta})$. Moreover, the limit entails that
the $\phi$-mediated amplitude in Fig.1 is the dominant one. The following is
its expression:

$\displaystyle\mathcal{M}\simeq\frac{\lambda_{11}(\mu_{21}+\lambda_{12}v_{\delta})}{(p_{1}+p_{2}+p_{3})^{2}-m^{2}_{\phi}}\simeq\frac{\lambda_{11}(\mu_{21}+\lambda_{12}v_{\delta})}{8m_{\phi}^{2}}.$
(8)

In the last step above, we have neglected the DM velocity, a justifiable
assumption for non-relativistic DM. One can approximate
$p_{i}\simeq(m_{\phi},\vec{0})$ in that case. The evolution of the DM comoving
number density $Y_{\phi}$ is described by the Boltzmann equation Kolb:1990vq

$\displaystyle\frac{dY_{\phi}}{dx}=-\frac{x~{}s^{2}}{3\mathcal{H}}\langle\sigma
v^{2}\rangle_{\phi\phi\phi\rightarrow\phi
h_{2}}\Big{(}Y^{3}_{\phi}-Y_{\phi}Y^{\text{eq}}_{\phi}Y^{\text{eq}}_{h_{2}}\Big{)},$
(9)

where $s=\frac{2\pi^{2}}{45}~{}g_{*s}T^{3}$ is the entropy density and the
Hubble parameter is
$\mathcal{H}=\sqrt{\frac{8\pi^{3}Gg_{*}T^{4}}{90}}$111$g_{*}$ is the effective
relativistic degrees of freedom.. Besides, $Y^{\text{eq}}_{P}$ denotes the
equilibrium comoving number density of a particle $P$ and $\langle\sigma
v^{2}\rangle_{\phi\phi\phi\rightarrow\phi h_{2}}$ is the thermally averaged
$3\to 2$ cross section Pierre:2018man . It can be related to $\sigma
v^{2}_{\phi\phi\phi\rightarrow\phi h_{2}}$ using modified Bessel functions
$K_{n}(x)$ as

$\displaystyle\braket{\sigma v^{2}}_{\phi\phi\phi\to\phi h_{2}}(T)$
$\displaystyle=$
$\displaystyle\Bigg{(}\frac{K_{1}\big{(}\frac{m_{\phi}}{T}\big{)}}{K_{2}\big{(}\frac{m_{\phi}}{T}\big{)}}\Bigg{)}^{3}(\sigma
v^{2})_{\phi\phi\phi\to\phi h_{2}},\text{with}$

$\displaystyle(\sigma v^{2})_{\phi\phi\phi\to\phi h_{2}}$ $\displaystyle=$
$\displaystyle\frac{|\mathcal{M}|^{2}}{72\pi m^{3}_{\phi}}$
$\displaystyle\times$
$\displaystyle\sqrt{\bigg{(}1-\frac{5}{16}\Big{(}\frac{m_{h_{2}}}{m_{\phi}}\Big{)}^{2}+\frac{1}{64}\Big{(}\frac{m_{h_{2}}}{m_{\phi}}\Big{)}^{4}\bigg{)}}.$
(10)

Further, Eq.(9) becomes the following upon taking $g_{*s}\simeq g_{*}$.

$\displaystyle\frac{dY_{\phi}}{dx}$ $\displaystyle=$
$\displaystyle-0.116~{}g^{3/2}_{*}~{}M_{\text{Pl}}\frac{m^{4}_{\phi}}{x^{5}}~{}\langle\sigma
v^{2}\rangle_{\phi\phi\phi\rightarrow\phi h_{2}}$ (11)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times\Big{(}Y^{3}_{\phi}-Y_{\phi}Y^{\text{eq}}_{\phi}Y^{\text{eq}}_{h_{2}}\Big{)}.$

We present an approximate analytical solution to Eq.(11) to gain insight on
the dynamics. We follow here the approach adopted in Bhattacharya:2019mmy . We
take $B\equiv
0.116~{}g^{3/2}_{*}~{}M_{\text{Pl}}~{}m^{4}_{\phi}~{}\langle\sigma
v^{2}\rangle_{\phi\phi\phi\rightarrow\phi h_{2}}$ and rewrite Eq.(11) using
$\Delta=Y_{\phi}-Y^{\text{eq}}_{\phi}$ as

$\displaystyle\frac{d(\Delta+Y^{\text{eq}}_{\phi})}{dx}$ $\displaystyle=$
$\displaystyle-\frac{B}{x^{5}}\Big{[}(\Delta+Y^{\text{eq}}_{\phi})^{3}$ (12)
$\displaystyle-$
$\displaystyle(\Delta+Y^{\text{eq}}_{\phi})~{}Y^{\text{eq}}_{\phi}~{}Y^{\text{eq}}_{h_{2}}\Big{]}.$

Before freeze-out, i.e. for $1\ll x\leq x_{f}$ ($x_{f}$ denotes freeze out of
DM), $\Delta\ll Y^{eq}_{\phi}$ and $\frac{d\Delta}{dx}\rightarrow 0$ and near
freeze-out, i.e. for $x\simeq x_{f}$, so
$\Delta(x_{f})=cY^{\text{eq}}_{\phi}(x_{f})$. One then writes

$\displaystyle\frac{dY^{\text{eq}}_{\phi}(x_{f})}{dx}$ $\displaystyle=$
$\displaystyle-\frac{B}{x^{5}}\Big{[}(c+1)^{3}(Y_{\phi}^{\text{eq}})^{3}(x_{f})$
(13) $\displaystyle-$
$\displaystyle(c+1)~{}(Y_{\phi}^{\text{eq}})^{2}(x_{f})(Y_{h_{2}}^{\text{eq}})(x_{f})\Big{]}.$

Using the equilibrium distributions

$\displaystyle Y^{\text{eq}}_{\phi}(x)$ $\displaystyle=$ $\displaystyle
0.145~{}\Big{(}\frac{g_{\phi}}{g_{*}}\Big{)}~{}x^{3/2}~{}e^{-x},$
$\displaystyle Y^{\text{eq}}_{h_{2}}(x)$ $\displaystyle=$ $\displaystyle
0.145~{}\Big{(}\frac{g_{h_{2}}}{g_{*}}\Big{)}~{}\Big{(}\frac{m_{h_{2}}}{m_{\phi}}\Big{)}~{}x^{3/2}~{}e^{-\frac{m_{h_{2}}}{m_{\phi}}x},$
(14)

Eq.(13) takes the form

$\displaystyle x^{2}_{f}-\frac{3}{2}x_{f}$ $\displaystyle=$
$\displaystyle(0.145)^{2}B(c+1)\Big{(}\frac{g_{\phi}}{g_{*}}\Big{)}^{2}e^{-2x_{f}}$
(15) $\displaystyle\times$
$\displaystyle\Big{[}(c+1)^{2}-\Big{(}\frac{m_{h_{2}}}{m_{\phi}}\Big{)}^{3/2}e^{-\Big{(}\frac{m_{h_{2}}}{m_{\phi}}+1\Big{)}x_{f}}\Big{]}.$

Eq.15 can be iteratively solved for $x_{f}$. Now
$Y^{\text{eq}}_{\phi}\ll\Delta$ for $x\gg x_{f}$. Eq.12 can then be integrated
from the freeze-out epoch $x=x_{f}$ to a later epoch $x\to\infty$ as

$\displaystyle\int_{\Delta(x_{f})}^{\Delta(x\to\infty)}\frac{d\Delta}{\Delta^{3}}=-B\int_{x_{f}}^{\infty}dx~{}\frac{1}{x^{5}}$
(16)

to finally yield
$Y_{\phi}(x\to\infty)=x_{f}^{2}\bigg{[}\frac{2}{0.12M_{\text{Pl}}m^{4}_{\phi}\braket{\sigma
v^{2}}}\bigg{]}^{1/2}$. The final relic abundance is related to the comoving
density $Y_{\phi}(x\to\infty)$ as $\Omega h^{2}=2.752\times
10^{8}~{}m_{\phi}Y_{\phi}(x\to\infty)$ Kolb:1990vq ; Bhattacharya:2016ysw .

##  III First order phase transitions and numerical analysis

We are interested in cosmological phase transitions in the direction of the
field $\delta$ given its crucial role in the generation of the DM relic
through the $3\to 2$ transitions. The tree-level potential in the
$\delta$-direction looks like
$V_{0}(\delta_{\text{cl}})=\frac{\mu_{\delta}^{2}}{2}\delta_{\text{cl}}^{2}+\frac{\mu_{22}}{3!}\delta_{\text{cl}}^{3}+\frac{\lambda_{22}}{4!}\delta_{\text{cl}}^{4}$.
Here $\delta_{\text{cl}}$ denotes the classical rolling field and is distinct
from the degree of freedom $\delta^{\prime}$. Next, the one-loop Coleman-
Weinberg correction to the tree-level potential reads

$\displaystyle\Delta V_{\rm CW}(\delta_{\text{cl}})$ $\displaystyle=$
$\displaystyle\frac{1}{64\pi^{2}}\sum_{a=\phi,\delta^{\prime},h,G^{+},G^{0}}n_{a}\bigg{[}m^{4}_{a}(\delta_{\text{cl}})\Big{(}\text{log}\frac{m_{i}^{2}(\delta_{\text{cl}})}{m_{a}^{2}(v_{\delta})}-\frac{3}{2}\Big{)}$
(17) $\displaystyle\hskip
79.6678pt+2m^{2}_{a}(v_{\delta})m^{2}_{a}(\delta_{\text{cl}})\bigg{]}.$

The sum in Eq.(17) runs over the scalars of the theory since these are the
fields that develop $\delta_{\text{cl}}$-dependent masses. Also, $n_{a}$
denotes number of degrees of freedom of the $a$th scalar. One finds
$n_{\phi}=n_{\delta}=n_{h}=n_{G^{0}}=1$ while $n_{G^{+}}=2$. The field-
dependent masses are

$\displaystyle m^{2}_{\phi}(\delta_{\text{cl}})$ $\displaystyle=$
$\displaystyle\mu_{\phi}^{2}+\frac{\lambda_{12}\hskip
1.42262pt\delta_{\text{cl}}^{2}}{2}+\mu_{21}\hskip
1.42262pt\delta_{\text{cl}},$ (18a) $\displaystyle
m^{2}_{\delta}(\delta_{\text{cl}})$ $\displaystyle=$
$\displaystyle\mu_{\delta}^{2}+\frac{\lambda_{22}\hskip
1.42262pt\delta_{\text{cl}}^{2}}{2}+\mu_{22}\hskip
1.42262pt\delta_{\text{cl}},$ (18b) $\displaystyle
m^{2}_{h}(\delta_{\text{cl}})$ $\displaystyle=$ $\displaystyle
m^{2}_{G^{0}}(\delta_{\text{cl}})=m^{2}_{G^{+}}(\delta_{\text{cl}}),$ (18c)
$\displaystyle=$ $\displaystyle\mu_{h}^{2}+\frac{\lambda_{23}\hskip
1.42262pt\delta_{\text{cl}}^{2}}{2}+\mu_{23}\hskip
1.42262pt\delta_{\text{cl}}.$

The one-loop correction to the scalar potential induced in presence of $T\neq
0$ reads

$\displaystyle\Delta V_{1}(\delta_{\text{cl}},T)$ $\displaystyle=$
$\displaystyle\frac{T^{4}}{2\pi^{2}}\sum_{a=\phi,\delta^{\prime},h,G^{+},G^{0}}n_{a}J_{B}[m^{2}_{a}(\delta_{\text{cl}})/T^{2}].$
(19)

Here, the function $J_{B}(\frac{m^{2}}{T^{2}})$ is defined as
$J_{B}(\frac{m^{2}}{T^{2}})=\int_{0}^{\infty}dx~{}x^{2}~{}\text{log}\big{[}1-e^{-\sqrt{x^{2}+\frac{m^{2}}{T^{2}}}}\big{]}$.
It can be approximated in the high temperature limit as

$\displaystyle J_{B}\Big{(}\frac{m^{2}}{T^{2}}\Big{)}$ $\displaystyle\simeq$
$\displaystyle-\frac{\pi^{4}}{45}+\frac{\pi^{2}}{12}~{}\frac{m^{2}}{T^{2}}-\frac{\pi}{6}~{}\frac{m^{3}}{T^{3}}$
(20)
$\displaystyle-\frac{1}{32}~{}\frac{m^{4}}{T^{4}}~{}\text{ln}\Big{(}\frac{m^{2}}{a_{B}T^{2}}\Big{)},$

where $a_{B}\equiv 16\pi^{2}~{}e^{\frac{3}{2}-2\gamma_{E}}$. Further, the
infrared effects are included using the daisy remmuation technique Weinberg ;
Dolan ; Kirzhnits:1976ts ; Carrington:1991hz ; Gross ; Fendley:1987ef ;
Kapusta ; Arnold:1992rz ; Parwani:1991gq . More precisely, the Arnold-Espinosa
prescription Arnold:1992rz in adopted in this study. That is

$\displaystyle\Delta V_{\text{ring}}(\delta_{\text{cl}},T)$ $\displaystyle=$
$\displaystyle-\frac{T}{12\pi}\sum_{\begin{subarray}{c}a=\phi,\delta^{\prime},\\\
h,G^{+},G^{0}\end{subarray}}\Big{[}(m_{a}^{2}(\delta_{\text{cl}},T))^{3/2}-(m_{a}^{2}(\delta_{\text{cl}}))^{3/2}\Big{]}.$

The total thermal scalar potential as a function of $\delta_{cl}$ and $T$ is
the sum of the individual components. That is

$\displaystyle V_{\text{total}}(\delta_{\text{cl}},T)$ $\displaystyle=$
$\displaystyle V_{0}(\delta_{\text{cl}})+\Delta V_{\rm
CW}(\delta_{\text{cl}})$ (22) $\displaystyle+$ $\displaystyle\Delta
V_{1}(\delta_{\text{cl}},T)+\Delta V_{\text{ring}}(\delta_{\text{cl}},T).$

We attempt to understand $V_{\text{total}}(\delta_{\text{cl}},T)$ more
intuitively with the aid of appropriate approximations. The thermal potential
can be cast as a polynomial in $\delta_{\text{cl}}$ using the high temperature
approximation for $J_{B}(\frac{m^{2}}{T^{2}})$ given in Eq.(20) and discarding
terms that do not contain $\delta_{\text{cl}}$. One then straightforwardly
derives

$\displaystyle V_{\text{total}}(\delta_{\text{cl}},T)$ $\displaystyle\simeq$
$\displaystyle\Big{[}\frac{\mu_{\delta}^{2}}{2}+\frac{1}{48}(\lambda_{12}+\lambda_{22}+\lambda_{23})~{}T^{2}\Big{]}~{}\delta_{\text{cl}}^{2}$
(23) $\displaystyle-$
$\displaystyle\frac{1}{24\sqrt{2}\pi}\Big{[}\Big{(}\lambda^{3/2}_{12}+\lambda^{3/2}_{22}+\lambda^{3/2}_{23}\Big{)}~{}T-\frac{\mu_{22}}{6}\Big{]}~{}\delta_{\text{cl}}^{3}$
$\displaystyle+$
$\displaystyle\frac{\overline{\lambda}(T)}{4}\delta_{\text{cl}}^{4}.$

where,
$\overline{\lambda}(T)=\frac{\lambda_{22}}{6}+\frac{1}{64\pi^{2}}\Big{[}\lambda^{2}_{12}~{}\text{ln}\Big{(}\frac{a_{B}~{}T^{2}}{e^{3/2}~{}m^{2}_{\phi}(v_{\delta})}\Big{)}+\lambda^{2}_{22}~{}\text{ln}\Big{(}\frac{a_{B}~{}T^{2}}{e^{3/2}~{}m^{2}_{\delta}(v_{\delta})}\Big{)}+\lambda^{2}_{23}~{}\text{ln}\Big{(}\frac{a_{B}~{}T^{2}}{e^{3/2}~{}m^{2}_{h}(v_{\delta})}\Big{)}\Big{]}$.
The form in Eq.(III) resembles the SM formula in Quiros:1999jp and therefore
permits an analytical treatment of the thermal evolution. We identify certain
temperature thresholds in Eq.(III). An extrema at $\delta_{\text{cl}}=0$ is
readily identified. Moreover, this extrema is a minima (maxima) for $T>T_{0}$
($T<T_{0}$) with
$T_{0}=\sqrt{-24\mu^{2}_{\delta}/(\lambda_{12}+\lambda_{22}+\lambda_{23})}$.
An inflection point is present for $T=T_{1}$ is encountered. The temperature
$T_{1}$ can be solved from
$\Big{[}\Big{(}\lambda^{3/2}_{12}+\lambda^{3/2}_{22}+\lambda^{3/2}_{23}\Big{)}~{}T_{1}-\frac{\mu_{22}}{6}\Big{]}^{2}=1152\pi^{2}\overline{\lambda}(T_{1})\Big{[}\mu_{\delta}^{2}+\frac{1}{24}(\lambda_{12}+\lambda_{22}+\lambda_{23})~{}T_{1}^{2}\Big{]}$.
For $T_{0}<T<T_{1}$, a maxima and a minima appear at
$\delta_{\text{cl}}=\delta_{\text{max}}(T),\delta_{\text{min}}(T)$ in addition
to the minima at $\delta_{\text{cl}}=0$. Thus, the $T_{0}<T<T_{1}$ temperature
band is the most interesting from the perspective of first order phase
transitions on account of the coexisting minima. The critical temperature
$T_{c}$ for this case is the temperature at which the minima at
$\delta_{\text{cl}}=0,\delta_{\text{min}}$ are degenerate, i.e.,
$V_{\text{total}}(0,T_{c})=V_{\text{total}}(\delta_{\text{min}}(T_{c}),T_{c})$.
A strong first order phase transition is identified by
$\frac{\delta_{\text{min}}(T_{c})}{T_{c}}\geq 1$.

We put forth two benchmark points in Table 1. The independent model parameters
and the corresponding critical temperatures are listed therein. These
benchmarks also predict relic densities within the Planck band. The shape of
$V_{\text{total}}(\delta_{\text{cl}},T)$ around the crtical temperatures is
shown in Fig.4.

Figure 4: The extrema of the thermal potential around $T=T_{c}$ for BP1 and
BP2. The color coding is shown in the legends.

Both BP1 and BP2 correspond to sub-GeV critical temperatures. One also
inspects in Fig.4 that $\delta(T_{c})\gtrsim$ 1 GeV for both thereby
signalling a SFOPT.

| $\lambda_{11}$ | $\lambda_{12}$ | $\lambda_{22}$ | $\lambda_{23}$ |  $\mu_{22}$ (GeV) |  $m_{h_{2}}$ (GeV) |  $m_{\phi}$ (GeV) |  $v_{\delta}$ (GeV) |  $\mu_{21}$ (GeV) |  $T_{c}$ (GeV) | $\frac{\delta(T_{c})}{T_{c}}$ | $\Omega_{\text{D}M}h^{2}$
---|---|---|---|---|---|---|---|---|---|---|---|---
BP1 | 0.1 | 4.69202 | 0.0060999 | 0.00001 | -0.208174 | 0.0647581 | 0.050095 | 66.8 | 0.756 | 0.975268 | 1.47608 | 0.114796
BP2 | 0.5 | 5.39401 | 0.0060999 | 0.00001 | -0.144512 | 0.116026 | 0.0830917 | 45.3 | 0.586 | 0.717489 | 1.46328 | 0.114441

Table 1: Benchmark model parameters along with the predicted relic densities
and the crtical temperatures.

|  $T_{n}$ (GeV) | $\alpha$ | $\beta/H$ |  $f_{\text{col}}$ (Hz) |  $f_{\text{sw}}$ (Hz) |  $f_{\text{tur}}$ (Hz)
---|---|---|---|---|---|---
BP1 | 0.821609 | 0.00931039 | 1119.22 |  4.16229 $\times 10^{-5}$ |  2.71898 $\times 10^{-4}$ |  3.86382 $\times 10^{-4}$
BP2 | 0.606572 | 0.009944 | 1149.43 |  3.15148 $\times 10^{-5}$ |  2.0544 $\times 10^{-4}$ |  2.91941 $\times 10^{-4}$

Table 2: The nucleation temperature and other GW parameters corresponding to
BP1 and BP2.

A combined analysis of the relic density and SFOPT is in order. We first fix
$\lambda_{13}=\lambda_{23}=10^{-5},~{}\lambda_{11}=0.1,~{}0.5,~{}1$ and vary
the rest of the parameters as
$|\sin\theta|<0.001,~{}|\lambda_{12}|<6,~{}|\lambda_{22}|<6,~{}|\mu_{21}|<1~{}\text{GeV},~{}1~{}\text{MeV}<m_{\phi}<1~{}\text{GeV},~{}1~{}\text{MeV}<m_{h_{2}}<2m_{\phi}$.
It can be noted that the small values taken for sin$\theta$ and $\lambda_{13}$
render the $\phi-\phi-h_{1(2)}$ interactions appropriately small so as to
conform with possible stringent direct detection bounds in the near future.
The small sin$\theta$ also implies safety from the Higgs signal strength
constraints and the direct search constraints involving $h_{2}$. The full
thermal potential is analysed using the publicly available tool PhaseTracer
Athron:2020sbe . We select the parameter points that lead to (a)
$0.1118<\Omega h^{2}<0.1280$ Planck:2018vyg and (b)
$\frac{\delta_{\text{min}}(T_{c})}{T_{c}}>1$. The results are shown as plots
in the $m_{\phi}-\lambda_{12}$ plane in Fig.5.

Figure 5: Parameter points in the $m_{\phi}-\lambda_{12}$ plane allowed by
$\Omega h^{2}$ and SFOPT criteria for $\lambda_{11}$ = 0.1 (top left), 0.5
(top right) and 1 (bottom). The color coding is described in the legends.

Fig.5 shows a $\mathcal{O}$(10)-$\mathcal{O}$(100)(MeV) mass range for the DM
thereby corroborating the results of Mohanty:2019drv . The correlation among
$m_{\phi}$, $\lambda_{11}$ and $\lambda_{12}$ can be understood from $(\sigma
v^{2})_{\phi\phi\phi\to\phi
h_{2}}\sim\frac{\lambda^{2}_{11}(\mu_{21}+\lambda_{12}v_{\delta})^{2}}{m^{7}_{\phi}}$
in Eq.(10). The maximum allowed value of $m_{\phi}$ increases with increasing
$\lambda_{11}$. And this is concurred by Fig.5 where one inspects that for the
aforesaid variation of $\lambda_{12}$, changing $\lambda_{11}=0.1$ to
$\lambda_{11}=1$ takes one from $m_{\phi}\lesssim$ 70 MeV to
$m_{\phi}\lesssim$ 190 MeV. Demanding a SFOPT restricts the parameter space
further.

Figure 6: The nucleation temperature and other GW parameters corresponding to
BP1 and BP2.

We now discuss the GW spectrum of such a scenario and a few pertinent
quantities are therefore introduced. The parameter $\beta$ is defined as

$\displaystyle\frac{\beta}{H}$ $\displaystyle=$
$\displaystyle\bigg{[}T\frac{d}{dT}\Big{(}\frac{S_{3}}{T}\Big{)}\bigg{]}_{T_{n}}.$
(24)

Here, $T_{n}$ is the nucleation temperature and $S_{3}$ is the Euclidean
action in three dimensions. Next, we define $\Delta V_{\text{tot}}(T)\equiv
V_{\text{tot}}(0,T)-V_{\text{tot}}(\delta_{\text{min}}(T),T)$ which measures
the difference in the depths of the potential at the two vacua at a
temperature $T$. The energy budget of the phase transition during the bubble
nucleation is then given by

$\displaystyle\epsilon$ $\displaystyle=$ $\displaystyle\Delta
V_{\text{tot}}(T_{n})-\bigg{[}T\frac{d\Delta
V_{\text{tot}}(T)}{dT}\bigg{]}_{T_{n}}.$ (25)

One subsequently defines $\alpha=\frac{\epsilon}{\rho_{n}}$, where
$\rho_{n}=\frac{g_{*}\pi^{2}}{30}T_{n}^{4}$, is the energy density during
nucleation. The sources of GWs associated with SFOPT are bubble collision
Kosowsky ; Turner ; Huber_2008 ; Watkins ; Marc ; Caprini_2008 , sound waves
Hindmarsh ; Leitao:2012tx ; Giblin:2013kea ; Giblin:2014qia ; Hindmarsh_2015 ,
and, magnetohydrodynamic turbulence Chiara ; Kahniashvili ;
Kahniashvili:2008pe ; Kahniashvili:2009mf ; Caprini:2009yp . The GW amplitudes
from these sources are functions of the frequency. These amplitudes attain
their maximum values at certain “peak frequencies” that we denote as
$f_{\text{col}},f_{\text{sw}}$ and $f_{\text{tur}}$. The expressions for these
frequencies and the corresponding GW amplitudes is given in the Appendix V.

The $\alpha$ and $\beta$ parameters are crucial in deciding the strength of
the GW signals stemming from a SFOPT. Table 2 displays $T_{n},\alpha,\beta$
for the chosen benchmarks. Bubble nucleation at the sub-GeV temperatures shown
can be understood as an artefact of the sub-GeV DM masses involved. This is an
important observation of our analysis. We also inspect in table 2 that
$\alpha$ and $\beta/H$ at nucleation are respectively $\mathcal{O}(10^{-3})$
and $\mathcal{O}(10^{3})$ for the BPs. Such ball-park values entail that the
magnitude of the total GW amplitude peaks at $\mathcal{O}(10^{-18})$, as is
seen in Fig.6. And the corresponding frequencies are
$\sim\mathcal{O}(10^{-4})$ Hz. This is expected since the total GW amplitude
is dominated by the sound wave component that has
$f_{\text{sw}}\sim\mathcal{O}(10^{-4})$ Hz. Overall, an
$\Omega_{\text{GW}}h^{2}\sim\mathcal{O}(10^{-18})$ with peak at
$f\sim\mathcal{O}(10^{-4})$ Hz is beyond the reach of the proposed GW
detectors. However, u-DECIGO Kudoh:2005as is expected to detect GW signals of
a similar magnitude, albeit at a different frequency range. In all, a future
detector with higher sensitivity in the $\sim\mathcal{O}$(0.1) mHz regime can
potentially detect GW signals of a sub-GeV DM species exhibiting number
changing annihilations.

##  IV Conclusions

In this study, we look at possible GW signals originating from a SIMP. The
SIMP species freezes-out via $3\to 2$ number changing processes that are in
turn driven by another scalar $\delta$. The observed relic density is attained
for sub-GeV DM masses. We have studied the shape of the thermal potential
along $\delta$ and identified the region of the parameter space permitting a
SFOPT. For a few representative benchmarks, we have then computed the GW
amplitude as a function of frequency. Our analysis reveals that a futuristic
detector with a sensitivity $\Omega_{\text{GW}}h^{2}\sim 10^{-18}$ in the mHz
regime can detect GWs stemming from such a SIMP framework.

##  V Appendix

The peak frequencies corresponding to bubble collisions, sound waves and
turbulence are expressed as

$\displaystyle f_{\text{col}}$ $\displaystyle=$ $\displaystyle 1.65\times
10^{-5}\Big{(}\frac{0.62}{1.8-0.1v_{b}+v^{2}_{b}}\Big{)}$ (26a)
$\displaystyle\frac{\beta}{H}\Big{(}\frac{T_{n}}{100}\Big{)}\Big{(}\frac{g_{*}}{100}\Big{)}^{\frac{1}{6}},$
$\displaystyle f_{\text{sw}}$ $\displaystyle=$ $\displaystyle 1.9\times
10^{-5}\Big{(}\frac{1}{v_{b}}\Big{)}\frac{\beta}{H}\Big{(}\frac{T_{n}}{100}\Big{)}\Big{(}\frac{g_{*}}{100}\Big{)}^{\frac{1}{6}},$
(26b) $\displaystyle f_{\text{tur}}$ $\displaystyle=$ $\displaystyle 2.7\times
10^{-5}\Big{(}\frac{1}{v_{b}}\Big{)}\frac{\beta}{H}\Big{(}\frac{T_{n}}{100}\Big{)}\Big{(}\frac{g_{*}}{100}\Big{)}^{\frac{1}{6}}.$
(26c)

We express below the GW amplitudes from the three sources as a function of
frequency $f$.

$\displaystyle\Omega_{\text{coll}}(f)$ $\displaystyle=$ $\displaystyle
1.67\times
10^{-5}\Big{(}\frac{\beta}{H}\Big{)}^{-2}\Big{(}\frac{0.11v^{3}_{b}}{0.42+v^{2}_{b}}\Big{)}\Big{(}\frac{\kappa\alpha}{1+\alpha}\Big{)}^{2}$
(27)
$\displaystyle\Big{(}\frac{g_{*}}{100}\Big{)}^{-\frac{1}{3}}\Big{(}\frac{3.8(f/f_{\text{coll}})^{2.8}}{1+2.8(f/f_{\text{coll}})^{3.8}}\Big{)}$
$\displaystyle\Omega_{\text{sw}}(f)$ $\displaystyle=$ $\displaystyle
2.65\times
10^{-6}\Big{(}\frac{\beta}{H}\Big{)}^{-2}v_{b}\Big{(}\frac{\kappa_{v}\alpha}{1+\alpha}\Big{)}^{2}$
(28)
$\displaystyle\Big{(}\frac{g_{*}}{100}\Big{)}^{-\frac{1}{3}}\Big{(}\frac{f}{f_{\text{sw}}}\Big{)}^{3}\Big{(}\frac{7}{4+3(f/f_{\text{sw}})^{2}}\Big{)}^{2},$
$\displaystyle\Omega_{\text{tur}}(f)$ $\displaystyle=$ $\displaystyle
3.35\times
10^{-4}\Big{(}\frac{\beta}{H}\Big{)}^{-2}v_{b}\Big{(}\frac{\epsilon\kappa_{v}\alpha}{1+\alpha}\Big{)}^{1.5}$
(29)
$\displaystyle\Big{(}\frac{g_{*}}{100}\Big{)}^{-\frac{1}{3}}\Big{(}\frac{f}{f_{\text{tur}}}\Big{)}^{3}\frac{(1+f/f_{\text{tur}})^{-11/3}}{1+8\pi
f/h_{s}}.$

##  VI Acknowledgements

NC acknowledges support from DST, India, under grant number IFA18-PH214
(INSPIRE Faculty Award). HR is supported by the Science and Engineering
Research Board, Government of India, under the agreement SERB/PHY/2016348
(Early Career Research Award). TS acknowledges the support from the Dr. D. S.
Kothari Postdoctoral scheme No. PH/20-21/0163.

## References

  * (1) V. C. Rubin and W. K. Ford, Jr., _Rotation of the Andromeda Nebula from a Spectroscopic Survey of Emission Regions_ , _Astrophys. J._ 159 (1970) 379–403.
  * (2) G. Bertone, D. Hooper and J. Silk, _Particle dark matter: Evidence, candidates and constraints_ , _Phys. Rept._ 405 (2005) 279–390, [hep-ph/0404175].
  * (3) W. Hu and S. Dodelson, _Cosmic Microwave Background Anisotropies_ , _Ann. Rev. Astron. Astrophys._ 40 (2002) 171–216, [astro-ph/0110414].
  * (4) WMAP collaboration, G. Hinshaw et al., _Nine-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Cosmological Parameter Results_ , _Astrophys. J. Suppl._ 208 (2013) 19, [1212.5226].
  * (5) Planck collaboration, N. Aghanim et al., _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6, [1807.06209].
  * (6) G. R. Blumenthal, S. M. Faber, J. R. Primack and M. J. Rees, _Formation of Galaxies and Large Scale Structure with Cold Dark Matter_ , _Nature_ 311 (1984) 517–525.
  * (7) E. W. Kolb and M. S. Turner, _The Early Universe_ , vol. 69.  1990, 10.1201/9780429492860.
  * (8) G. Jungman, M. Kamionkowski and K. Griest, _Supersymmetric dark matter_ , _Phys. Rept._ 267 (1996) 195–373, [hep-ph/9506380].
  * (9) J. L. Feng, _Dark Matter Candidates from Particle Physics and Methods of Detection_ , _Ann. Rev. Astron. Astrophys._ 48 (2010) 495–545, [1003.0904].
  * (10) G. Arcadi, M. Dutra, P. Ghosh, M. Lindner, Y. Mambrini, M. Pierre et al., _The waning of the WIMP? A review of models, searches, and constraints_ , _Eur. Phys. J. C_ 78 (2018) 203, [1703.07364].
  * (11) L. Roszkowski, E. M. Sessolo and S. Trojanowski, _WIMP dark matter candidates and searches—current status and future prospects_ , _Rept. Prog. Phys._ 81 (2018) 066201, [1707.06277].
  * (12) XENON collaboration, E. Aprile et al., _Dark Matter Search Results from a One Ton-Year Exposure of XENON1T_ , _Phys. Rev. Lett._ 121 (2018) 111302, [1805.12562].
  * (13) LUX collaboration, D. S. Akerib et al., _Limits on spin-dependent WIMP-nucleon cross section obtained from the complete LUX exposure_ , _Phys. Rev. Lett._ 118 (2017) 251302, [1705.03380].
  * (14) D. Abercrombie et al., _Dark Matter benchmark models for early LHC Run-2 Searches: Report of the ATLAS/CMS Dark Matter Forum_ , _Phys. Dark Univ._ 27 (2020) 100371, [1507.00966].
  * (15) N. Bernal and X. Chu, _$\mathbb{Z}_{2}$ SIMP Dark Matter_, _JCAP_ 1601 (2016) 006, [1510.08527].
  * (16) S.-M. Choi and H. M. Lee, _Resonant SIMP dark matter_ , _Phys. Lett._ B758 (2016) 47–53, [1601.03566].
  * (17) S.-M. Choi, Y.-J. Kang and H. M. Lee, _On thermal production of self-interacting dark matter_ , _JHEP_ 12 (2016) 099, [1610.04748].
  * (18) N. Bernal, X. Chu and J. Pradler, _Simply split strongly interacting massive particles_ , _Phys. Rev._ D95 (2017) 115023, [1702.04906].
  * (19) B. Chauhan, _Sub-MeV Self Interacting Dark Matter_ , _Phys. Rev._ D97 (2018) 123017, [1711.02970].
  * (20) S. Mohanty, A. Patra and T. Srivastava, _MeV scale model of SIMP dark matter, neutrino mass and leptogenesis_ , _JCAP_ 03 (2020) 027, [1908.00909].
  * (21) S. Bhattacharya, P. Ghosh and S. Verma, _SIMPler realisation of Scalar Dark Matter_ , _JCAP_ 01 (2020) 040, [1904.07562].
  * (22) H. M. Lee and M.-S. Seo, _Communication with SIMP dark mesons via Z’ -portal_ , _Phys. Lett. B_ 748 (2015) 316–322, [1504.00745].
  * (23) N. Yamanaka, S. Fujibayashi, S. Gongyo and H. Iida, _Dark Matter in the Nonabelian Hidden Gauge Theory_ , in _2nd Toyama International Workshop on Higgs as a Probe of New Physics_ , 4, 2015.  1504.08121.
  * (24) Y. Hochberg, E. Kuflik and H. Murayama, _SIMP Spectroscopy_ , _JHEP_ 05 (2016) 090, [1512.07917].
  * (25) Y. Hochberg, E. Kuflik, H. Murayama, T. Volansky and J. G. Wacker, _Model for Thermal Relic Dark Matter of Strongly Interacting Massive Particles_ , _Phys. Rev. Lett._ 115 (2015) 021301, [1411.3727].
  * (26) S.-M. Choi, Y. Hochberg, E. Kuflik, H. M. Lee, Y. Mambrini, H. Murayama et al., _Vector SIMP dark matter_ , _JHEP_ 10 (2017) 162, [1707.01434].
  * (27) S.-M. Choi, H. M. Lee, P. Ko and A. Natale, _Unitarizing SIMP scenario with dark vector resonances_ , _PoS_ ICHEP2018 (2019) 036, [1811.02751].
  * (28) J. Herms, A. Ibarra and T. Toma, _A new mechanism of sterile neutrino dark matter production_ , _JCAP_ 06 (2018) 036, [1802.02973].
  * (29) E. Ma, N. Pollard, R. Srivastava and M. Zakeri, _Gauge $B-L$ Model with Residual $Z_{3}$ Symmetry_, _Phys. Lett. B_ 750 (2015) 135–138, [1507.03943].
  * (30) E. Thrane and J. D. Romano, _Sensitivity curves for searches for gravitational-wave backgrounds_ , _Phys. Rev. D_ 88 (2013) 124032, [1310.5300].
  * (31) B. Sathyaprakash et al., _Scientific Objectives of Einstein Telescope_ , _Class. Quant. Grav._ 29 (2012) 124013, [1206.0331].
  * (32) LISA collaboration, P. Amaro-Seoane et al., _Laser Interferometer Space Antenna_ , 1702.00786.
  * (33) J. Crowder and N. J. Cornish, _Beyond LISA: Exploring future gravitational wave missions_ , _Phys. Rev. D_ 72 (2005) 083005, [gr-qc/0506015].
  * (34) N. Seto, S. Kawamura and T. Nakamura, _Possibility of direct measurement of the acceleration of the universe using 0.1-Hz band laser interferometer gravitational wave antenna in space_ , _Phys. Rev. Lett._ 87 (2001) 221103, [astro-ph/0108011].
  * (35) S. Sato et al., _The status of DECIGO_ , _J. Phys. Conf. Ser._ 840 (2017) 012010.
  * (36) L. Lentati et al., _European Pulsar Timing Array Limits On An Isotropic Stochastic Gravitational-Wave Background_ , _Mon. Not. Roy. Astron. Soc._ 453 (2015) 2576–2598, [1504.03692].
  * (37) G. Janssen et al., _Gravitational wave astronomy with the SKA_ , _PoS_ AASKA14 (2015) 037, [1501.00127].
  * (38) M. Fairbairn, E. Hardy and A. Wickens, _Hearing without seeing: gravitational waves from hot and cold hidden sectors_ , _JHEP_ 07 (2019) 044, [1901.11038].
  * (39) M. Breitbach, J. Kopp, E. Madge, T. Opferkuch and P. Schwaller, _Dark, Cold, and Noisy: Constraining Secluded Hidden Sectors with Gravitational Waves_ , _JCAP_ 07 (2019) 007, [1811.11175].
  * (40) P. Schwaller, _Gravitational Waves from a Dark Phase Transition_ , _Phys. Rev. Lett._ 115 (2015) 181101, [1504.07263].
  * (41) J. Jaeckel, V. V. Khoze and M. Spannowsky, _Hearing the signal of dark sectors with gravitational wave detectors_ , _Phys. Rev. D_ 94 (2016) 103519, [1602.03901].
  * (42) A. Addazi, _Limiting First Order Phase Transitions in Dark Gauge Sectors from Gravitational Waves experiments_ , _Mod. Phys. Lett. A_ 32 (2017) 1750049, [1607.08057].
  * (43) I. Baldes, _Gravitational waves from the asymmetric-dark-matter generating phase transition_ , _JCAP_ 05 (2017) 028, [1702.02117].
  * (44) K. Tsumura, M. Yamada and Y. Yamaguchi, _Gravitational wave from dark sector with dark pion_ , _JCAP_ 07 (2017) 044, [1704.00219].
  * (45) I. Baldes and C. Garcia-Cely, _Strong gravitational radiation from a simple dark matter model_ , _JHEP_ 05 (2019) 190, [1809.01198].
  * (46) D. Croon, V. Sanz and G. White, _Model Discrimination in Gravitational Wave spectra from Dark Phase Transitions_ , _JHEP_ 08 (2018) 203, [1806.02332].
  * (47) M. Kamionkowski, A. Kosowsky and M. S. Turner, _Gravitational radiation from first order phase transitions_ , _Phys. Rev. D_ 49 (1994) 2837–2851, [astro-ph/9310044].
  * (48) E. Witten, _Cosmic Separation of Phases_ , _Phys. Rev. D_ 30 (1984) 272–285.
  * (49) C. J. Hogan, _Gravitational radiation from cosmological phase transitions_ , _Mon. Not. Roy. Astron. Soc._ 218 (1986) 629–636.
  * (50) C. Grojean and G. Servant, _Gravitational Waves from Phase Transitions at the Electroweak Scale and Beyond_ , _Phys. Rev. D_ 75 (2007) 043507, [hep-ph/0607107].
  * (51) M. M. Pierre, _Dark Matter phenomenology : from simplified WIMP models to refined alternative solutions_.  PhD thesis, U. Paris-Saclay, 2018.  1901.05822.
  * (52) S. Bhattacharya, P. Poulose and P. Ghosh, _Multipartite Interacting Scalar Dark Matter in the light of updated LUX data_ , _JCAP_ 04 (2017) 043, [1607.08461].
  * (53) S. Weinberg, _Gauge and global symmetries at high temperature_ , _Phys. Rev. D_ 9 (Jun, 1974) 3357–3378.
  * (54) L. Dolan and R. Jackiw, _Symmetry behavior at finite temperature_ , _Phys. Rev. D_ 9 (Jun, 1974) 3320–3341.
  * (55) D. A. Kirzhnits and A. D. Linde, _Symmetry Behavior in Gauge Theories_ , _Annals Phys._ 101 (1976) 195–238.
  * (56) M. E. Carrington, _The Effective potential at finite temperature in the Standard Model_ , _Phys. Rev. D_ 45 (1992) 2933–2944.
  * (57) D. J. Gross, R. D. Pisarski and L. G. Yaffe, _Qcd and instantons at finite temperature_ , _Rev. Mod. Phys._ 53 (Jan, 1981) 43–80.
  * (58) P. Fendley, _The Effective Potential and the Coupling Constant at High Temperature_ , _Phys. Lett. B_ 196 (1987) 175–180.
  * (59) J. Kapusta, _Finite temperature Field Theory_ , _Cambridge University Press_ (1989) .
  * (60) P. B. Arnold and O. Espinosa, _The Effective potential and first order phase transitions: Beyond leading-order_ , _Phys. Rev. D_ 47 (1993) 3546, [hep-ph/9212235].
  * (61) R. R. Parwani, _Resummation in a hot scalar field theory_ , _Phys. Rev. D_ 45 (1992) 4695, [hep-ph/9204216].
  * (62) M. Quiros, _Finite temperature field theory and phase transitions_ , in _ICTP Summer School in High-Energy Physics and Cosmology_ , pp. 187–259, 1, 1999.  hep-ph/9901312.
  * (63) P. Athron, C. Balázs, A. Fowlie and Y. Zhang, _PhaseTracer: tracing cosmological phases and calculating transition properties_ , _Eur. Phys. J. C_ 80 (2020) 567, [2003.02859].
  * (64) A. Kosowsky, M. S. Turner and R. Watkins, _Gravitational radiation from colliding vacuum bubbles_ , _Phys. Rev. D_ 45 (Jun, 1992) 4514–4535.
  * (65) A. Kosowsky and M. S. Turner, _Gravitational radiation from colliding vacuum bubbles: Envelope approximation to many-bubble collisions_ , _Phys. Rev. D_ 47 (May, 1993) 4372–4391.
  * (66) S. J. Huber and T. Konstandin, _Gravitational wave production by collisions: more bubbles_ , _Journal of Cosmology and Astroparticle Physics_ 2008 (Sep, 2008) 022.
  * (67) A. Kosowsky, M. S. Turner and R. Watkins, _Gravitational waves from first-order cosmological phase transitions_ , _Phys. Rev. Lett._ 69 (Oct, 1992) 2026–2029.
  * (68) M. Kamionkowski, A. Kosowsky and M. S. Turner, _Gravitational radiation from first-order phase transitions_ , _Phys. Rev. D_ 49 (Mar, 1994) 2837–2851.
  * (69) C. Caprini, R. Durrer and G. Servant, _Gravitational wave generation from bubble collisions in first-order phase transitions: An analytic approach_ , _Physical Review D_ 77 (Jun, 2008) .
  * (70) M. Hindmarsh, S. J. Huber, K. Rummukainen and D. J. Weir, _Gravitational waves from the sound of a first order phase transition_ , _Phys. Rev. Lett._ 112 (Jan, 2014) 041301.
  * (71) L. Leitao, A. Megevand and A. D. Sanchez, _Gravitational waves from the electroweak phase transition_ , _JCAP_ 10 (2012) 024, [1205.3070].
  * (72) J. T. Giblin, Jr. and J. B. Mertens, _Vacuum Bubbles in the Presence of a Relativistic Fluid_ , _JHEP_ 12 (2013) 042, [1310.2948].
  * (73) J. T. Giblin and J. B. Mertens, _Gravitional radiation from first-order phase transitions in the presence of a fluid_ , _Phys. Rev. D_ 90 (2014) 023532, [1405.4005].
  * (74) M. Hindmarsh, S. J. Huber, K. Rummukainen and D. J. Weir, _Numerical simulations of acoustically generated gravitational waves at a first order phase transition_ , _Physical Review D_ 92 (Dec, 2015) .
  * (75) C. Caprini and R. Durrer, _Gravitational waves from stochastic relativistic sources: Primordial turbulence and magnetic fields_ , _Phys. Rev. D_ 74 (Sep, 2006) 063521.
  * (76) T. Kahniashvili, A. Kosowsky, G. Gogoberidze and Y. Maravin, _Detectability of gravitational waves from phase transitions_ , _Phys. Rev. D_ 78 (Aug, 2008) 043003.
  * (77) T. Kahniashvili, L. Campanelli, G. Gogoberidze, Y. Maravin and B. Ratra, _Gravitational Radiation from Primordial Helical Inverse Cascade MHD Turbulence_ , _Phys. Rev. D_ 78 (2008) 123006, [0809.1899].
  * (78) T. Kahniashvili, L. Kisslinger and T. Stevens, _Gravitational Radiation Generated by Magnetic Fields in Cosmological Phase Transitions_ , _Phys. Rev. D_ 81 (2010) 023004, [0905.0643].
  * (79) C. Caprini, R. Durrer and G. Servant, _The stochastic gravitational wave background from turbulence and magnetic fields generated by a first-order phase transition_ , _JCAP_ 12 (2009) 024, [0909.0622].
  * (80) H. Kudoh, A. Taruya, T. Hiramatsu and Y. Himemoto, _Detecting a gravitational-wave background with next-generation space interferometers_ , _Phys. Rev. D_ 73 (2006) 064006, [gr-qc/0511145].