# MPIX Stream: An Explicit Solution to Hybrid
MPI+X Programming

Hui Zhou  Argonne National LaboratoryLemont, IL 60439, USA ,  Ken Raffenetti
Argonne National LaboratoryLemont, IL 60439, USA ,  Yanfei Guo  Argonne
National LaboratoryLemont, IL 60439, USA and  Rajeev Thakur  Argonne National
LaboratoryLemont, IL 60439, USA

(2022)

###### Abstract.

The hybrid MPI+X programming paradigm, where X refers to threads or GPUs, has
gained prominence in the high-performance computing arena. This corresponds to
a trend of system architectures growing more heterogeneous. The current MPI
standard only specifies the compatibility levels between MPI and threading
runtimes. No MPI concept or interface exists for applications to pass thread
context or GPU stream context to MPI implementations explicitly. This lack has
made performance optimization complicated in some cases and impossible in
other cases. We propose a new concept in MPI, called MPIX stream, to represent
the general serial execution context that exists in X runtimes. MPIX streams
can be directly mapped to threads or GPU execution streams. Passing thread
context into MPI allows implementations to precisely map the execution
contexts to network endpoints. Passing GPU execution context into MPI allows
implementations to directly operate on GPU streams, lowering the CPU/GPU
synchronization cost.

MPI+X, MPI+Threads, MPI+GPUs, Network Endpoints, MPIX Stream, GPU Stream

††copyright: licensedusgovmixed††conference: ; ††price: ;

##  1\. Introduction

Modern high-performance computing applications are more and more dependent on
additional runtimes besides MPI to manage the increased number of cores per
node and myriad limited on-node resources such as shared memory, network
interfaces, and computational accelerators. Increasingly, applications are
being deployed by using a hybrid MPI+X model, where X refers to a threading
runtime such as OpenMP or an accelerator runtime such as CUDA.

The first step in MPI+X was to make MPI compatible with a threading runtime.
Since 1997, MPI has introduced four thread compatibility levels:
MPI_THREAD_SINGLE, MPI_THREAD_FUNNELED, MPI_THREAD_SERIALIZED, and
MPI_THREAD_MULTIPLE. When the appropriate thread level is chosen, threaded
applications can work correctly with MPI without MPI specifically
acknowledging the runtime. MPI implementations also can be made GPU aware.
Recent MPICH (Laboratory, 2022), MVAPICH (University, 2022), and Open MPI
(Team, 2022) releases are able to detect GPU buffers without hints from users
and make MPI work without a GPU-specific MPI interface. The application can
benefit from an MPI+GPU compatibility level similar to the MPI thread levels.
Currently, the GPU compatibility level is simply assumed.

While MPI+Threads is successful on the compatibility side of MPI+X, the
performance side has been a multi-decade struggle. With MPI+Threads—in
particular, at the MPI_THREAD_MULTIPLE thread level—applications today are
still likely to meet dismal performance. This performance is due to the extra
critical sections introduced by MPI communications. Much research has been
done on both the application side (Wang and Chandramowlishwaran, 2019) and
implementation side (Amer et al., 2015; Patinyasakdikul et al., 2019; Zambre
et al., 2020) to address the performance of MPI+Threads. To reach good
performance, applications need to make sure that the communications can happen
concurrently, and the implementations need to map the communications to
multiple communication channels to allow the communication to proceed in
parallel. Without an explicit MPI interface, making the latter mapping to
match the application layer concurrency remains an art. Mismatch will result
in either incorrect results or the introduction of extra thread contention and
bad performance.

The performance story of MPI+GPUs is different from that of MPI+Threads.
Accelerators typically require special runtime to coordinate between CPU and
accelerator executions. The launching and synchronization between CPU context
and accelerator context are carried out by the accelerator runtime. A key
performance factor here is how to minimize the launching and synchronization
cost. To optimize the performance, we need MPI operations to be enqueued to an
accelerator execution context and then let the accelerator runtime manage its
actual execution. In order to realize this new mode of MPI operations, new MPI
interfaces that work directly with accelerator execution context are needed.

A common theme from the pursuit of performance in MPI+X is the need for MPI to
have the concept of execution context. In this paper we survey the current
status of MPI+Threads and MPI+GPUs and propose a new MPI concept, called MPIX
stream, that can be used to represent execution context from other runtimes.
MPIX stream allows explicit coordination for MPI+Threads and enables direct
GPU runtime operation for MPI communications.

##  2\. Background

###  2.1. MPI’s Execution Model

Among the four MPI thread levels, only MPI_THREAD_MULTIPLE requires
clarification on its execution model. The other three are all serial execution
model enforced by applications. For MPI_THREAD_MULTIPLE, the current guiding
principle is as follows: “When a thread is executing one of these (MPI)
routines, if another concurrently running thread also makes an MPI call, the
outcome will be as if the calls executed in some order” (Message Passing
Interface Forum, 2021). In another words, MPI is an outcome-dictated serial
execution model.

A naïve implementation of the serial execution model is to impose a global
critical section for every MPI call and yield only during its progress loop.
With the global critical section, if two threads concurrently call MPI
communication functions, both threads not only are serialized at the point of
the MPI operation but also incur a significant cost from the synchronization,
resulting in performance worse than if all communications are called from a
single thread. However, an implementation is allowed to “optimize,” making
some parts or whole communications parallel as long as the outcome is not
affected.

What is an outcome in MPI? Unfortunately, this is not clearly specified in the
MPI standard and thus is a source of ambiguity and debate. Nevertheless, some
consensus has been reached on what is and what is not an outcome. For example,
a message delivery order is not an MPI outcome. A second sequentially issued
message is allowed to be delivered before the first one. On the other hand, a
message matching order is an MPI-defined outcome. Two sequentially issued
sends that both match the same receive are guaranteed to match the first one
before the second one. If we ignore outcomes that are outside the MPI
standard, then an MPI implementation can execute some of the communications in
parallel as long as there is no matching order between them. For example,
messages issued from different communicators are matched independently. If
assertions that no wildcard tag matching will be used, then messages using
different tags can also be carried out in parallel.

This strategy raises two issues. First, the consensus is not universal.
Certainly outcomes occur outside those specified by the MPI standard. For
example, what if the serialization or parallelization of the communication is
an application-intended outcome? Second, communicators or tags are not perfect
identification of concurrent execution context. Using communicators or tags to
sideload the expression of parallelism may result in convoluted code yet still
not necessarily be able to achieve perfect parallelization.

###  2.2. Network Endpoints

Modern high-speed interconnection fabrics are designed with capabilities to
support communication by multiple execution threads. Generally, this is done
by allowing separate fabric resources to be allocated. In this paper we refer
to these allocated fabric resources as network endpoints. Communications can
be carried out concurrently from separate endpoints. The network endpoints are
abstractions over hardware capability and may include software contexts such
as address table, message queues, and completion event queues. In libfabric
(Open Fabric Alliance, 2022), a network endpoint may be represented by a
domain, an endpoint, or a scalable endpoint. In UCX (Open UCX, 2022), a
network endpoint is typically represented by a UCP worker.

MPI implementations that utilize network endpoints also need to allocate their
own internal communication context to isolate global states that are needed
during a communication. For best performance, these implementation-level
contexts need to be matched to the network endpoints. Both MPICH and Open MPI
have implemented such communication contexts. In MPICH, it is referred to as
the virtual communication interface (VCI) (Zambre et al., 2021). In Open MPI,
it is referred to as the communication resources instance (CRI)
(Patinyasakdikul et al., 2019). In this paper we generally refer to these
contexts as network endpoints.

Network endpoints are a finite resource. More endpoints beyond a hardware’s
capability will be serialized at the hardware anyway and will incur more
overhead in managing the multiple endpoints. A limit is often imposed by a
network library and sometimes by a network driver. It is common to have a
limit matching the number of cores in a node.

Concurrent access to a single network endpoint is not allowed, or it will
result in data race and state corruption. Thus, a critical section around the
access of each network endpoint is necessary unless it can be guaranteed that
concurrent usages will not occur.

###  2.3. Nonlocal Nature of Communication

While both thread contexts and network endpoints are local process concepts, a
communication necessarily involves a pair of network endpoints from both the
local process and the remote process. When one does not specify a network
endpoint in a communication, as is the case with the current MPI standard, the
implementation chooses a default network endpoint for both the local process
and remote process. If this default choice is a constant choice, then all
communications are serialized on both the sender side and the receiver side.
Implicit schemes or semi-explicit schemes via hints can be used to hash the
network endpoints’ choice in order to achieve parallel communications. The
hashing algorithm must be deterministic and consistent for both the sender
side and receiver side. If the sender side sends a message to a remote
endpoint that does not expect to receive it, either the message will get lost
or heavy synchronization will occur in order to move the message to the
receiving context.

Implicit hashing schemes often employ certain enforcement policies. A typical
policy is to enforce a one-to-one mapping of endpoints. If we assign a
sequential id to each network endpoint in a process, then a one-to-one mapping
allows communications only between network endpoints with the same id. With
this policy, network endpoints can be easily determined by hashing common
values between sender and receiver, such as communicator id, sender and
receiver ranks, and tags.

Another policy is to allow the sender to send from any network endpoint but
receive only from a default network endpoint. With this policy, the sender
side can easily achieve concurrent sends by hashing local information or even
by random assignment. Since messages are all received by a single network
endpoint, however, the overall message rates are limited by the single
receiving thread.

The two policies match to the two common communication patterns illustrated in
Figure 1: the one-to-one pattern and the N-to-1 pattern. In a one-to-one
pattern, one thread from one process communicates only to one thread in other
processes. An example of a one-to-one pattern is the stencil application,
where a partition for a single thread only shares a halo region with another
thread in the neighbor processes (see Figure 2). Note that the pairing of
threads may depend on the geometry layout and may not correspond to the
ordering of the thread numbers. In an N-to-1 pattern, multiple threads may
send messages to other processes, but a single thread in each process is
dedicated to receive all messages. An example of an N-to-1 pattern is a task-
based application, where multiple threads run tasks and generate events and
where a single progress thread polls and responds to events.

Figure 1. Typical MPI+Threads communication patterns: (a) one-to-one pairwise
mapping; (b) N-to-1 mapping. Figure 2. Communication patterns in a 2-D stencil
partition.

###  2.4. GPU Queuing Stream

Graphics processing units (GPUs) are common components today in high-
performance computing clusters. GPU programming imposes its own execution
model. A common concept in GPU programming is an abstract execution queue. In
NVIDIA’s CUDA runtime, this is referred to as a CUDA stream. Unlike a thread,
calls do not directly run on the execution queue. Instead, the operations are
enqueued, and the GPU runtime will dispatch the operations to GPU kernels
asynchronously. The execution queue can be generalized as an execution graph
with explicit dependency. An execution graph may allow the GPU runtime more
flexibility in optimizations. A sequential queue, on the other hand, often
contains explicit synchronizations that may prevent optimization.

With GPU-aware MPI implementations (Wang et al., 2014), MPI functions can be
called to send from or receive into GPU memory directly. However, MPI
functions are still implemented to execute entirely in a CPU context. Thus,
full CPU/GPU synchronizations are necessary. When the application relies on
the GPU running computations, the full CPU/GPU synchronization imposes a great
performance penalty.

To break the performance bottleneck, MPI needs to enable partial CPU/GPU
synchronization, directly embed synchronization into the GPU kernel, or
explicitly expose dependency by using a facility such as a CUDA graph. One way
to do so is to pass the GPU runtime execution context to MPI so that the MPI
implementation can directly operate under the GPU context. In the case of
CUDA, this may mean passing the CUDA stream object into MPI send and receive
functions.

In the GPU execution context, examples including CUDA stream and SYCL queue
are often asynchronous. This poses another challenge on working with MPI’s
execution model. How do we assess matching order when the issuing order of the
operation itself is not deterministic? Unlike in the MPI+Threads case where a
default serialization model makes sense, MPI+GPUs will need a new concept to
accommodate the foreign execution context.

##  3\. Proposal of MPIX Stream

We propose a new MPI object, called MPIX stream, and a set of new APIs to
allow users to explicitly identify their thread context and GPU queuing
context. We use the prefix MPIX instead of MPI to refer to objects and
functions that are proposed in this paper and not (yet) officially part of
MPI.

###  3.1. MPIX Stream

First, we introduce MPIX stream as an abstract concept to facilitate a direct
mapping from user-level runtime execution contexts to MPI execution contexts.
To MPI, an MPIX stream represents a local serial execution context. Any
runtime execution contexts outside MPI, as long as the serial semantic is
strictly followed, can be associated to an MPIX stream. Examples include
kernal threads, user-level threads, GPU queuing streams, or even code across
multiple threads with serialized synchronizations.

To illustrate, in the following two pseudocode listings, we use
MPIR_SEND_ON_STREAM to represent an arbitrary MPI operation, and stream_1
represents a specific MPIX stream object that is associated with the MPI
operation. Because in both listings, the MPI operations are strictly
serialized, both are valid usages of MPIX stream.

Listing 1: Pseudocode using MPIX stream within a single thread

⬇

{

/* within a single thread */

MPIR_SEND_ON_STREAM(stream_1, msg1);

MPIR_SEND_ON_STREAM(stream_1, msg2);

}

Listing 2: Pseudocode using MPIX stream from two threads with explicit thread
synchronization

⬇

{

/* Thread 1 */

MPIR_SEND_ON_STREAM(stream_1, msg1);

THREAD_BARRIER();

}

{

/* Thread 2 */

THREAD_BARRIER();

MPIR_SEND_ON_STREAM(stream_1, msg2);

}

Note that pseudocode is used in these examples, and in particular,
MPIR_SEND_ON_STREAM is not a proposed API. We will describe the actual APIs
later in this section. The examples also illustrate that it is the users’
responsibility to map their thread context to the MPIX stream. This relieves
the burden of MPI having to synchronize between threads, which can be thread
runtime dependent.

To create an MPIX stream, one calls MPIX_Stream_create.

⬇

int MPIX_Stream_create(MPI_Info info, MPIX_Stream *stream)

Info hints can be used to create implementation-supported special streams, for
example, a CUDA stream. Otherwise, MPI_INFO_NULL can be used. Whether the
returned stream is backed by distinct network endpoints is implementation
dependent. For MPI thread contexts, we recommend allocating unique network
endpoints for each new stream. This approach allows programming the
application with a predictable performance. The implementation should return
failure if it runs out of network endpoints. With unique endpoints and the
strict serial execution context, the implementation may safely skip critical
sections in the communication path. At the extreme end of the strong scaling,
even an uncontended critical section can be too expensive.

The implementation may also assign a single network endpoint to multiple MPIX
streams, allowing applications to create more streams than available network
endpoints. For example, network endpoints can be assigned to a newly created
stream in a round-robin fashion. This may provide flexibility for some
applications. We note that a per-endpoint critical section is necessary to
prevent concurrent access to network endpoints.

Because network endpoints are finite resources, users should free the stream
to make the resource available for future allocation.

⬇

int MPIX_Stream_free(MPIX_Stream *stream)

The network resource can be deallocated only when all the operations using the
stream have been completed. We imagine that whether a deallocation is
successful is important feedback for users. In particular, a failed or delayed
deallocation may prevent a future MPIX_Stream_create from succeeding. Thus,
MPIX_Stream_free may fail with an appropriate error code if the internal
resource deallocation cannot be completed.

GPU streams emphasize lightweight synchronization between the CPU and
accelerator. Thus, having concurrent CPU communications may not be as
important as in multithreaded programming. In addition, an implementation may
choose to use a dedicated CPU thread to progress all GPU stream
communications. Thus, GPU streams are likely to be assigned with duplicate
network endpoints. On the other hand, having dedicated network endpoints for
GPU stream-related communications is likely beneficial so that the progress
thread needs to poll only a few network endpoints rather than polling global
progress and contending with traffic on other CPU threads.

For backward compatibility, MPIX_STREAM_NULL is defined. Operations on
MPIX_STREAM_NULL have the same semantics as do conventional operations without
explicit streams.

###  3.2. Passing Opaque Binary Info Hints

Currently, an MPI_Info object supports values only as strings. A GPU queuing
object not only is not a string but is often opaque to the user. For example,
is a CUDA stream an integer or a pointer, or could it be neither? To pass an
opaque binary as a string, we need some encoding scheme that users can use to
encode and implementations can consistently decode. We propose a new MPI
function for this purpose.

⬇

int MPIX_Info_set_hex(MPI_Info info, const char *key, void *value, int vallen)

An implementation can choose any binary to ASCII encoding to implement this
function. For completeness, there should be a corresponding MPIX_Info_get_hex
function. For this proposal, however, we focus on how to let users pass the
GPU queuing object into MPI. Since getting the object back from MPI is not
critical, we are not proposing the retrieving function here.

###  3.3. MPIX Stream Communicator

Once an MPIX stream is created, it is possible to define MPI operation APIs
that directly use the stream as an argument to each operation. This has two
drawbacks, however. First, we will need to create a new API for every MPI
operation, from MPI_Send, MPI_Isend, to MPI collectives and MPI remote memory
access (RMA) operations. Not only does it take considerable effort to maintain
an inflated standard, but it is also a burden for users to learn these APIs.
Second, it is not sufficient to add only a local stream to an operations
argument. We must add an argument for a remote stream as well unless we want
to restrict to an arbitrary policy. For collectives, this may require an array
of stream arguments, one for every participating process. Even for one-sided
RMA operations, a specific targeting stream may be critical if the user wants
to dedicate a progress thread to drive passive progress. Unlike local streams,
identifying remote streams per operation is cumbersome for the user.

Thus we propose the MPIX stream communicator.

⬇

int MPIX_Stream_comm_create(MPI_Comm parent_comm, MPIX_Stream stream, MPI_Comm
*stream_comm)

This is a collective operation. Stream information from all processes or its
network endpoint address can be Allgathered and stored locally. All
conventional MPI operations can be issued to a stream communicator without
additional parameter changes.

If the parent_comm is also a stream communicator, it is treated as a normal
communicator. That is, the stream attached to the parent_comm is discarded in
the new communicator.

Each stream is still local to each process, and streams do not need to agree
in any aspects between processes. In particular, any process is allowed to use
MPIX_STREAM_NULL in constructing the stream communicator. The operation
context is dictated by each process’s attached stream. Unless it is
MPIX_STREAM_NULL, a strict serial context per stream is assumed.

One also can mix normal streams with GPU streams. Once again, the operation
context is dictated by its local stream. It is possible to have a collective
where one process is posting immediate operations when another process is
enqueuing its corresponding operation, and it is up to the GPU runtime to
asynchronously execute them.

###  3.4. MPIX Enqueue APIs

Because an attached local stream dictates the operation mode, it is feasible
to use the same conventional MPI operations without syntax change for GPU
enqueue functions. However, there are subtle semantic differences between an
enqueued operation and a non-enqueued operation. To an advocate of explicit
coding style, having, for example, MPI_Send enqueue to a GPU stream can seem a
very bad code. Thus we propose the following specific enqueue APIs.

⬇

int MPIX_Send_enqueue(const void *buf, int count, MPI_Datatype datatype, int
dest, int tag, MPI_Comm comm)

int MPIX_Recv_enqueue(void *buf, int count, MPI_Datatype datatype, int source,
int tag, MPI_Comm comm, MPI_Status *status)

int MPIX_Isend_enqueue(const void *buf, int count, MPI_Datatype datatype, int
dest, int tag, MPI_Comm comm, MPI_Request *request)

int MPIX_Irecv_enqueue(void *buf, int count, MPI_Datatype datatype, int
source, int tag, MPI_Comm comm, MPI_Request *request)

int MPIX_Wait_enqueue(MPI_Request *request, MPI_Status *status)

int MPIX_Waitall_enqueue(int count, MPI_Request array_of_requests[],
MPI_Status array_of_statuses[])

These routines have identical signatures as their conventional counterparts.
It is an error to call the enqueue functions if the communicator is not a
stream communicator or does not have a local GPU stream attached. We
intentionally omitted MPIX_Waitsome and MPIX_Waitany, as well as the various
test functions, because the nondeterministic nature of these functions does
not work with the enqueue environment. Moreover, MPIX_Waitall_enqueue must
have requests all issued on the same local stream.

There are semantic differences between enqueuing APIs and conventional APIs.
For example, MPIX_Send_enqueue, as with all enqueuing APIs, returns
immediately after registering the operation. A separate progress thread, which
may be the GPU runtime thread, will initiate and complete the communication
asynchronously. This is different from the conventional nonblocking API, e.g.,
MPI_Isend. An implementation of MPI_Isend may process the message buffer and
initiate the communication immediately before return, while both
MPIX_Send_enqueue or MPIX_Isend_enqueue will only process the buffer and
initiate communication after the previous enqueued operations complete. Both
MPIX_Send_enqueue and MPIX_Isend_enqueue will return immediately after
registering the operation, but MPIX_Isend_enqueue allows the following
enqueued functions to proceed before the communication completes with a
corresponding MPIX_Wait_enqueue. It may appear confusing since both the
traditional non-blocking APIs and the new enqueue APIs are asynchronous. It is
worth noting that they are dealing with two orthogonal kinds of synchronicity.
Traditional non-blocking APIs are mostly concerned with synchronizing data in
the message buffer, while the enqueue APIs are concerned with synchronization
between execution contexts. The former is synchronized by calls such as
MPI_Wait, while the latter is sychronized by calls such as
cudaStreamSynchronize. In particular, with the addition of the enqueue APIs,
GPU synchronization calls, such as cudaStreamSynchronize, are no longer needed
for message data or communication synchronizations.

The enqueue APIs can be extended to collectives and RMA functions. All the
extended enqueue functions will have identical function signatures as their
conventional counterparts. They are not listed here due to the large number of
them. For collectives, if some of the processes are not associated with an
enqueuing stream, then those processes should call the conventional non-
enqueue API. For RMA, enqueuing operations include window synchronizations.

###  3.5. MPIX Multiplex Stream Communicator

The stream communicator works well for codes that communicate only to a single
thread on a remote process. The communicator effectively constructs a thread
communication group and, similar to legacy MPI code, can achieve concurrent
multithread performance without any change to the code.

On the other hand, when a thread needs to communicate with multiple threads of
a remote process, it may be necessary to create and manage multiple stream
communicators. Doing so can be cumbersome and become unmanageable quickly. For
example, two processes each with 4 threads will need 16 stream communicators
to achieve an all-to-all communications.

To address this usability issue, we propose the MPIX multiplex stream
communicator.

⬇

int MPIX_Stream_comm_create_multiple(MPI_Comm parent_comm, int count,
MPIX_Stream local_streams[], MPI_Comm *stream_comm)

With a multiplex stream communicator, each process can attach multiple local
streams. To use the MPIX multiplex stream communicator, we propose the
following point-to-point APIs.

⬇

int MPIX_Stream_send(const void *buf, int count, MPI_Datatype datatype, int
dest, int tag, MPI_Comm comm, int src_idx, int dst_idx)

int MPIX_Stream_recv(void *buf, int count, MPI_Datatype datatype, int source,
int tag, MPI_Comm comm, int src_idx, int dst_idx, MPI_Status *status)

int MPIX_Stream_isend(const void *buf, int count, MPI_Datatype datatype, int
dest, int tag, MPI_Comm comm, int src_idx, int dst_idx, MPI_Request *request)

int MPIX_Stream_irecv(void *buf, int count, MPI_Datatype datatype, int source,
int tag, MPI_Comm comm, int src_idx, int dst_idx, MPI_Request *request)

These APIs allow users to explicitly address local and remote streams via an
index. This index can be thought of as a rank for threads. Unlike the source
rank in send and destination rank in receive, both src_idx and dst_idx are
needed in the argument list since the communicator here does not uniquely
identify the local thread. MPIX_ANY_INDEX can be used to support a wildcard
receive.

The multiplex stream communicator is especially useful for an N-to-1
communication pattern, where one or a few polling threads are responsible for
receiving messages sent from any other thread. Without multiplex stream
communicators, one must create multiple single-stream communicators and have
the polling thread poll each communicator in turn. With multiplex stream
communicators, the polling thread needs to poll only a single communicator.

##  4\. Comparison with Previous Efforts

Researchers have been addressing the paradigm of MPI+Threads for two decades.
Many previous efforts have been reported, and this proposal draws experience
from these past efforts. In this section we compare our work with the previous
proposals and alternative solutions.

###  4.1. Implicit Method

The implicit method (Patinyasakdikul et al., 2019; Zambre et al., 2020) builds
on the thesis that MPI already has a sufficient mechanism to allow users to
express the inherent parallelism in their communications. Based on the
outcome-dictated serial execution model, any operations that do not affect
MPI-specified outcomes are candidates for deserialization. In particular,
users can express parallelism using separate communicators.

Using distinct communicators for different thread communication groups is
similar to the usage of stream communicators in our proposal. However, to
write code to fully take advantage of implicit methods requires knowledge of
the particular MPI implementation, and the behavior is never guaranteed. For
example, for a one-to-one pattern, if the implementation uses a different
policy other than per-communicator mapping, the performance will not be ideal.
In contrast, with our proposal, the mapping to network endpoints is explicitly
identified with the MPIX stream, and its behavior is guaranteed.

When the communication does not fit into a one-to-one pattern, one thread may
need to use more than one communicator. In a typical implementation of the
implicit method, where network endpoints are assigned to communicators in a
round-robin fashion, when each thread uses multiple communicators, two threads
may still be assigned with the same network endpoint despite using different
communicators. Using MPIX stream, however, network endpoint assignments are
explicit. Thus, using multiple communicators per thread will not be an issue.

One must use fine-grained critical sections with the implicit method. When
there is no thread contention, there will be a slight performance penalty
compared with using a global critical section because of the overhead of using
more locks. On the other hand, correct usage of MPIX stream allows
implementation to skip locking altogether, resulting in a performance gain.

###  4.2. MPI Endpoints Proposal

The MPIX stream proposal has its roots in the MPI endpoints proposal (Dinan et
al., 2014). The multiplex stream communicator is nearly identical to an
endpoints communicator. Both proposals allows direct addressing of individual
endpoint or thread context.

A key flaw in the endpoint proposal, from our view, is the inflation of thread
context into virtual processes. To a multithread programmer, a process and a
thread are separate concepts. Between threads the memory is shared, and thus
there is no need for explicit data exchange. Between processes, explicit
messaging is needed, and that is where MPI is used. Thus, the process concept
and its identification using ranks are important. With the endpoints proposal,
the process become less identifiable. Users may have to manage their own
endpoint ranks to process the rank translation table. The endpoints proposal
also makes interthread messages equally accessible as interprocess messages.
Since users rarely need interthread messages, this inflation makes MPI more
difficult to understand and use.

In contrast, a multiplex stream communicator maintains the address via process
ranks plus the thread index. This fits the model of multithreaded application
naturally, and thus it is easier to learn and use.

###  4.3. MPI-4 Partitioned Communication

Partitioned communication is a new addition to the MPI-4 standard. One of its
motivations is to provide simpler and more effective multithread optimization
(Grant et al., 2019). Partitioned communication has an explicit init stage
where implementations can set up strategy and decide network endpoints mapping
to partitions. The actual communications can be triggered by MPI_Pready calls,
which can occur concurrently or out of order.

This should be compared with sending multiple messages from multiple threads,
each message corresponding to a single partitioned data. Using explicit
streams, users always can control the thread-to-message mapping, thus
achieving the desired parallelism. Partitioned communication, on the other
hand, is still an implicit mapping mechanism, although through the init stage
implementations may be able to achieve better mapping than implicit static
mapping can.

However, partitioned communication focuses on the optimization of a single
parallel region with a single coordinated functionality where individual
partitioned data is part of bigger data that can be described with a single
message. It does not solve the concurrency issue when there are other messages
besides the partitioned message or when the message cannot be equally
partitioned. For example, an application may need to exchange data on an
irregular region, which cannot use a partitioned scheme.

MPIX stream, on the other hand, lets users explicitly control thread mapping
and thus orchestrate communications across multiple areas and even orchestrate
beyond single parallel regions. The latter is critical when applications have
dedicated threads taking care of some of the communication needs.

###  4.4. MPI-4 Sessions

Another new addition to the MPI-4 standard is MPI Sessions (Hjelm et al.,
2019). Intuitively from the name, an MPI session is intended as a local
context to encapsulate all MPI local objects. It can support independent,
including concurrent, MPI operations from different MPI sessions.

Most general-purpose objects can be specialized for a particular purpose.
During proposal development, the MPI Forum quickly discovered that MPI
sessions can be used to solve all kinds of issues. Simply by attaching special
attributes or hints to a session, it can become specific enough to address
almost any issue that requires specific context.

For example, an MPI session can be used as an MPIX stream. To do so, we need
to create an MPI session with MPI_THREAD_SERIALIZED thread level, then
identify the session with a local thread context. We can also attach a GPU
queueing object to the session using info hints. The session can proceed to
create communicators, which essentially are stream communicators. From here,
all the semantics proposed for stream communicators can proceed. The only
missing functionality is the multiple stream communicator.

Although this solution via MPI Sessions works, it is not intuitive and is
convoluted. It is difficult to connect the name “session” as an equivalent to
a thread, stream, or network endpoint. To educate users on such session usages
will be difficult. To use sessions in this way, each session needs to go
through a bootstrapping stage until arriving at a communicator. Compared with
the pattern that we initialize, setting up the common part including the
initial collective communications, and then creating specialized communicators
for special code regions, the code using an MPI session can be difficult to
manage because the code is less separated by purposes.

Unless the MPI standard specifies special usage, an implementation is unlikely
to optimize the usage of an MPI session as an MPIX stream.

###  4.5. NVIDIA Collective Communication Library (NCCL)

NCCL (NVIDIA, 2022) is a library providing selected point-to-point and
collective communications between GPU devices. NCCL’s design is optimized for
dense multi-GPU systems, while MPI is geared toward interprocess
communications across many nodes in a cluster. NCCL uses a stream-based
enqueue-only API, which our proposed enqueue APIs are directly modeled after.
NCCL being a CUDA-specific library can directly use CUDA stream in its
interface. This tie to a specific external runtime is undesirable for MPI,
which need be usable for a variety of GPU runtimes. With MPIX_Stream, the
specific tie is limited to the info hints during stream creation, while the
rest of the APIs are portable across different GPU runtimes.

Beside the enqueuing streams, NCCL focuses on communications between GPU
devices, and its communicator is formed by a collection of GPU devices. On the
other hand, MPI communications are interprocess by default, although self
messages are also allowed. An MPI communicator is always formed by a group of
processes, and each rank addresses a single process.

Lastly, NCCL, being a specialized library, only supports contiguous buffers
with intrinsic datatypes, and its operations are limited to a selected set of
operations. Our proposed enqueue APIs, on the other hand, work for MPI
datatypes, and can be readily extended to all MPI operations.

###  4.6. Alternative Proposal for GPU Enqueues

Alternative proposals to add GPU enqueue operations to MPI are to directly add
the GPU queue objects to point-to-point operations as extra arguments.

⬇

int MPIX_Send_enqueue(const void *buf, int count, MPI_Datatype datatype, int
dest, int tag, MPI_Comm comm, enum MPIX_QUEUE_TYPE type, void *stream)

int MPIX_Recv_enqueue(void *buf, int count, MPI_Datatype datatype, int source,
int tag, MPI_Comm comm, enum MPIX_QUEUE_TYPE type, void *stream, MPI_Status
*status)

int MPIX_Isend_enqueue(const void *buf, int count, MPI_Datatype datatype, int
dest, int tag, MPI_Comm comm, enum MPIX_QUEUE_TYPE type, void *stream,
MPI_Request *request)

int MPIX_Irecv_enqueue(void *buf, int count, MPI_Datatype datatype, int
source, int tag, MPI_Comm comm, enum MPIX_QUEUE_TYPE type, void *stream,
MPI_Request *request)

int MPIX_Wait_enqueue(MPI_Request *request, enum MPIX_QUEUE_TYPE type, void
*stream, MPI_Status *status)

int MPIX_Waitall_enqueue(int count, MPI_Request array_of_requests[], enum
MPIX_QUEUE_TYPE type, void *stream, MPI_Status array_of_statuses[])

This is essentially the same as skipping the MPIX stream creation and stream
communicator creation and moving the info hint from the stream directly to
each enqueue operation. If we look at only the listed functions, this is a
simpler and more direct way of achieving it.

The MPIX stream uses separate steps to create the stream and then the stream
communicators. Thus it has more opportunities for implementations to validate
and optimize. It is also easily extensible by accepting more info hints. It
readily extends the functionality to collectives and one-sided communications
without extra APIs. Also, the stream creation and communicator construction
provide error check opportunities so users can choose to use fallback
algorithms if necessary. Cosmetically, passing an opaque object via a void
pointer reference is less desirable.

Namashivayam et al. in a recent study (Namashivayam et al., 2022) proposed a
new data object, MPIX_Queue, which is similar to MPIX_Stream but limited as an
abstraction over GPU stream objects. They proposed a similar set of APIs using
MPIX_Queue as a direct argument. Compared to our proposal, it shares the
benefit of MPIX_Stream abstraction, but lacks the extensibility from the
stream communicator construction.

##  5\. Prototype Implementation

We have implemented a prototype of the stream APIs proposed in this paper. The
prototype is available in the MPICH 4.1a1 release.

###  5.1. Mapping VCI to MPIX Stream

MPICH internally already maintains a pool of virtual communication interfaces.
With the per-VCI critical section model, each VCI uses separate mutexes and
accesses dedicated network endpoints. Communications from separate VCIs can be
fully concurrent. For a detailed discussion of VCIs, see (Zambre et al.,
2020).

MPICH currently supports implicit VCI hashing for MPI_THREAD_MULTIPLE using
traditional APIs. For our prototype implementation, we separate the pool of
VCIs into an implicit pool and an explicit pool. The size of each pool can be
controlled by the user via MPI tool interface control variables. The total
number of VCIs is limited by both available network endpoints and software
limits. More VCIs will demand larger tables in various internal data
structures, as well as heavier address exchange during initialization, so we
advise users to set the size of VCI pools according to their application
usage. For example, if the application is not using the stream APIs proposed
in this paper, users should leave the reserved VCI pool size at 0 and set the
implicit VCI pool size to match the number of threads or number of cores they
are using. On the other hand, if the application is using the stream APIs, we
expect users want to control their stream mapping explicitly throughout the
application. Thus, they should leave the implicit VCI pool size at the
default, 1, and set the reserved VCI pool size according to the total number
of allocated streams.

All the functions proposed in this paper are implemented; however, not all
functionality is complete. In particular, one-sided operations are not
explicitly stream-aware. A window created by using a stream communicator will
behave like a conventional communicator with implicit VCI assignment. Point-
to-point functions and collective functions, including nonblocking and
persistent variations, are fully stream-aware.

###  5.2. GPU Enqueue APIs

The GPU enqueue functionality is implemented only for CUDA and only through
the explicit point-to-point enqueue functions. The work to extend the
functionality to collectives and one-sided communication is ongoing.

The current implementation uses CUDA’s cudaLaunchHostFunc to enqueue the MPI
operation to the CUDA stream. Not all GPU runtimes provide host enqueue
functions; and even with CUDA, this is not optimal. The current CUDA
implementation incurs a heavy switching cost for cudaLaunchHostFunc.

A better implementation may use a dedicated host thread to the progress
operation queue and enqueue only the event triggers or event synchronizations
to the kernel queues. Having MPI launch separate progress threads is another
area worth exploring. Having each runtime spawn hidden progress threads behind
the user’s knowledge is often not optimal.

###  5.3. Results

To verify our implementation, we conducted a microbenchmark measurement on an
Intel Skylake cluster in the Joint Laboratory for System Evaluation (JLSE) at
Argonne National Laboratory. The nodes in the cluster are connected by
Mellanox InfiniBand EDR. The microbenchmark launches a number of threads, and
each thread then sends 8-byte messages to a corresponding thread on another
process. Each thread communicates using a per-thread communicator. The results
are shown in Figure 3. Three measurements were taken. The first was with MPICH
configured to use the global critical section. Shown in the red curve, we see
that the total message rate drops as soon as more threads start to compete for
the critical section.

The green curve shows the message rate with MPICH configured to use the per-
VCI critical section. Here traditional MPI communicators are used, and thus
MPICH is implicitly hashing the communications using different VCIs. The
microbenchmark is designed to achieve perfect implicit hashing, and we see
good scaling as we increase the number of threads. Note that the message rate
with a single thread is actually smaller than the corresponding message rate
with the global critical section. The reason is that the per-VCI critical
sections are finer grained and it often takes multiple critical sections along
the communication path—in particular, the receive path and progress engine—for
each message to complete. Even without contention, the extra locking and
unlocking hurt the performance.

The third curve, shown in blue, is from rewriting the benchmark to use the
proposed stream communicators. Each stream communicator is attached with a
unique MPIX stream object per thread. Because the semantics of MPIX stream
guarantees a serial execution context, our implementation is able to
completely remove locking, resulting in around $20\%$ gain in the total
message rate up to 20 threads.

In our current implementation, atomic variables and atomic operations are
still used to reference count request objects and completion flags. Even
uncontended atomics hurt performance in these microbenchmarks. Unfortunately,
the current MPICH code structure made it difficult to switch off the atomic
operations. This issue remains on our to-do list.

Figure 3. Multithread message rate on 8-byte messages using
MPI_Isend/MPI_Irecv.The message rate using MPIX_Stream is around $20\%$ higher
than with implicit VCIs.

##  6\. Examples

In this section we show two examples to illustrate how the proposed APIs will
be used. In Listing 3 we show an example hybrid MPI+OpenMP program for a one-
to-one thread communication pattern. Each thread is represented by a unique
MPIX_Stream and uses a dedicated stream communicator to send and receive
messages. A good implementation can ensure these communications happen
concurrently without incurring extra locking cost.

Listing 3: Example using MPIX stream communicator with OpenMP

⬇

#define NT 4

int main(void) {

int rank;

int tl;

MPI_Init_thread(NULL, NULL, MPI_THREAD_MULTIPLE, &tl);

MPI_Comm_rank(MPI_COMM_WORLD, &rank);

MPIX_Stream streams[NT];

MPI_Comm comms[NT];

for (int i = 0; i < NT; i++) {

MPIX_Stream_create(MPI_INFO_NULL, &streams[i]);

MPIX_Stream_comm_create(MPI_COMM_WORLD, streams[i], &comms[i]);

}

#pragma omp parallel num_threads(NT)

{

int id = omp_get_thread_num();

char buf[100];

int tag = 0;

if (rank == 0) {

MPI_Send(buf, 100, MPI_CHAR, 1, tag, comms[id]);

} else if (rank == 1) {

MPI_Recv(buf, 100, MPI_CHAR, 0, tag, comms[id], MPI_STATUS_IGNORE);

}

}

for (int i = 0; i < NT; i++) {

MPIX_comm_free(&comms[i]);

MPIX_Stream_free(&streams[i]);

}

MPI_Finalize();

return 0;

}

In Listing 4 we show an example MPI+CUDA program using CUDA’s asynchronous
stream enqueue facility. It is a simple vector computation, SAXPY. Process 0
generates a portion of the data and sends it to process 1, which launches the
kernel to do the computation after receiving the data. All memory copies, MPI
send/receive, and computation kernels are asynchronously launched to a user-
supplied CUDA stream.

Listing 4: Example using MPIX stream for CUDA stream enqueue operations

⬇

const float a_val = 2.0;

const float x_val = 1.0;

const float y_val = 2.0;

__global__

void saxpy(int n, float a, float *x, float *y)

{

int i = blockIdx.x*blockDim.x + threadIdx.x;

if (i < n) y[i] = a_val * x[i] + y[i];

}

int main(void)

{

cudaStream_t stream;

cudaStreamCreate(&stream);

int rank;

MPI_Init(NULL, NULL);

MPI_Comm_rank(MPI_COMM_WORLD, &rank);

float *x, *y, *d_x, *d_y;

MPI_Info info;

MPI_Info_create(&info);

MPI_Info_set(info, "type", "cudaStream_t");

MPIX_Info_set_hex(info, "value", &stream, sizeof(stream));

MPIX_Stream mpi_stream;

MPIX_Stream_create(info, &mpi_stream);

MPI_Info_free(&info);

MPI_Comm stream_comm;

MPIX_Stream_comm_create(MPI_COMM_WORLD, mpi_stream, &stream_comm);

/* Rank 0 sends x data to Rank 1, Rank 1 performs a * x + y and checks result
*/

if (rank == 0) {

x = (float*)malloc(N*sizeof(float));

for (int i = 0; i < N; i++) {

x[i] = x_val;

}

MPIX_Send_enqueue(x, N, MPI_FLOAT, 1, 0, stream_comm);

free(x);

} else if (rank == 1) {

y = (float*)malloc(N*sizeof(float));

cudaMalloc(&d_x, N*sizeof(float));

cudaMalloc(&d_y, N*sizeof(float));

for (int i = 0; i < N; i++) {

y[i] = y_val;

}

cudaMemcpyAsync(d_y, y, N*sizeof(float), cudaMemcpyHostToDevice, stream);

MPIX_Recv_enqueue(d_x, N, MPI_FLOAT, 0, 0, stream_comm, MPI_STATUS_IGNORE);

saxpy<<<(N+255)/256, 256, 0, stream>>>(N, a, d_x, d_y);

cudaMemcpyAsync(y, d_y, N*sizeof(float), cudaMemcpyDeviceToHost, stream);

cudaFree(d_x);

cudaFree(d_y);

free(y);

}

MPI_Comm_free(&stream_comm);

MPIX_Stream_free(&mpi_stream);

cudaStreamDestroy(stream);

MPI_Finalize();

return 0;

}

##  7\. Summary

We have surveyed the current status of MPI+Threads and MPI+GPUs. Both will
need an explicit interface in MPI to allow better arrangements between non-MPI
runtimes and MPI. We proposed a new MPI concept, called MPIX stream, and a set
of new APIs that allow users to communicate their external execution context
to MPI implementations in a general and reliable way. The proposed APIs are
implemented in the MPICH 4.1a1 release. Example codes for typical application
patterns are provided for reference.

###### Acknowledgements.

Special thanks to Jim Dinan from NVIDIA for insightful discussions and his
effort in leading the MPI Forum hybrid working group. We thank MPI developer
teams from Intel corporation and Hewlett Packard Enterprise for feedback and
the members of the MPI Forum for past efforts and on-going discussions related
to this work. We gratefully acknowledge the computing resources provided and
operated by the Joint Laboratory for System Evaluation (JLSE) at Argonne
National Laboratory. This research was supported by the Exascale Computing
Project (17-SC-20-SC), a collaborative effort of the U.S. Department of Energy
Office of Science and the National Nuclear Security Administration, and by the
U.S. Department of Energy, Office of Science, under Contract DE-
AC02-06CH11357.

## References

  * (1)
  * Amer et al. (2015) Abdelhalim Amer, Huiwei Lu, Yanjie Wei, Pavan Balaji, and Satoshi Matsuoka. 2015.  MPI+ threads: Runtime contention and remedies.  _ACM SIGPLAN Notices_ 50, 8 (2015), 239–248.
  * Dinan et al. (2014) James Dinan, Ryan E Grant, Pavan Balaji, David Goodell, Douglas Miller, Marc Snir, and Rajeev Thakur. 2014.  Enabling communication concurrency through flexible MPI endpoints.  _The International Journal of High Performance Computing Applications_ 28, 4 (2014), 390–405.
  * Grant et al. (2019) Ryan E Grant, Matthew GF Dosanjh, Michael J Levenhagen, Ron Brightwell, and Anthony Skjellum. 2019\.  Finepoints: Partitioned multithreaded MPI communication. In _International Conference on High Performance Computing_. Springer, 330–350.
  * Hjelm et al. (2019) Nathan Hjelm, Howard Pritchard, Samuel K. Gutiérrez, Daniel J. Holmes, Ralph Castain, and Anthony Skjellum. 2019.  MPI Sessions: Evaluation of an Implementation in Open MPI. In _2019 IEEE International Conference on Cluster Computing (CLUSTER)_. 1–11.
  * Laboratory (2022) Argonne National Laboratory. 2022\.  MPICH.  https://www.mpich.org/.
  * Message Passing Interface Forum (2021) Message Passing Interface Forum. 2021\.  MPI: A Message-Passing Interface Standard, Version 4.0.  https://www.mpi-forum.org/docs/.
  * Namashivayam et al. (2022) Naveen Namashivayam, Krishna Kandalla, Trey White, Nick Radcliffe, Larry Kaplan, and Mark Pagel. 2022.  Exploring GPU Stream-Aware Message Passing using Triggered Operations.  https://doi.org/10.48550/ARXIV.2208.04817
  * NVIDIA (2022) NVIDIA. 2022.  NCCL.  https://developer.nvidia.com/nccl.
  * Open Fabric Alliance (2022) Open Fabric Alliance. 2022\.  libfabric.  https://github.com/ofiwg/libfabric.
  * Open UCX (2022) Open UCX. 2022\.  UCX.  https://github.com/openucx/ucx.
  * Patinyasakdikul et al. (2019) Thananon Patinyasakdikul, David Eberius, George Bosilca, and Nathan Hjelm. 2019.  Give MPI threading a fair chance: A study of multithreaded MPI designs. In _2019 IEEE International Conference on Cluster Computing (CLUSTER)_. IEEE, 1–11.
  * Team (2022) The Open MPI Team. 2022\.  Open MPI.  https://www.open-mpi.org/.
  * University (2022) The Ohio State University. 2022\.  MVAPICH.  https://mvapich.cse.ohio-state.edu/.
  * Wang and Chandramowlishwaran (2019) Hengjie Wang and Aparna Chandramowlishwaran. 2019.  Multi-criteria partitioning of multi-block structured grids. In _Proceedings of the ACM International Conference on Supercomputing_. 261–271.
  * Wang et al. (2014) Hao Wang, Sreeram Potluri, Devendar Bureddy, Carlos Rosales, and Dhabaleswar K. Panda. 2014\.  GPU-Aware MPI on RDMA-Enabled Clusters: Design, Implementation and Evaluation.  _IEEE Transactions on Parallel and Distributed Systems_ 25, 10 (2014), 2595–2605.  https://doi.org/10.1109/TPDS.2013.222
  * Zambre et al. (2020) Rohit Zambre, Aparna Chandramowliswharan, and Pavan Balaji. 2020\.  How I learned to stop worrying about user-visible endpoints and love MPI. In _Proceedings of the 34th ACM International Conference on Supercomputing_. 1–13.
  * Zambre et al. (2021) Rohit Zambre, Damodar Sahasrabudhe, Hui Zhou, Martin Berzins, Aparna Chandramowlishwaran, and Pavan Balaji. 2021.  Logically Parallel Communication for Fast MPI+ Threads Applications.  _IEEE Transactions on Parallel and Distributed Systems_ 32, 12 (2021), 3038–3052.


# Implications for Cosmic Domain Walls from LIGO-Virgo First Three Observing
Runs

Yang Jiang<EMAIL_ADDRESS>CAS Key Laboratory of Theoretical Physics,
Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190,
China  School of Physical Sciences, University of Chinese Academy of Sciences,
No. 19A Yuquan Road, Beijing 100049, China  Qing-Guo Huang  Corresponding
author<EMAIL_ADDRESS>CAS Key Laboratory of Theoretical Physics, Institute
of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190, China
School of Physical Sciences, University of Chinese Academy of Sciences, No.
19A Yuquan Road, Beijing 100049, China  School of Fundamental Physics and
Mathematical Sciences Hangzhou Institute for Advanced Study, UCAS, Hangzhou
310024, China

###### Abstract

We put constraints on the normalized energy density in gravitaional waves from
cosmic domain walls (DWs) by searching for the stochastic gravitational-wave
background (SGWB) in the data of Advanced LIGO and Virgo’s first three
observing runs. By adopting a phenomenological broken power-law model, we
obtain the upper limit of normalized energy density of SGWB generated by DWs
in the peak frequency band $10\sim 200$ Hz, and get the most stringent
limitation at the peak frequency $f_{*}=35$ Hz, namely
$\Omega_{\text{DW}}(f_{*}=35\,\text{Hz})<1.4\times 10^{-8}$ at $95\%$
confidence level (CL). Subsequently, we work out the constraints on the
parameter space in the appealing realization of DW structure – the heavy axion
model which can avoid the so-called quality problem.

Introduction.  The successful observation of gravitational-wave (GW) signal
from compact binary coalescence Abbott  _et al._ (2021a) by Advanced LIGO Aasi
_et al._ (2015) and Advanced Virgo Acernese  _et al._ (2015) has open up a new
era of GW astronomy. The superposition of a large number of not loud enough GW
incidences can produce a stochastic gravitational-wave background (SGWB) which
encodes information about the physics of astrophysical sources. Until the O3
observing period, the LIGO and Virgo Scientific Collaboration has not detected
such a SGWB, and therefore they set the upper limits of SGWB Abbott  _et al._
(2021b).

Besides the origin of astronomy like compact binary coalescence, SGWB can also
have the cosmological origins, including phase transitions Kibble (1980);
Witten (1984); Hogan (1986); Mazumdar and White (2019), primordial
perturbations during inflation Starobinskiǐ (1979); Turner (1997); Bar-Kana
(1994), cosmic string Kibble (1976); Sarangi and Tye (2002); Damour and
Vilenkin (2005); Siemens  _et al._ (2007); Abbott  _et al._ (2021c) etc.
Detecting SGWB signals from these cosmological origins is of great
significance for understanding the history of the early Universe. In this
letter, we focus on SGWB signal from cosmic domain wall (DW) Vilenkin and
Everett (1982); Sikivie (1982) network – laminar topological defects that form
when a discrete symmetry is spontaneously broken. Such discrete symmetries
arise in various particle physics frameworks beyond the Standard Model (SM),
such as Higgs models Battye  _et al._ (2020); Chacko  _et al._ (2006); Di
Luzio  _et al._ (2019), grand unification Lazarides  _et al._ (1982); Everett
and Vilenkin (1982), supersymmetry Abel  _et al._ (1995); Dvali and Shifman
(1997); Kovner  _et al._ (1997), non-Abelian discrete symmetries Riva (2010),
axion models Sikivie (1982); Grilli di Cortona  _et al._ (2016) and so on.
Once the breaking of symmetry occurs after inflation, DWs begin to form and
leave different Hubble patches in different degenerate vacua. Although the
formation of DW networks cause potential conflict with cosmological
observations if they dominate the energy and then overclose the universe
Zel’dovich  _et al._ (1974), one can expect that the discrete symmetry is only
approximate Larsson  _et al._ (1997) and the energy bias of the potential can
induce the annihilation of the DWs. Along with the motion and annihilation of
DWs, GWs are radiated and form the SGWB for present observers. Through the
observation and investigation of these signals, it is possible to deduce
constraints on related theoretical models and explore the physics at the
energy scales inaccessible by accelerator laboratories.

An appealing mechanism forming such a DW network is the spontaneous breaking
of $U(1)$ Peccei-Quinn (PQ) symmetry Peccei and Quinn (1977a, b). The
resulting QCD axion Weinberg (1978) explains the undetectable CP violation in
strong interaction. However, PQ mechanism depends particularly on the axion
potential from QCD instantons and it is susceptible to additional global
symmetry breaking. It is widely accepted that there is no continuous global
symmetry in quantum gravity Banks and Seiberg (2011) and thus PQ solution is
spoiled. This so-called quality problem has always been perplexing. To avoid
the problem, heavy axion Holdom and Peskin (1982); Holdom (1985); Rubakov
(1997); Dimopoulos  _et al._ (2016); Gherghetta and Nguyen (2020); Zambujal
Ferreira  _et al._ (2022) is taken to be consideration. Its heavy sector
provides a larger contribution to the axion potential and lift the quality of
PQ symmetry. In general, the heavy axion decays more easily and leave few
detectable objects. Therefore, detecting GW radiation becomes an important way
to understand this model.

In this letter, we provide a first search for the SGWB from DWs using the data
of Advanced LIGO and Advanced Virgo’s first three observing runs. First of
all, we take a model-independent analysis and figure out the constraints on
the parameters characterizing the SGWB from DWs. Then, like Zambujal Ferreira
_et al._ (2022), we focus on the heavy axion Holdom and Peskin (1982); Holdom
(1985) for explaining strong CP problems and avoiding the vulnerable PQ
mechanism. In the standard QCD axion model, the DW network makes up a very
tiny fraction of energy density of the Universe and leave SGWB undetectable
Hiramatsu  _et al._ (2013). But GWs from the DW network in the heavy axion
model can be much stronger because of the large tension carried by heavy axion
DWs.

Recently, based on the evidence for a stochastic common-spectrum process in
Pulsar Timing Arrays (PTAs) dataset Arzoumanian  _et al._ (2020); Chen  _et
al._ (2021a, b); Goncharov  _et al._ (2021); Chen  _et al._ (2021c);
Antoniadis  _et al._ (2022), the energy spectrum of SGWB from DWs has been
searched and used to put constraint on axionlike particles in Ferreira  _et
al._ (2022). In our analysis, instead, the data from Advanced LIGO and
Advanced Virgo’s first three observing runs is adopted to probe the mechanics
of DWs at much higher energy scale because the frequency band of LIGO and
Virgo is much higher than that of PTA.

SGWB from DWs.  The energy spectrum of SGWB normalized by the critical density
of Universe is defined as

$\Omega_{\text{gw}}(f)=\frac{1}{\rho_{\text{c}}}\frac{\mathrm{d}\rho_{\text{gw}}}{\mathrm{d}\ln
f},$ (1)

where $\rho_{\text{c}}$ is present critical density, $\rho_{\text{gw}}$ is the
energy density of GWs, and $f$ is the frequency of GWs. In the absence of
significant friction from surrounding plasma, the DW network quickly achieves
a scaling regime Kibble (1976, 1980) with the energy density
$\rho_{\text{DW}}=c\sigma H$, where $c\sim\mathcal{O}(1)$ is a model-dependent
prefactor, $\sigma$ is the tension or surface energy density of DW, and $H$ is
the Hubble rate. This corresponds to the fractional energy density:

$\alpha=\frac{\rho_{\text{DW}}}{3H^{2}M_{p}^{2}},$ (2)

where $M_{p}=(8\pi G)^{-1/2}$ is the reduced Planck energy scale. Notice that
DWs can possess large energy and become dangerous to cosmology Zel’dovich  _et
al._ (1974); Zel’Dovich  _et al._ (1975). Therefore, a mechanism inducing the
annihilation of DW network at certain time is usually needed in practice.

Due to the time-varying quadrupole, DW network emits GWs and the simple
estimations show that $\rho_{\text{gw}}\sim\rho_{\text{DW}}^{2}/(32\pi
H^{2}M_{p}^{2})$ Vilenkin (1981); Kawasaki and Saikawa (2011); Hiramatsu  _et
al._ (2014); Preskill  _et al._ (1991); Chang  _et al._ (1999); Gleiser and
Roberts (1998). Most GWs are radiated at the frequency $\tilde{f}_{*}\simeq H$
corresponding to the time of annihilation and then redshifted to today with
the peak frequency:

$f_{*}\simeq
1.5\left(\frac{g_{*}}{106.75}\right)^{1/6}\left(\frac{T_{*}}{10^{7}\,\text{GeV}}\right)\,\text{Hz},$
(3)

where $g_{*}$ is the number of relativistic degrees of freedom at
annihilation. We take the number of (entropy) relativistic degrees of freedom
$g_{*,s}=g_{*}=106.75$ in the whole letter. By assuming annihilation happens
at radiation domination period, numerical simulations Hiramatsu  _et al._
(2013, 2014) show that the normalized energy density of SGWB takes the form

$\Omega_{\text{DW}}(f)h^{2}\simeq 4.7\times
10^{-7}\tilde{\epsilon}\alpha_{*}^{2}\left(\frac{106.75}{g_{*}}\right)^{1/3}\mathcal{S}\left(\frac{f}{f_{*}}\right).$
(4)

Here $\alpha_{*}$ to represent the fractional energy density at annihilation.
$\tilde{\epsilon}\simeq 0.1-1$ is a numerical efficiency factor which is
suggested in Hiramatsu  _et al._ (2014), and its value is fixed to be $0.7$ in
our analysis. The causality ensures $\Omega_{\text{DW}}\sim f^{3}$ when
$f<f_{*}$. Besides, $\Omega_{\text{DW}}\sim f^{-1}$ when $f\to\infty$ until
the cutoff depicted by the inverse of wall width. The order of
$\mathcal{O}(1)$ width around the maximum is suggested by numerical studies
Hiramatsu  _et al._ (2014). Therefore, the shape of SGWB spectrum can be
parameterized as a broken power law:

$\mathcal{S}(x)=\frac{4}{x^{-3}+3x}.$ (5)

We adopt cross-correlation spectrum of isotropic SGWB during advanced LIGO’s
and advanced Virgo’s O1$\sim$O3 observing runs Abbott  _et al._ (2021b) and
follow the method described in Mandic  _et al._ (2012); Romano and Cornish
(2017) to do Bayesian analysis and estimate the probability of models. The
likelihood is given by

$p(\hat{C}|\bm{\theta};\lambda)\propto\exp\left[-\frac{1}{2}\sum_{IJ}\sum_{f}\frac{\left(\hat{C}_{IJ}(f)-\Omega_{\text{gw}}(f;\bm{\theta})\right)^{2}}{\sigma^{2}_{IJ}(f)}\right].$
(6)

The model of energy spectrum $\Omega_{\text{gw}}(f;\bm{\theta})$ is
parameterized by a series parameters $\bm{\theta}$. $\hat{C}_{IJ}$ is the
cross-correlation spectrum and $\sigma_{IJ}$ denotes its variance. Likelihoods
obtained by each detector pair are multiplied together to obtain Eq. (6). The
Bayes factor between SGWB signal and pure noise is used to show the fitness of
the model. Here $\bm{\theta}=(\alpha_{*},f_{*})$ are free parameters and their
priors are listed in Table 1.

Parameter |  Prior
---|---
$\alpha_{*}$ |  $\text{Uniform}\>[0,\,0.3]$
$f_{*}\,\text{(Hz)}$ |  $\text{LogUniform}\>[10,\,200]$
Table 1: Priors of the parameters adopted in broken power law model.

The Bayesian analysis makes use of python bilby package Ashton  _et al._
(2019) and dynamic nested sampling package dynesty Koposov  _et al._ (2022).
The posterior distributions of parameters in the broken power law model is
presented in Fig. 1. The Bayes factor is
$\log\mathcal{B}_{\text{noise}}^{\text{DW}}=-0.27$, indicating that no
evidence of such signal in the strain data. In fact, the posterior of
$\alpha_{*}$ allows us to put constraint on the amplitude of the energy
spectrum of GWs from DW network. In Fig. 2, the upper limit of
$\Omega_{\text{DW}}(f_{*})$ at $95\%$ CL is illustrated for different peak
frequency $f_{*}$. The most stringent limitation is
$\Omega_{\text{DW}}(f_{*}=35\,\text{Hz})<1.4\times 10^{-8}$ at $95\%$ CL.

Figure 1: Posterior distribution of $\alpha_{*}$ and $f_{*}$ for the broken
power law model. The $68\%$ and $95\%$ CL. exclusion contours are shown with
blue shaded region. Figure 2: Upper limit of $\Omega_{\text{DW}}(f_{*})$ at
95% CL for different peak frequency $f_{*}$.

After the annihilation, DWs may decay into SM particles or dark radiation
(DR). Usually the abundance of DR is described by the effective number of
neutrino species $\Delta N_{\text{eff}}:=\rho_{\text{DR}}/\rho_{\nu}$. If DWs
decay into DR entirely, $\rho_{DR}\simeq\rho_{\text{DW}}$ and then

$\Delta N_{\text{eff,decay}}\simeq 13.6g_{*}^{-1/3}\alpha_{*},$ (7)

In addition, GWs emitted by the DW network also make a contribution to DR,
namely

$\Delta N_{\mathrm{eff,gw}}\simeq
0.09\alpha_{*}^{2}\left(\frac{g_{*}}{106.75}\right)^{-1/3}.$ (8)

In this scenario, $\Delta N_{\mathrm{eff}}=\Delta
N_{\mathrm{eff,decay}}+\Delta N_{\mathrm{eff,gw}}$. It is worth mentioning
that BBN (CMB+BAO) also put constraints on the abundance of DR with $\Delta
N_{\text{eff}}\leq 0.39$ Fields  _et al._ (2020) ($0.29$ Aghanim  _et al._
(2020)) at $95\%$ confidence level (CL). This may also provide a constraint on
the energy density of DW network. In Fig. 3, we transform the posterior of
$(\alpha_{*},f_{*})$ to $(\Delta N_{\mathrm{eff}},T_{*})$ using Eqs. (3), (7)
and (8). It implies that the constraints from LIGO-Virgo is comparable to
those from BBN and CMB+BAO in the scenario where DWs are all assumed to decay
into DR.

Figure 3: Posterior distribution of $\Delta N_{\text{eff}}$ and $T_{*}$ (GeV)
for the broken power law model. The $68\%$ and $95\%$ CL. exclusion contours
are shown with blue shaded region. The vertical lines denote the constraint
from BBN (CMB+BAO).

Heavy QCD Axion Model.  From now on, we will use the former results to
constrain the parameter space in the heavy axion model. The composite
potential of heavy axion periodic potential and the contribution from QCD is
given by

$V_{a}=\left(\kappa_{\text{Q}}^{2}\Lambda_{\text{Q}}^{4}+\kappa_{\text{H}}^{2}\Lambda_{\text{H}}^{4}\right)\left[1-\cos\left(\frac{a}{F_{a}}\right)\right],$
(9)

where $F_{a}$ is axion decay constant, $\Lambda_{\text{Q,H}}$ are coupling
scales of QCD and heavy sector. The factors $\kappa_{\text{Q,H}}\leq 1$ are
related to mass of fermion. After the spontaneous breaking of $U(1)$ PQ
symmetry, axionic string-DW network forms and then GWs are emitted. In the
form of Eq. (9), periodic potential from the heavy sector is aligned with that
of QCD. It is typically ensured by the $\mathbb{Z}_{2}$ symmetry Hook  _et
al._ (2020). To that end, QCD cannot induce the annihilation of DWs. However,
regardless of its specific physical origin (See some models in Dvali (2005);
Dine and Seiberg (1986); Svrcek and Witten (2006); Kamionkowski and March-
Russell (1992)), we could expect a term misaligned with $V_{a}$ as follows

$V_{b}\simeq-\mu_{b}^{4}\cos\left(\frac{N_{b}}{N_{\text{DW}}}\frac{a}{F_{a}}-\delta\right),$
(10)

where $\mu_{b}$ is related to the specific model under consideration, $N_{b},\
N_{\text{DW}}$ are integers describing discrete symmetry subgroup. Note that
energy bias of potential only appears when $N_{b}=1$ or co-prime with
$N_{\text{DW}}$. $\delta$ is the generic misaligned phase. CP violation caused
by $V_{b}$ is depicted by

$\Delta\theta=\theta-\theta_{\text{Q}}\simeq
r^{4}\left(\frac{N_{b}}{N_{\text{DW}}}\right)\left(\frac{\sin\delta}{\kappa_{H}^{2}}\right),$
(11)

where $r=\mu_{b}/\Lambda_{\text{H}}$. Once $N_{\text{DW}}>1$, a long-lasting
DW network could form, and the fractional energy density is

$\alpha_{*}\simeq
0.1\left(\frac{\sqrt{d}c\kappa_{\text{H}}}{\sin(N_{b}\pi/N_{\text{DW}})}\right)^{2}\left(\frac{F_{a}}{10^{12}\,\text{GeV}}\right)^{2}\left(\frac{0.002}{r}\right)^{4},$
(12)

and the temperature at annihilation reads

$T_{*}\simeq
10^{8}\frac{\sin(N_{b}\pi/N_{\text{DW}})}{\sqrt{dc\kappa_{\text{H}}}}\left(\frac{106.75}{g_{*}}\right)^{1/4}\left(\frac{r}{0.005}\right)^{2}\\\
\left(\frac{10^{12}\,\text{GeV}}{F_{a}}\right)^{1/2}\left(\frac{\Lambda_{\text{H}}}{10^{10}\,\text{GeV}}\right)\,\text{GeV}.$
(13)

Different from the standard QCD axion, the heavy axion tends to decay into SM
particles: gluons, photons or fermions, depending on the energy scale, and
then the heavy axion avoids the constraints on the abundance of DR from BBN
and CMB+BAO. If the heavy axion decays not fast enough, there is a temporary
period of matter domination (MD) and hence Eq. (3) and Eq. (4) should be
modified. To evaluate whether the situation takes place, we compare the
characteristic temperature of efficient decay with that of MD, where,
according to Zambujal Ferreira  _et al._ (2022), the effective characteristic
temperature of decay is

$T_{\text{d}}\simeq
10^{6}\left(\frac{\sqrt{\kappa_{\text{H}}}\Lambda_{\text{H}}}{10^{10}\,\text{GeV}}\right)^{3}\left(\frac{10^{12}\,\text{GeV}}{F_{a}}\right)^{5/2}\text{GeV},$
(14)

and the temperature of MD is

$T_{\text{MD}}=8\times
10^{6}\left(\frac{\sqrt{d}(c\kappa_{\text{H}})^{3/2}}{\sin(N_{b}\pi/N_{\text{DW}})}\right)\left(\frac{F_{a}}{10^{12}\,\text{GeV}}\right)^{3/2}\\\
\left(\frac{\Lambda_{\text{H}}}{10^{10}\,\text{GeV}}\right)\left(\frac{0.001}{r}\right)^{2}\left(\frac{106.75}{g_{*}}\right)^{1/4}\text{GeV}.$
(15)

If $T_{\text{d}}<T_{\text{MD}}$, there is a MD phase, and the peak frequency
and the energy spectrum of SGWB from DWs become Zambujal Ferreira  _et al._
(2022),

$\displaystyle f_{*}\to
f_{*}\left(\frac{T_{\text{d}}}{T_{\text{MD}}}\right)^{1/3},$ (16)
$\displaystyle\Omega_{\text{DW}}\to\Omega_{\text{DW}}\left(\frac{g_{*}(T_{\text{d}})}{g_{*}(T_{\text{MD}})}\right)^{1/3}\left(\frac{T_{\text{d}}}{T_{\text{MD}}}\right)^{4/3}.$
(17)

In fact, there are also some inherent limitations which should be taken into
account in our analysis. Firstly, the characteristic temperature
$(T_{\text{DW-dom}})$ of string-wall dominating the Universe can be estimated
by $\alpha_{*}\simeq 1$ and

$T_{\text{DW-dom}}\simeq 5\times
10^{6}\sqrt{\frac{c\kappa_{\text{H}}F_{a}}{10^{12}\,\text{GeV}}}\\\
\left(\frac{\Lambda_{\text{H}}}{10^{10}\,\text{GeV}}\right)\left(\frac{106.75}{g_{*}}\right)^{1/4}\>\text{GeV}.$
(18)

In the case of $T_{\text{DW-dom}}>T_{*}$, DWs will dominate our Universe,
which may conflict with the standard cosmology. Even if the contradiction
might be avoided Kawasaki and Takahashi (2005), there is a lack of knowledge
about the evolution of DWs in the situation. Secondly, the validity of the
axion effective theory requires
$\sqrt{\kappa_{\text{H}}}\Lambda_{\text{H}}<F_{a}$. At last, Neutron Electric
Dipole Moment (nEDM) measurements Abel  _et al._ (2020) require
$\Delta\theta\lesssim 10^{-10}$.

In this letter we mainly focus on constraining the three physical parameters
(e.g. the coupling scale $\Lambda_{\text{H}}$ of heavy sector, axion decay
constant $F_{a}$ and CP violation coefficient $\Delta\theta$) in the heavy
axion model. From Eq. (12) and Eq. (13), the former constraints on
$\alpha_{*}$ and $T_{*}$ can be used to constrain these three parameters. Here
we set the factor $\kappa_{\text{H}}\simeq 1$ providing that the heavy sector
contains a light fermion. Similar to Kawasaki  _et al._ (2015), we take
$N_{b}=1,N_{\text{DW}}=6$, and then $c=4.48,\ d=6.28$. The dependence of
$\delta$ only appears in Eq. (11) as a factor between $\Delta\theta$ and $r$.
Therefore, without losing generality, for example, we set $\delta=0.3$. Since
there are three free parameters, our strategy of analysis is to constrain the
other two parameters by keeping $\Lambda_{\text{H}}$ or $\Delta\theta$ fixed,
respectively. The priors adopted in our analysis are listed in Table 2, and
the excluded parameter spaces in the heavy axion model are illustrated by the
blue shaded regions in Fig. 4.

$\Lambda_{\text{H}}=10^{10},10^{10.5},10^{11}$ GeV
---
Parameter |  Prior
$F_{a}\,\text{(GeV)}$ |  $\text{LogUniform}\>[10^{10},\,10^{12.2}]$
$\Delta\theta$ |  $\text{LogUniform}\>[10^{-14.5},\,10^{-10}]$
$\Delta\theta=10^{-13},10^{-12},10^{-11}$
Parameter |  Prior
$F_{a}\,\text{(GeV)}$ |  $\text{LogUniform}\>[10^{10},\,10^{12.2}]$
$\Lambda_{\text{H}}\,\text{(GeV)}$ |  $\text{LogUniform}\>[10^{9},\,10^{12.2}]$
Table 2: Priors of the parameters adopted in the heavy axion model.

Figure 4: Parameter space excluded by Advanced LIGO and Advanced Virgo
O1$\sim$O3 observing runs. The blue shaded regions denote the 95% exclusion
contours. The orange-red and orange shaded regions correspond to DW domination
and violation of axion effective theory, respectively.

Conclusions. By adopting the cross-correlation spectrum of Advanced LIGO and
Advanced Virgo’s first three observing runs, we search for the GW signal from
DW network in the peak frequency range of $10\sim 200$ Hz. Since no such
signal in the strain data has been observed, we place constraints on the
energy spectrum of SGWB from cosmic DW network. The most stringent limitation
is $\Omega_{\text{DW}}(f_{*}=35\,\text{Hz})<1.4\times 10^{-8}$ at $95\%$ CL.
The posterior of $f_{*}$ dose not show obvious preference for any certain peak
frequency, and the fractional energy density of DWs with $\alpha_{*}>0.2$ is
excluded at $95\%$ CL for $f_{*}=20\sim 100$ Hz corresponding to the most
sensitive band of LIGO and Virgo.

Furthermore, for fixed values of $\Lambda_{\text{H}}=10^{10}$, $10^{10.5}$ or
$10^{11}$ GeV, the CP violation coefficient $\Delta\theta$ in the blue shaded
regions in the upper panel of Fig. 4 is excluded; for fixed values of
$\Delta\theta=10^{-13}$, $10^{-12}$ or $10^{-11}$, the coupling scale
$\Lambda_{\text{H}}$ of the heavy sector in the blue shaded regions in the
bottom panel of Fig. 4 is excluded. In all, even though SGWB has not been
detected by Advanced LIGO and Advanced Virgo, the non-detection of the SGWB
from cosmic DWs already places constraints on the heavy axion model.

Our results have shown the ability of LVK interferometers to understand
physics related to DWs and we expect more elaborate results with the
improvement of the detector sensitivity in the near future.

Acknowledgments.  We acknowledge the use of HPC Cluster of ITP-CAS and HPC
Cluster of Tianhe II in National Supercomputing Center in Guangzhou. This work
is supported by the National Key Research and Development Program of China
Grant No.2020YFC2201502, grants from NSFC (grant No. 11975019, 11991052,
12047503), Key Research Program of Frontier Sciences, CAS, Grant NO. ZDBS-
LY-7009, CAS Project for Young Scientists in Basic Research YSBR-006, the Key
Research Program of the Chinese Academy of Sciences (Grant NO. XDPB15).

## References

  * Abbott  _et al._ (2021a) R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), “GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run,” (2021a), arXiv:2111.03606 [gr-qc] .
  * Aasi  _et al._ (2015) J. Aasi _et al._ (LIGO Scientific), “Advanced LIGO,” Class. Quant. Grav. 32, 074001 (2015), arXiv:1411.4547 [gr-qc] .
  * Acernese  _et al._ (2015) F. Acernese _et al._ (VIRGO), “Advanced Virgo: a second-generation interferometric gravitational wave detector,” Class. Quant. Grav. 32, 024001 (2015), arXiv:1408.3978 [gr-qc] .
  * Abbott  _et al._ (2021b) R. Abbott _et al._ (KAGRA, Virgo, LIGO Scientific), “Upper limits on the isotropic gravitational-wave background from Advanced LIGO and Advanced Virgo’s third observing run,” Phys. Rev. D 104, 022004 (2021b), arXiv:2101.12130 [gr-qc] .
  * Kibble (1980) T.W.B. Kibble, “Some implications of a cosmological phase transition,” Physics Reports 67, 183–199 (1980).
  * Witten (1984) Edward Witten, “Cosmic separation of phases,” Phys. Rev. D 30, 272–285 (1984).
  * Hogan (1986) C. J. Hogan, “Gravitational radiation from cosmological phase transitions,” Monthly Notices of the Royal Astronomical Society 218, 629–636 (1986), https://academic.oup.com/mnras/article-pdf/218/4/629/3299141/mnras218-0629.pdf .
  * Mazumdar and White (2019) Anupam Mazumdar and Graham White, “Review of cosmic phase transitions: their significance and experimental signatures,” Rept. Prog. Phys. 82, 076901 (2019), arXiv:1811.01948 [hep-ph] .
  * Starobinskiǐ (1979) A. A. Starobinskiǐ, “Spectrum of relict gravitational radiation and the early state of the universe,” Soviet Journal of Experimental and Theoretical Physics Letters 30, 682 (1979).
  * Turner (1997) Michael S. Turner, “Detectability of inflation-produced gravitational waves,” Phys. Rev. D 55, R435–R439 (1997).
  * Bar-Kana (1994) Rennan Bar-Kana, “Limits on direct detection of gravitational waves,” Phys. Rev. D 50, 1157–1160 (1994).
  * Kibble (1976) T W B Kibble, “Topology of cosmic domains and strings,” Journal of Physics A: Mathematical and General 9, 1387–1398 (1976).
  * Sarangi and Tye (2002) Saswat Sarangi and S.-H.Henry Tye, “Cosmic string production towards the end of brane inflation,” Physics Letters B 536, 185–192 (2002).
  * Damour and Vilenkin (2005) Thibault Damour and Alexander Vilenkin, “Gravitational radiation from cosmic (super)strings: Bursts, stochastic background, and observational windows,” Phys. Rev. D 71, 063510 (2005).
  * Siemens  _et al._ (2007) Xavier Siemens, Vuk Mandic, and Jolien Creighton, “Gravitational-wave stochastic background from cosmic strings,” Phys. Rev. Lett. 98, 111101 (2007).
  * Abbott  _et al._ (2021c) R. Abbott _et al._ (LIGO Scientific, Virgo, KAGRA), “Constraints on Cosmic Strings Using Data from the Third Advanced LIGO–Virgo Observing Run,” Phys. Rev. Lett. 126, 241102 (2021c), arXiv:2101.12248 [gr-qc] .
  * Vilenkin and Everett (1982) Alexander Vilenkin and Allen E. Everett, “Cosmic Strings and Domain Walls in Models with Goldstone and Pseudo-Goldstone Bosons,” Phys. Rev. Lett. 48, 1867–1870 (1982).
  * Sikivie (1982) P. Sikivie, “Axions, domain walls, and the early universe,” Phys. Rev. Lett. 48, 1156–1159 (1982).
  * Battye  _et al._ (2020) Richard A. Battye, Apostolos Pilaftsis, and Dominic G. Viatic, “Domain wall constraints on two-Higgs-doublet models with $Z_{2}$ symmetry,” Phys. Rev. D 102, 123536 (2020), arXiv:2010.09840 [hep-ph] .
  * Chacko  _et al._ (2006) Z. Chacko, Hock-Seng Goh, and Roni Harnik, “The Twin Higgs: Natural electroweak breaking from mirror symmetry,” Phys. Rev. Lett. 96, 231802 (2006), arXiv:hep-ph/0506256 .
  * Di Luzio  _et al._ (2019) Luca Di Luzio, Michele Redi, Alessandro Strumia, and Daniele Teresi, “Coset Cosmology,” JHEP 06, 110 (2019), arXiv:1902.05933 [hep-ph] .
  * Lazarides  _et al._ (1982) G. Lazarides, Q. Shafi, and T.F. Walsh, “Cosmic strings and domains in unified theories,” Nuclear Physics B 195, 157–172 (1982).
  * Everett and Vilenkin (1982) Allen E. Everett and Alexander Vilenkin, “Left-right symmetric theories and vacuum domain walls and strings,” Nuclear Physics B 207, 43–53 (1982).
  * Abel  _et al._ (1995) S. A. Abel, Subir Sarkar, and P. L. White, “On the cosmological domain wall problem for the minimally extended supersymmetric standard model,” Nucl. Phys. B 454, 663–684 (1995), arXiv:hep-ph/9506359 .
  * Dvali and Shifman (1997) G. R. Dvali and Mikhail A. Shifman, “Domain walls in strongly coupled theories,” Phys. Lett. B 396, 64–69 (1997), [Erratum: Phys.Lett.B 407, 452 (1997)], arXiv:hep-th/9612128 .
  * Kovner  _et al._ (1997) A. Kovner, Mikhail A. Shifman, and Andrei V. Smilga, “Domain walls in supersymmetric Yang-Mills theories,” Phys. Rev. D 56, 7978–7989 (1997), arXiv:hep-th/9706089 .
  * Riva (2010) Francesco Riva, “Low-Scale Leptogenesis and the Domain Wall Problem in Models with Discrete Flavor Symmetries,” Phys. Lett. B 690, 443–450 (2010), arXiv:1004.1177 [hep-ph] .
  * Grilli di Cortona  _et al._ (2016) Giovanni Grilli di Cortona, Edward Hardy, Javier Pardo Vega, and Giovanni Villadoro, “The QCD axion, precisely,” JHEP 01, 034 (2016), arXiv:1511.02867 [hep-ph] .
  * Zel’dovich  _et al._ (1974) Ia.B. Zel’dovich, I. Iu. Kobzarev, and Lev Borisovich Okun, “Cosmological consequences of a spontaneous breakdown of a discrete symmetry,” Journal of Experimental and Theoretical Physics (1974).
  * Larsson  _et al._ (1997) Sebastian E. Larsson, Subir Sarkar, and Peter L. White, “Evading the cosmological domain wall problem,” Phys. Rev. D 55, 5129–5135 (1997), arXiv:hep-ph/9608319 .
  * Peccei and Quinn (1977a) R.D. Peccei and Helen Quinn, “Cp conservation in the presence of pseudoparticles,” Physical Review Letters - PHYS REV LETT 38, 1440–1443 (1977a).
  * Peccei and Quinn (1977b) Roberto D. Peccei and Helen R. Quinn, “Constraints imposed by cp conservation in the presence of pseudoparticles,” Physical Review D 16, 1791–1797 (1977b).
  * Weinberg (1978) Steven Weinberg, “A new light boson?” Phys. Rev. Lett. 40, 223–226 (1978).
  * Banks and Seiberg (2011) Tom Banks and Nathan Seiberg, “Symmetries and Strings in Field Theory and Gravity,” Phys. Rev. D 83, 084019 (2011), arXiv:1011.5120 [hep-th] .
  * Holdom and Peskin (1982) Bob Holdom and Michael E. Peskin, “Raising the axion mass,” Nuclear Physics B 208, 397–412 (1982).
  * Holdom (1985) Bob Holdom, “Strong qcd at high energies and a heavy axion,” Physics Letters B 154, 316–319 (1985).
  * Rubakov (1997) V. A. Rubakov, “Grand unification and heavy axion,” JETP Lett. 65, 621–624 (1997), arXiv:hep-ph/9703409 .
  * Dimopoulos  _et al._ (2016) Savas Dimopoulos, Anson Hook, Junwu Huang, and Gustavo Marques-Tavares, “A collider observable QCD axion,” JHEP 11, 052 (2016), arXiv:1606.03097 [hep-ph] .
  * Gherghetta and Nguyen (2020) Tony Gherghetta and Minh D. Nguyen, “A Composite Higgs with a Heavy Composite Axion,” JHEP 12, 094 (2020), arXiv:2007.10875 [hep-ph] .
  * Zambujal Ferreira  _et al._ (2022) Ricardo Zambujal Ferreira, Alessio Notari, Oriol Pujolàs, and Fabrizio Rompineve, “High Quality QCD Axion at Gravitational Wave Observatories,” Phys. Rev. Lett. 128, 141101 (2022), arXiv:2107.07542 [hep-ph] .
  * Hiramatsu  _et al._ (2013) Takashi Hiramatsu, Masahiro Kawasaki, Ken’ichi Saikawa, and Toyokazu Sekiguchi, “Axion cosmology with long-lived domain walls,” JCAP 01, 001 (2013), arXiv:1207.3166 [hep-ph] .
  * Arzoumanian  _et al._ (2020) Zaven Arzoumanian _et al._ (NANOGrav), “The NANOGrav 12.5 yr Data Set: Search for an Isotropic Stochastic Gravitational-wave Background,” Astrophys. J. Lett. 905, L34 (2020), arXiv:2009.04496 [astro-ph.HE] .
  * Chen  _et al._ (2021a) Zu-Cheng Chen, Yu-Mei Wu, and Qing-Guo Huang, “Searching for Isotropic Stochastic Gravitational-Wave Background in the International Pulsar Timing Array Second Data Release,” (2021a), arXiv:2109.00296 [astro-ph.CO] .
  * Chen  _et al._ (2021b) Zu-Cheng Chen, Chen Yuan, and Qing-Guo Huang, “Non-tensorial gravitational wave background in NANOGrav 12.5-year data set,” Sci. China Phys. Mech. Astron. 64, 120412 (2021b), arXiv:2101.06869 [astro-ph.CO] .
  * Goncharov  _et al._ (2021) Boris Goncharov _et al._ , “On the Evidence for a Common-spectrum Process in the Search for the Nanohertz Gravitational-wave Background with the Parkes Pulsar Timing Array,” Astrophys. J. Lett. 917, L19 (2021), arXiv:2107.12112 [astro-ph.HE] .
  * Chen  _et al._ (2021c) S. Chen _et al._ , “Common-red-signal analysis with 24-yr high-precision timing of the European Pulsar Timing Array: inferences in the stochastic gravitational-wave background search,” Mon. Not. Roy. Astron. Soc. 508, 4970–4993 (2021c), arXiv:2110.13184 [astro-ph.HE] .
  * Antoniadis  _et al._ (2022) J. Antoniadis _et al._ , “The International Pulsar Timing Array second data release: Search for an isotropic gravitational wave background,” Mon. Not. Roy. Astron. Soc. 510, 4873–4887 (2022), arXiv:2201.03980 [astro-ph.HE] .
  * Ferreira  _et al._ (2022) Ricardo Z. Ferreira, Alessio Notari, Oriol Pujolas, and Fabrizio Rompineve, “Gravitational Waves from Domain Walls in Pulsar Timing Array Datasets,” (2022), arXiv:2204.04228 [astro-ph.CO] .
  * Zel’Dovich  _et al._ (1975) Ya. B. Zel’Dovich, I. Yu. Kobzarev, and L. B. Okun’, “Cosmological consequences of a spontaneous breakdown of a discrete symmetry,” Soviet Journal of Experimental and Theoretical Physics 40, 1 (1975).
  * Vilenkin (1981) Alexander Vilenkin, “Gravitational field of vacuum domain walls and strings,” Phys. Rev. D 23, 852–857 (1981).
  * Kawasaki and Saikawa (2011) Masahiro Kawasaki and Ken’ichi Saikawa, “Study of gravitational radiation from cosmic domain walls,” JCAP 09, 008 (2011), arXiv:1102.5628 [astro-ph.CO] .
  * Hiramatsu  _et al._ (2014) Takashi Hiramatsu, Masahiro Kawasaki, and Ken’ichi Saikawa, “On the estimation of gravitational wave spectrum from cosmic domain walls,” JCAP 02, 031 (2014), arXiv:1309.5001 [astro-ph.CO] .
  * Preskill  _et al._ (1991) John Preskill, Sandip P. Trivedi, Frank Wilczek, and Mark B. Wise, “Cosmology and broken discrete symmetry,” Nuclear Physics B 363, 207–220 (1991).
  * Chang  _et al._ (1999) Sanghyeon Chang, C. Hagmann, and P. Sikivie, “Studies of the motion and decay of axion walls bounded by strings,” Phys. Rev. D 59, 023505 (1999), arXiv:hep-ph/9807374 .
  * Gleiser and Roberts (1998) Marcelo Gleiser and Ronald Roberts, “Gravitational waves from collapsing vacuum domains,” Phys. Rev. Lett. 81, 5497–5500 (1998), arXiv:astro-ph/9807260 .
  * Mandic  _et al._ (2012) V. Mandic, E. Thrane, S. Giampanis, and T. Regimbau, “Parameter estimation in searches for the stochastic gravitational-wave background,” Phys. Rev. Lett. 109, 171102 (2012).
  * Romano and Cornish (2017) Joseph D. Romano and Neil J. Cornish, “Detection methods for stochastic gravitational-wave backgrounds: a unified treatment,” Living Rev. Rel. 20, 2 (2017), arXiv:1608.06889 [gr-qc] .
  * Ashton  _et al._ (2019) Gregory Ashton, Moritz Hübner, Paul D. Lasky, Colm Talbot, Kendall Ackley, Sylvia Biscoveanu, Qi Chu, Atul Divakarla, Paul J. Easter, Boris Goncharov, Francisco Hernandez Vivanco, Jan Harms, Marcus E. Lower, Grant D. Meadors, Denyz Melchor, Ethan Payne, Matthew D. Pitkin, Jade Powell, Nikhil Sarin, Rory J. E. Smith, and Eric Thrane, “Bilby: A user-friendly bayesian inference library for gravitational-wave astronomy,” The Astrophysical Journal Supplement Series 241, 27 (2019).
  * Koposov  _et al._ (2022) Sergey Koposov, Josh Speagle, Kyle Barbary, Gregory Ashton, Johannes Buchner, Carl Scheffler, Ben Cook, Colm Talbot, James Guillochon, Patricio Cubillos, Andrés Asensio Ramos, Ben Johnson, Dustin Lang, Ilya, Matthieu Dartiailh, Alex Nitz, Andrew McCluskey, Anne Archibald, Christoph Deil, Dan Foreman-Mackey, Danny Goldstein, Erik Tollerud, Joel Leja, Matthew Kirk, Matt Pitkin, Patrick Sheehan, Phillip Cargile, ruskin23, Ruth Angus, and Tansu Daylan, “joshspeagle/dynesty: v1.2.3,”  (2022).
  * Fields  _et al._ (2020) Brian D. Fields, Keith A. Olive, Tsung-Han Yeh, and Charles Young, “Big-Bang Nucleosynthesis after Planck,” JCAP 03, 010 (2020), [Erratum: JCAP 11, E02 (2020)], arXiv:1912.01132 [astro-ph.CO] .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Hook  _et al._ (2020) Anson Hook, Soubhik Kumar, Zhen Liu, and Raman Sundrum, “High Quality QCD Axion and the LHC,” Phys. Rev. Lett. 124, 221801 (2020), arXiv:1911.12364 [hep-ph] .
  * Dvali (2005) Gia Dvali, “Three-form gauging of axion symmetries and gravity,” (2005), arXiv:hep-th/0507215 .
  * Dine and Seiberg (1986) Michael Dine and Nathan Seiberg, “String theory and the strong cp problem,” Nuclear Physics B 273, 109–124 (1986).
  * Svrcek and Witten (2006) Peter Svrcek and Edward Witten, “Axions In String Theory,” JHEP 06, 051 (2006), arXiv:hep-th/0605206 .
  * Kamionkowski and March-Russell (1992) Marc Kamionkowski and John March-Russell, “Planck scale physics and the Peccei-Quinn mechanism,” Phys. Lett. B 282, 137–141 (1992), arXiv:hep-th/9202003 .
  * Kawasaki and Takahashi (2005) Masahiro Kawasaki and Fuminobu Takahashi, “Late-time entropy production due to the decay of domain walls,” Phys. Lett. B 618, 1–6 (2005), arXiv:hep-ph/0410158 .
  * Abel  _et al._ (2020) C. Abel _et al._ , “Measurement of the Permanent Electric Dipole Moment of the Neutron,” Phys. Rev. Lett. 124, 081803 (2020), arXiv:2001.11966 [hep-ex] .
  * Kawasaki  _et al._ (2015) Masahiro Kawasaki, Ken’ichi Saikawa, and Toyokazu Sekiguchi, “Axion dark matter from topological defects,” Phys. Rev. D 91, 065014 (2015), arXiv:1412.0789 [hep-ph] .