# Spatial distribution of two symmetric four-wave mixing signals induced by
Gaussian beams

M. R. L. da Motta  A. A. C. de Almeida  S. S. Vianna<EMAIL_ADDRESS>Departamento de Física, Universidade Federal de Pernambuco, 50670-901, Recife,
Pernambuco, Brazil

###### Abstract

We present a theoretical analysis of the spatial shape of two symmetric
signals of degenerate four-wave mixing induced by Gaussian beams in a thin
sample of two-level atoms. Our calculations take into account the full spatial
and spectral dependencies of the relevant nonlinear susceptibilities that
govern the two processes. This reveals two interesting effects. The first one
is that the total power of incident beams affects the transverse profile of
the four-wave mixing signals at the medium exit and their free propagation.
The second one is the influence of the spectral characteristics of the medium
on the longitudinal profile of both generated signals upon free propagation.
We argue that the first effect can be seen as the saturation of the medium in
regions of higher intensity, while the second can be understood as the result
of a nonlinear contribution to the refractive index inside the atomic medium.
These effects can be symmetric between the two signals, with asymmetries
induced by different detunings from resonance of the incident fields.

Usage

Structure

††preprint: APS/123-QED

##  I Introduction

The fundamental solution to the paraxial wave equation (PWE), the Gaussian
beam, encompasses important properties of the actual output of laser sources,
such as their characteristic intensity profile and phase distribution, both of
which evolve as the beam propagates. One exact higher-order solution to the
PWE is the Hermite-Gaussian (HG) mode, which is given in cartesian coordinates
and is characterized by lobes of light disposed in a rectangular grid. HG
modes are important in the study of laser cavities, as they describe the
spatial distribution of the output beams from these light sources [1]. The
Laguerre-Gaussian (LG) mode is another exact solution to the PWE. It is given
in cylindrical coordinates and carries a ring-shaped intensity distribution.
Most importantly, LG modes carry well-defined orbital angular momentum (OAM)
in the propagation direction per photon, as first demonstrated by Allen et al.
[2]. HG and LG modes can be represented in terms of one another and this
correspondence is well-known [2, 3]. Also, they are limiting cases of the so-
called Ince-Gaussian paraxial mode [4], that carries elliptical symmetry.
There are many other higher-order solutions to the PWE that present vastly
diverse characteristics and interesting properties [5].

The above discussion considers the radiation field of light, in its many
accessible spatial distributions or modes, propagating in free space. Here, we
are interested in the nonlinear light-matter interactions taking place in
atomic media and in understanding how these interactions affect the spatial
properties of the light field. We will be concerned with a particular
configuration of four-wave mixing (FWM), a third-order nonlinear optical
process that can take place in a variety of systems, such as atomic vapors,
cold atom samples, optical fibers, and others. FWM has allowed the
investigation of several optical phenomena, for instance, AC Stark shift,
phase conjugation and electromagnetically induced transparency [6, 7, 8, 9,
10]. It has also been widely employed to study the interaction of structured
optical modes with matter.

One setting that is commonly used to explore the spatial degrees of freedom of
light is that of FWM induced by amplified spontaneous emission in a hot atomic
vapor, with a 3-level cascade system [11, 12, 13, 14, 15]. In cold atomic
samples, FWM was employed to transfer OAM from incident to generated beams in
nondegenerate [16] and degenerate [17] atomic systems, to transfer more
complicated phase structures (obtained by superimposing LG modes of different
orders) [18], and to store the information carried by the spatial structure of
light in the ensemble of atoms, and later retrieve it [19, 20].

When the primary objective is the study of the spatial shape of the FWM beam,
the usual approach is based on the overlap integral of four paraxial modes
[12, 21, 14, 15]. In these cases, the calculations are performed regarding the
interaction medium as a channel for the nonlinear process to take place, not
possessing degrees of freedom that can affect the output mode superposition.
In other words, the spatial distribution of the FWM field is fully determined
solely by the distributions of the participating beams. With this approach,
the theoretical predictions are remarkably accurate [12, 15, 21]. The role of
the spatially dependent nonlinear coherence in the signal generation process
was discussed in Ref. [22]. In Ref. [23], the full spatial dependence of
medium quantities is taken into account in calculations, and effects of
detunings from resonance and phase mismatch on the phase distribution of the
FWM beam are evidenced. No connection is established, however, between the
medium quantities and the spatial properties of the beam outside the
interaction medium. This is the point our work seeks to highlight.

This work is a theoretical study of the spatial shape of two symmetric signals
of degenerate four-wave mixing generated in a sample of cold two-level atoms.
Our focus is on the influence of the spatially dependent nonlinear
susceptibility, induced as a result of the nonlinear process, on the overall
shape of the two FWM signals. We investigate the spatial properties of the
generated beams in various configurations, including different detunings from
the atomic resonances and several intensities. In particular, we show that the
total power of incident beams affects not only the transverse profile of both
FWM fields at the medium exit, but also the evolution of this profile under
free propagation. We show that the so-called root mean square (rms) parameters
of the generated beam, which serve as an effective measure of the longitudinal
profile, are intuitively affected by the frequency degrees of freedom of the
atomic medium, in a way that resembles the Kerr effect. These effects indicate
that even in a thin medium, the power and frequency of pump beams may lead to
significant modifications on the transverse and longitudinal profiles of a
nonlinear signal.

Our results focus on the thin-medium regime, characterized by a medium length
much smaller than the Rayleigh ranges of the participating fields. This regime
is in agreement with the experimental conditions in a cold atom sample
obtained with a magneto-optical trap (MOT). On the other hand, in the thick-
medium regime, the requirement of the Gouy phase-matching condition [15] can
lead to different results. This regime is more easily achievable in a hot
atomic vapor cell, for example.

This paper is divided as follows. In Sec. II, we present details on the
calculations performed, more specifically the semi-classical modeling of the
nonlinear interaction, the solution to the wave equation for the FWM beam, and
the properties of this solution. In Sec. III we present and discuss the main
results of this work. Finally, we establish our concluding remarks in Sec. IV.

##  II Four wave mixing in a two-level system

The theoretical model used to describe the generated FWM beam can be divided
in two main parts. The first one is the semi-classical description of the
interaction between the atomic medium and the radiation field of the laser
beams via optical Bloch’s equations. In particular, we are interested in the
case where both incident beams are strong and possess the same power. The
second part consists on solving the nonhomogeneous wave equation for the FWM
field $\mathcal{E}_{s}$ with the source term given by the nonlinear
polarization vector $\mathbf{P}_{\mathrm{NL}}$, related to the nonlinear
coherence obtained in the first part.

###  II.1 Nonlinear light-atom interaction

Figure 1: (a) Spatial orientation of incident and FWM fields. (b) Depiction of
the parametric processes that generate signals $S_{1}$ and $S_{2}$ in a two-
level atom.

We consider that both nonlinear signals are induced by two almost co-
propagating Gaussian laser beams, with wave-vectors $\mathbf{k}_{a}$ and
$\mathbf{k}_{b}$, and detected in the $2\mathbf{k}_{a}-\mathbf{k}_{b}$ and
$2\mathbf{k}_{b}-\mathbf{k}_{a}$ directions, as shown in Fig. 1(a). We employ
the density operator formalism to calculate the atomic medium response related
to the nonlinear processes represented in Fig. 1(b). We follow closely the
development of Ref. [7]. The total Hamiltonian is
$\hat{H}=\hat{H}_{o}+\hat{H}_{\mathrm{int}}$, where $\hat{H}_{o}$ is the free-
atom Hamiltonian, such that for the ground $\ket{1}$ and excited $\Ket{2}$
states, we have $\hat{H}_{o}\Ket{j}=\mathbb{E}_{j}\Ket{j}$, where
$\mathbb{E}_{j}$ is the energy eigenvalue of $\ket{j}$, $j\in\\{1,2\\}$. The
interaction Hamiltonian is
$\hat{H}_{\mathrm{int}}=-\hat{\bm{\mu}}\cdot\mathbf{E}\left(\mathbf{r},t\right)$,
where $\hat{\bm{\mu}}=e\hat{\mathbf{r}}$ is the electric dipole operator and
the total electric field is given by
$\mathbf{E}(\mathbf{r},t)=\sum_{\upsilon}\mathbf{E}_{\upsilon}(\mathbf{r},t)$,
$\upsilon\in\\{a,b\\}$,

$\begin{split}\mathbf{E}_{\upsilon}\left(\mathbf{r},t\right)&=\frac{1}{2}\bm{\epsilon}_{\upsilon}\mathcal{E}_{\upsilon}\left(\mathbf{r}\right)e^{-i\left(\mathbf{k}_{\upsilon}\cdot\mathbf{r}-\omega_{\upsilon}t\right)}+c.c.,\\\
&=\bm{\epsilon}_{\upsilon}E_{\upsilon}(\mathbf{r},t)+c.c.,\end{split}$ (1)

where $\mathbf{k}_{\upsilon}$ is the wave vector, $\bm{\epsilon}_{\upsilon}$
is the polarization direction, the amplitudes
$\mathcal{E}_{\upsilon}(\mathbf{r})$ carry the transverse dependence of the
fields and $c.c.$ means the complex conjugate. We consider a quasi co-
propagating configuration, where the angle $\vartheta$ between
$\mathbf{k}_{a}$ and $\mathbf{k}_{b}$ (Fig. 1(a)) is very small, making
$\mathbf{k}_{\upsilon}\cdot\mathbf{r}\simeq k_{\upsilon}z$, where
$k_{\upsilon}=|\mathbf{k}_{\upsilon}|=\omega_{\upsilon}n_{\upsilon}/c$ and
$n_{\upsilon}$ is the index of refraction at frequency $\omega_{\upsilon}$.

The polarization directions $\bm{\epsilon}_{a}$ and $\bm{\epsilon}_{b}$ of
input beams $\mathbf{E}_{a}$ and $\mathbf{E}_{b}$, respectively, determine the
number of atomic states involved in the nonlinear process. We are interested
in the case where $\bm{\epsilon}_{a}$ and $\bm{\epsilon}_{b}$ are parallel
circular polarizations. Thus, the light-atom interaction can be described in
terms of a two-level system. We consider the
$\Ket{5S_{1/2},F=2}\rightarrow\Ket{5P_{3/2},F=3}$ hyperfine transition of
${}^{87}\mathrm{Rb}$. The nonlinear interaction leads to the generation of
signal $S_{1}$, due to the absorption of two photons from beam $E_{a}$ and the
stimulated emission of one photon from beam $E_{b}$; and signal $S_{2}$, due
to the absorption of two photons from $E_{b}$ and the stimulated emission of
one photon from $E_{a}$. Figure 1(b) illustrates these processes
schematically.

We define the Rabi frequency

$\begin{split}\Omega_{\upsilon}(\mathbf{r})=\dfrac{\mu_{jk}\mathcal{E}_{\upsilon}\left(\mathbf{r}\right)}{2\hbar},\end{split}$
(2)

where $\mu_{jk}=\bra{j}(\hat{\bm{\mu}}\cdot\bm{\epsilon}_{\upsilon})\ket{k}$,
with $j\neq k$, is the electric dipole matrix element (due to parity,
$\mu_{jj}=0$). The matrix elements of $\hat{H}_{\mathrm{int}}$ can be writtten

$H_{\mathrm{int},jk}=-\hbar\sum_{\upsilon}\Omega_{\upsilon}e^{-i(k_{\upsilon}z-\omega_{\upsilon}t)}+c.c.$
(3)

The density operator $\hat{\rho}=\sum_{jk}\rho_{jk}\ket{j}\bra{k}$ describes
the state of the atomic ensemble and satisfies $\sum_{j}\rho_{jj}=1$, and
$\rho_{jk}=\rho_{kj}^{*}$, where $({}^{*})$ means complex conjugation. Its
time evolution is given by Liouville’s equation with a relaxation term
$\hat{\mathcal{L}}$ associated with spontaneous decay from $\ket{2}$ to
$\ket{1}$

$\frac{d\hat{\rho}}{dt}=\dfrac{i}{\hbar}[\hat{\rho},\hat{H}]+\hat{\mathcal{L}}.$
(4)

In our model, populations and coherences decay at rates $\Gamma$ and
$\Gamma/2$, respectively, where $\Gamma/2\pi$ is the natural decay rate of the
excited state. For the $D_{2}$ line of 87Rb, $\Gamma/2\pi\approx 6$ MHz [24].
With this, we obtain the optical Bloch’s equations (OBEs) for the two-level
system

$\begin{split}\dot{(\Delta\rho)}=&-\dfrac{2i}{\hbar}[\rho_{12}H_{\mathrm{int},21}-c.c.]-\Gamma\left[\Delta\rho-(\Delta\rho)^{0}\right],\\\
\dot{\rho}_{12}=&-\dfrac{i}{\hbar}\left[H_{\mathrm{int},12}\Delta\rho-\rho_{12}\left(\mathbb{E}_{2}-\mathbb{E}_{1}\right)\right]-\dfrac{\Gamma}{2}\rho_{12},\end{split}$
(5)

where $\Delta\rho=(\rho_{22}-\rho_{11})$ is the population difference and
$(\Delta\rho)^{0}$ is the population difference far from the region of
interaction with fields $E_{a}$ and $E_{b}$. We assume that the coherence
$\rho_{12}$ oscillates with frequencies $\omega_{a}$, $\omega_{b}$ and
$2\omega_{a}-\omega_{b}$ [7, 6]

$\begin{split}\rho_{12}&=\sigma_{12}^{a}e^{i\omega_{a}t}+\sigma_{12}^{b}e^{i\omega_{b}t}+\sigma_{12}^{2a-b}e^{i(2\omega_{a}-\omega_{b})t},\end{split}$
(6)

where $\sigma_{ij}$ are the slowly varying coherences. The
$2\omega_{a}-\omega_{b}$ component is responsible for the FWM process that
generates the signal $S_{1}$. The population difference $\Delta\rho$ has a
stationary component and one oscillating at $|\omega_{a}-\omega_{b}|$ [6],

$\Delta\rho=(\Delta\rho)^{\mathrm{dc}}+\left[(\Delta\rho)^{a-b}e^{i(\omega_{a}-\omega_{b})t}+c.c.\right].$
(7)

We now substitute equations (3), (6) and (7) into Eqs. (5), perform the
rotating wave approximation and collect terms that oscillate with the same
frequency. Then, in the steady state regime we arrive at the set of algebraic
equations shown in Appendix A for the slowly varying coherences and population
differences. Defining
$\widetilde{\Omega}_{\upsilon}\equiv\Omega_{\upsilon}e^{-ik_{\upsilon}z}$,
$\delta_{\upsilon}\equiv\omega_{\upsilon}-\omega_{o}$ as the detuning from
resonance of field $\Omega_{\upsilon}$ and
$\omega_{o}\equiv(\mathbb{E}_{2}-\mathbb{E}_{1})/\hbar$ as the resonance
frequency, we obtain

$\sigma_{12}^{2a-b}=\dfrac{-2i\widetilde{\Omega}_{a}^{2}\widetilde{\Omega}_{b}^{*}(\Delta\rho)^{\mathrm{dc}}(1/\Delta_{a}+1/\Delta_{b})}{(\Delta_{a}+\Delta_{b})\left(2i\delta_{a}-i\delta_{b}+\Gamma/2\right)+2|\Omega_{a}|^{2}},$
(8)

where $\Delta_{a}=i\delta_{a}+\Gamma/2$, $\Delta_{b}=-i\delta_{b}+\Gamma/2$
and

$(\Delta\rho)^{\mathrm{dc}}=\dfrac{(\Delta\rho)^{0}}{1+\dfrac{2|\Omega_{a}|^{2}}{\delta_{a}^{2}+\Gamma^{2}/4}+\dfrac{2|\Omega_{b}|^{2}}{\delta_{b}^{2}+\Gamma^{2}/4}}.$
(9)

The coherences $\sigma_{12}^{a}$ and $\sigma_{12}^{b}$ can be found as

$\displaystyle\sigma_{12}^{a}$
$\displaystyle=\dfrac{-i\widetilde{\Omega}_{a}\rho_{11}^{0}}{\Delta_{a}+2|\Omega_{a}|^{2}/\Gamma+2|\Omega_{b}|^{2}/\Delta_{ab}},$
(10) $\displaystyle\sigma_{12}^{b}$
$\displaystyle=\dfrac{-i\widetilde{\Omega}_{b}\rho_{11}^{0}}{\Delta^{*}_{b}+2|\Omega_{b}|^{2}/\Gamma+2|\Omega_{a}|^{2}/\Delta_{ab}^{*}},$
(11)

where $\Delta_{ab}=i\delta_{a}-i\delta_{b}+\Gamma/2$ and $\rho_{11}^{0}=1$ is
the population of the ground state far from the region of interaction.

We rewrite the coherences given by Eqs. (8), (10) and (11) in the forms

$\displaystyle\sigma_{12}^{2a-b}$
$\displaystyle=X^{2a-b}\widetilde{\Omega}_{a}^{2}\widetilde{\Omega}_{b}^{*},$
(12) $\displaystyle\sigma_{12}^{a}$
$\displaystyle=X^{a}\widetilde{\Omega}_{a},$ (13)
$\displaystyle\sigma_{12}^{b}$ $\displaystyle=X^{b}\widetilde{\Omega}_{b},$
(14)

where $X^{2a-b}$, $X^{a}$ and $X^{b}$ are the couplings associated with the
processes in directions $(2\mathbf{k}_{a}-\mathbf{k}_{b})$, $\mathbf{k}_{a}$
and $\mathbf{k}_{b}$, respectively. These coupling factors carry all the
spectral response of the medium. They also depend on the field amplitudes
$|\Omega_{a,b}|^{2}$, and thus the corresponding susceptibilities contain
information regarding higher order processes in the same direction. We define
the effective susceptibilities

$\displaystyle\chi^{2a-b}$
$\displaystyle=\dfrac{\mathcal{N}|\mu_{12}|^{4}}{\varepsilon_{o}\hbar^{3}}X^{2a-b},$
(15) $\displaystyle\chi^{ab}$
$\displaystyle=\dfrac{\mathcal{N}|\mu_{12}|^{2}}{\varepsilon_{o}\hbar}(X^{a}+X^{b}),$
(16)

where $\mathcal{N}$ is the atomic density. Since the atoms are considered
stationary, we do not need to include the effect of Doppler broadening.

###  II.2 Wave equation

We are interested in the field distributions of the generated signals, and
thus seek a solution to the wave equation for field $\mathbf{E}_{s}$. We focus
only on $S_{1}$, since the equations for $S_{2}$ are obtained and solved in
the same way. The wave equation for the FWM electric field is [25]

$\nabla^{2}\mathbf{E}_{s}-\dfrac{1}{c^{2}}\dfrac{\partial^{2}\mathbf{E}_{s}}{\partial
t^{2}}=\mu_{o}\dfrac{\partial^{2}\mathbf{P}}{\partial t^{2}},$ (17)

where $\mathbf{E}_{s}$ is written as in Eq. (1). The macroscopic polarization
$\mathbf{P}$ can be divided in two components: $\mathbf{P}_{2a-b}$, which
describes the generation process of the FWM field $E_{s}$; and
$\mathbf{P}_{ab}$, associated with the propagation of the generated field
inside the medium affected by fields $E_{a}$ and $E_{b}$. We can write the
total polarization as $\mathbf{P}=\mathbf{P}_{ab}+\mathbf{P}_{2a-b}$, and its
projection onto the oscillation direction of the generated field is

$(\mathbf{P}\cdot\bm{\epsilon}^{*}_{s})=\varepsilon_{o}\chi^{ab}E_{s}+\varepsilon_{o}\chi^{2a-b}E_{a}^{2}E_{b}^{*}.$
(18)

In the calculations that follow, and throughout the rest of this work, we
neglect the first term on the right-hand side of Eq. (18). This is done under
the thin-medium regime assumption, such that the generated field is not
influenced by the medium response associated with $\chi^{ab}$. Under the
rotating wave and paraxial approximations, we obtain

$\left(\frac{i}{2k_{s}}\nabla^{2}_{\perp}+\frac{\partial}{\partial
z}\right)\mathcal{E}_{s}=\kappa\Omega_{a}^{2}\Omega_{b}^{*}e^{-i\Delta kz},$
(19)

where $\nabla^{2}_{\perp}$ is the transverse Laplacian, $\Delta
k=|2\mathbf{k}_{a}-\mathbf{k}_{b}-\mathbf{k}_{s}|\simeq
2k_{a}-(k_{b}+k_{s})\cos\vartheta$ is the phase mismatch and $\kappa$ is the
nonlinear coupling, given by

$\displaystyle\kappa(r,z;\delta)=-i\dfrac{\omega_{s}\hbar^{3}}{2c\mu_{12}|\mu_{12}|^{2}}\chi^{2a-b}(r,z;\delta).$
(20)

In the last equation, $\delta$ represents $\delta_{a}$ and $\delta_{b}$.

We highlight that there is an implicit position dependence in $\kappa$. It
exists only inside the interaction region, between $z=-L/2$ and $z=L/2$, and
is zero everywhere else in the $z$-axis. Thus, Eq. (19) describes the
nonlinear signal generation process inside the sample. There is no nonlinear
signal at positions $z\leq-L/2$, and so the boundary condition is
$\mathcal{E}_{s}(\mathbf{r}_{\perp},-L/2)=0$. For $z>L/2$, where the generated
beam propagates in free space, $\mathcal{E}_{s}$ must satisfy the homogeneous
paraxial wave equation, $(i\nabla^{2}_{\perp}/2k_{s}+\partial/\partial
z)\mathcal{E}_{s}=0$, with the boundary condition given by the solution of Eq.
(19) at $z=L/2$, $\mathcal{E}_{s}(\mathbf{r}_{\perp},L/2)$.

For the solution of Eq. (19), we first note that the incident beams are
strong, and we expect them to undergo little extinction and transverse
structure variation along the interaction region. Thus, we treat Eq. (19)
uncoupled from the wave equations for $\mathcal{E}_{a}$ and $\mathcal{E}_{b}$.
This simplifies the mathematical work of our problem and is shown to bring
results with good agreement with experimental measurements [21, 15, 12].

###  II.3 Solution to the FWM field wave equation

It is well known that both the HG and LG paraxial modes form complete
orthonormal sets of functions on the transverse plane. This property allows to
write any scalar optical field as a superposition of the form
$U(\mathbf{r}_{\perp},z)=\sum_{m,n}c_{mn}(z)u_{mn}(\mathbf{r}_{\perp},z)$,
where indices $(m,n)$ characterize the modes of the chosen basis and, in
analogy with a quantum-mechanical system, the expansion coefficients
$c_{mn}(z)$ can be seen as probability amplitudes of finding the system
$\Ket{U}$ in the state $\ket{m,n}$.

We consider only incident Gaussian beams, and therefore, the LG basis is most
convenient due to its cylindrical symmetry. The LG mode is denoted as

$\begin{split}u_{\ell p}(r,\phi,z)&=\frac{C_{\ell
p}}{w(z)}\left(\frac{\sqrt{2}r}{w(z)}\right)^{|\ell|}L_{p}^{|\ell|}\left[\frac{2r^{2}}{w^{2}(z)}\right]e^{i\ell\phi}\\\
&\times
e^{-\frac{r^{2}}{w^{2}(z)}}\exp\left[-i\frac{kr^{2}}{2R(z)}+i\Psi_{\mathrm{G}}(z)\right],\end{split}$
(21)

where $C_{\ell p}=\sqrt{2p!/\pi(p+|\ell|)!}$ is the normalization constant,
$L_{p}^{|\ell|}(\cdot)$ is the associated Laguerre polynomial,
$w(z)=w_{o}\sqrt{1+(z/z_{R})^{2}}$ is the beam waist,
$R(z)=z\left[1+(z_{R}/z)^{2}\right]$ is the curvature radius,
$\Psi_{\mathrm{G}}(z)=(N_{\ell p}+1)\tan^{-1}(z/z_{R})$ is the Gouy phase
shift, with the total mode order defined as $N_{\ell p}=2p+|\ell|$,
$z_{R}=kw_{o}^{2}/2$ is the Rayleigh range and $w_{o}$ is the minimum beam
waist.

A light beam described by an LG mode carries well-defined OAM in the
$z$-direction, which is related to the azimuthal phase factor $e^{i\ell\phi}$,
where integer $\ell\in(-\infty,\infty)$, called the topological charge,
defines the OAM per photon in the beam [2]. The other index characterizing the
mode, $p\in[0,\infty)$, is called the radial index. It is related to the
number of dark rings in the intensity profile of $u_{\ell p}$, but does not
have a straightforward connection with a physical quantity as is the case for
$\ell$. In recent years, however, the radial index has been the subject of
theoretical works [26, 27] that have enlightened its significance.

We write the generated field amplitude $\mathcal{E}_{s}$ as the superposition

$\mathcal{E}_{s}(\mathbf{r})=\sum_{\ell=-\infty}^{\infty}\sum_{p=0}^{\infty}\mathcal{A}_{\ell
p}(z)u_{\ell p}(\mathbf{r}).$ (22)

The problem becomes that of finding the set of relevant coefficients
$\\{\mathcal{A}_{\ell p}\\}$. Substituting Eq. (22) into Eq. (19) and
employing the orthogonality relation of $u_{\ell p}$, we obtain an equation
for $\mathcal{A}_{\ell p}(z)$.

$\dfrac{\partial\mathcal{A}_{\ell p}(z)}{\partial
z}=\Lambda_{p}^{\ell}(z)e^{-i\Delta kz},$ (23)

where

$\Lambda_{p}^{\ell}(z)=\int_{0}^{2\pi}\int_{0}^{\infty}\kappa(\mathbf{r})\Omega_{a}^{2}\Omega_{b}^{*}u_{\ell
p}^{*}rdrd\phi,$ (24)

is the projection of the spatially dependent nonlinear source term onto the LG
function space, called the transverse overlap integral.

As we consider only Gaussian incident beams, that carry no OAM and thus
possess $\ell=0$, it will be useful to define the single index mode
$u_{p}(r,z)\equiv u_{0p}(r,z)$, which is azimuthally symmetric. Since
$\kappa(\mathbf{r})$ only contains the squared modulus of fields
$\Omega_{a,b}$, it has no $\phi$ dependence, even for $\ell_{a,b}\neq 0$. The
azimuthal integral in Eq. (24) yields
$\int^{2\pi}_{0}e^{i\ell\phi}d\phi=2\pi\delta_{\ell,0}$, and
$\Lambda^{\ell}_{p}(z)$ becomes
$\Lambda_{p}(z)\equiv\Lambda^{0}_{p}(z)=2\pi\int^{\infty}_{0}\kappa\Omega^{2}_{a}\Omega^{*}_{b}u^{*}_{p}rdr$.
We then see that, as the incident fields have $\ell=0$, the FWM field does not
carry OAM content, as anticipated, and can be written as a superposition of
radial modes $\mathcal{E}_{s}=\sum_{p}\mathcal{A}_{p}(z)u_{p}(r,z)$, where
$\mathcal{A}_{p}\equiv\mathcal{A}_{\ell=0,p}$. The presence of modes with
different orders $N_{p}=N_{0p}=2p$ can lead to the interference of the various
Gouy phase factors $\exp{[i(1+N_{p})\tan^{-1}(z/z_{R})]}$ that affects the
radial structure of $|\mathcal{E}_{s}|^{2}$ upon propagation [28, 29]. When
incident fields carry single topological charges or an arbitrary superposition
of topological charges, the azimuthal integral imposes the conservation of OAM
[12, 14, 15]. This is the subject of future work currently underway.

We then solve Eq. (23) to find the expansion coefficients as

$\mathcal{A}_{p}(z)=\int_{-L/2}^{z_{<}}\Lambda_{p}(z^{\prime})e^{-i\Delta
kz^{\prime}}dz^{\prime},$ (25)

where $z_{<}=\min[z,L/2]$. This solution is suitable for both regions of
space: $-L/2<z\leq L/2$ and $z>L/2$. We are interested in the FWM beam outside
the medium, where it can be detected, and thus seek to evaluate
$\mathcal{A}_{p}(L/2)$. It is important to note that, although not explicit,
$\mathcal{A}_{p}(z)$ is also a function of the detunings $\delta_{a}$ and
$\delta_{b}$.

The characteristic length of $\Lambda_{p}(z)$ is given by the Rayleigh range
of the beams that participate in the FWM process, $z_{R}$. For a thin-medium,
characterized by $L\ll z_{R}$, we can neglect the variation of
$\Lambda_{p}(z)$ inside the interaction region and take its value at $z=0$.
This allows to remove the transverse overlap from the $z$ integral in Eq. (25)
and write the approximate form

$\begin{split}\mathcal{A}_{p}(L/2)&\simeq\Lambda_{p}(0)\int_{-L/2}^{L/2}e^{-i\Delta
kz^{\prime}}dz^{\prime},\\\ &=\Lambda_{p}(0)T(L).\end{split}$ (26)

In the above expression, $T(L)=2\sin(\Delta kL/2)/\Delta k$ can be regarded as
an efficiency measure of the signal generation process inside the medium. It
takes into account the phase-mismatch and is a common factor to all
$\mathcal{A}_{p}$. Thus, in the thin-medium regime, all the information of the
nonlinear wave mixing process is contained in the transverse overlap integral
evaluated at $z=0$, $\Lambda_{p}(0)$. It is therefore the main quantity in our
calculations, and completely determines the distribution of modes of the
generated fields.

The orthogonality of LG modes allows to write the total power of field
$\mathcal{E}_{s}$ as $P=\sum_{p}|\mathcal{A}_{p}|^{2}$, which gives the
normalization factor of $\mathcal{E}_{s}$ at fixed $\delta$. We define the
mode purity or mode weight as $\eta_{p}\equiv|\mathcal{A}_{p}|^{2}/P$, a
measure of the relative contribution of the mode $u_{p}$ to $\mathcal{E}_{s}$.
The phase angle of $\mathcal{A}_{p}$ is $\Phi_{p}$, such that we can write the
expansion coefficient in the normalized form
$\mathcal{A}_{p}=\sqrt{P\eta_{p}}\,e^{i\Phi_{p}}$. Note that the relative
phases between the various expansion coefficients can also be responsible for
changes on the output beam superposition.

The nonlinear coupling $\kappa(r,z;\delta)\propto\chi^{2a-b}(r,z;\delta)$ has
a complicated dependence on the input fields, and, rigorously, on the position
$\mathbf{r}$. However, even though this coupling can promote sensible
modifications to the FWM beam in the near-field, it does not contribute
substantially to the transverse shape of the _far-field_ FWM beam, which is
mainly determined by the mode components of fields $\Omega_{a}$ and
$\Omega_{b}$. This suggests that the FWM beam profile is dictated by the
overlap of input beams. Indeed, this is usually assumed in the description of
nonlinear processes involving beams with OAM or arbitrary transverse
structures. In references [30, 28, 31, 32], that focus on second-order
nonlinearities (parametric oscillation and second-harmonic generation), and in
references [12, 14, 15], that treat FWM, the quantities analogous to that of
Eq. (24) are overlap integrals of three ($\chi^{(2)}$) and four ($\chi^{(3)}$)
LG modes. In all of these cases, the relevant nonlinear susceptibility,
$\chi^{(2)}$ or $\chi^{(3)}$, is a uniform quantity that factors out of the
overlap integral and does not affect the value of the mode expansion
coefficients. With the present study, we seek to understand the influence of
the atomic medium on the FWM process. More specifically, the role of the full
spatial and spectral dependencies of $\chi^{2a-b}(\mathbf{r};\delta)$ on the
spatial features of signals $S_{1}$ and $S_{2}$ in the particular FWM
configuration considered. Our results, in Sec. III, show that by varying the
power and the frequency of input fields, one can produce changes on the
distributions of mode weights $\eta_{p}$ and phases $\Phi_{p}$, which in turn
generate modifications on the FWM beams upon free space propagation. This is
because in our case, the effective susceptibilty can not be factored out from
the overlap integral, and is taken fully into account in the calculations.

###  II.4 Longitudinal profile of the FWM field

The spot size of the FWM beam is well described by the rms radius, defined as
$r_{\mathrm{rms}}(z)\equiv[\frac{1}{P}\iint
r^{2}|\mathcal{E}_{s}|^{2}rdrd\phi]^{\frac{1}{2}}$ [33]. Substituting Eq.
(22), at $z>L/2$ we can obtain the form [34]

$r_{\mathrm{rms}}(z)=\sqrt{r^{2}_{m}+\theta^{2}_{\mathrm{rms}}(z-z_{m})^{2}},$
(27)

where $r_{m}$ is the minimum spot size, $z_{m}$ is the position where it
occurs and $\theta_{\mathrm{rms}}$ is the divergence angle. These three
parameters are calculated from the set of coefficients $\\{\mathcal{A}_{p}\\}$
(see Appendix B) and completely determine the longitudinal profile of the
beam. We can also introduce the dimensionless quantity
$\mathcal{M}^{2}=k\theta_{\mathrm{rms}}r_{m}$, called the beam quality factor,
which is proven to satisfy the bound $\mathcal{M}^{2}\geq
1+\langle|\ell|\rangle$ [34], where $\langle{|\ell|}\rangle$ is the mean value
of the topological charge magnitude in the field expansion. In our case,
$\langle|\ell|\rangle=0$ and the quality factor must therefore satisfy
$\mathcal{M}^{2}\geq 1$.

We highlight that $\theta_{\mathrm{rms}}$, $r_{m}$, $z_{m}$, and
$\mathcal{M}^{2}$, the so-called rms quantities, depend on the detunings
$\delta_{a,b}$. In order to highlight that this dependence arises in the
theory due to the spatial distribution of $\chi^{2a-b}$, we consider the
spatially uniform nonlinear coupling
$\overline{\kappa}(\delta)\propto\overline{\chi}(\delta)$, an averaged measure
of the nonlinear response over the interaction volume. Then,
$\overline{\kappa}(\delta)$ factors out of the integral given by Eq. (24),
which becomes the integral of a product of four LG modes [12, 14, 15]. It can
be seen that in this case all of the dependence on the frequency degrees of
freedom factor out from the coefficients $\mathcal{A}_{p}$ and the rms
quantities do not vary with $\delta$.

##  III Results and Discussion

In this section, we present and discuss our main results. Together, they
comprise a theoretical investigation of various aspects of the FWM signal
generation (induced by Gaussian incident beams) and the free space propagation
after leaving the nonlinear medium.

Our calculations consist of evaluating the expansion coefficients
$\mathcal{A}_{p}$ for the two symmetric signals $S_{1}$ and $S_{2}$. With the
set $\\{\mathcal{A}_{p}\\}$, all of the quantities we are interested in can be
obtained, such as the intensity distribution $|\Omega_{s}|^{2}$ and its
propagation outside the medium, the mode components $\eta_{p}$, and others. We
do this for several configurations. Here, each configuration is defined by the
amplitudes $\Omega_{a}^{0}$ and $\Omega_{b}^{0}$, related to the total power
of the incident beams
($\Omega_{a,b}(\mathbf{r})=\Omega^{0}_{a,b}u_{0}(\mathbf{r})$); and the
detunings from resonance $\delta_{a}$ and $\delta_{b}$. In the calculations,
all beams possess the same wavelength $\lambda$ and minimum waist $w_{o}$,
giving the same Rayleigh range $z_{R}$.

We consider that the waist of all beams is $w_{o}=1$ mm near the interaction
region. For the wavelength $\lambda=780$ nm, used to excite the $D_{2}$ line
of 87Rb, the Rayleigh range is $z_{R}=\pi w_{o}^{2}/\lambda\approx 4$ m. A
cold atom cloud usually obtained with a MOT has a size $L$ of a few
millimeters. The condition $L/z_{R}\ll 1$ is satisfied, thus justifying our
focus on the thin medium regime.

As will be seen, for purely Gaussian incident beams, the transverse shapes of
the FWM beams are mainly Gaussian but can have significant contributions from
modes with $p>0$. We show how the transverse profiles of the generated FWM
beams and their free space propagation are affected by the total power
contained in the incident beams and by the detunings from resonance
$\delta_{a}$ and $\delta_{b}$.

###  III.1 FWM induced by Gaussian beams

For Gaussian input beams, the two signals $S_{1}$ and $S_{2}$ carry null
topological charges, $\ell_{1}=\ell_{2}=0$. Also, for equal amplitudes
$\Omega_{a}^{0}=\Omega_{b}^{0}$ and detunings $\delta_{a}=\delta_{b}=\delta$,
they possess completely symmetric transverse shapes and free space propagation
characteristics. This is because the nonlinear coherence, given by Eq. (8),
remains the same with an exchange of labels $a\leftrightarrow b$.

####  III.1.1 Effect of pump intensity

Figure 2: (a) Mode weights ${\eta}_{p}$, (b) propagation of the FWM beam
intensity profile outside the interaction medium from $z/z_{R}=L/2z_{R}\approx
0$ to $z/z_{R}=1$ (insets show the intensity profile at $z=L/2,z_{R}/2,z_{R}$)
and (c) normalized radial distribution of intensity at positions
$z=L/2,z_{R}/2$ and $z_{R}$, for incident Gaussian beams with Rabi frequencies
$\Omega_{a,b}^{0}=\Omega=0.5\Gamma,\Gamma,2\Gamma$ and $4\Gamma$ (from top to
bottom).

We now investigate the effect of the total power contained in the incident
beams on the overall shape of the generated beams. The weights ${\eta}_{p}$,
that quantify the contribution of each mode $u_{p}$ to the generated field
$\mathcal{E}_{s}$, is represented in Fig. 2(a), and the radial distribution of
intensity at different longitudinal positions is discriminated in Fig. 2(c).
Figure 2(b) shows the FWM beam propagation for $\delta_{a}=\delta_{b}=0$ and
different values of $\Omega_{a}^{0}=\Omega_{b}^{0}=\Omega$. We note that as
the amplitudes increase, the overall shape of the generated beam at the
nonlinear medium exit changes significantly. However, keeping
$\Omega_{a}^{0}=\Omega_{b}^{0}$, the symmetry between both generated beams is
preserved. Near the medium exit, $z/z_{R}=L/2z_{R}\approx 0$, the beam becomes
ring-shaped. We attribute this behavior to a spatial saturation effect, which
can be understood by inspecting equation (8). In all cases, the near-field
intensity profile is determined by the nonlinear coherence,
$|\sigma^{2a-b}_{12}|^{2}\propto|\chi^{2a-b}|^{2}I^{2}_{a}I_{b}$, and due to
the Gaussian distribution of the fields, for greater amplitudes, the
denominator in Eq. (8) is larger at the center, making
$|\sigma_{12}^{2a-b}|^{2}$ smaller in this region. Nonetheless, after
propagating distances of the order of $z_{R}$, the beams acquire a shape that
corresponds to the dominant modes, as we can see in the radial profile at
$z=z_{R}$ (green line) in Fig. 2(c). This is similar to what is verified in
Ref. [28].

For the lower intensities ($\Omega=0.5\Gamma,\Gamma$) the dominant mode is
$u_{0}$, and the far-field profile is mainly Gaussian. In this case, the
saturation effects are very small, making the influence of the nonlinear
coupling $\kappa(\mathbf{r})$ in Eq. (24) negligible. The values of $\eta_{p}$
approach those obtained by calculating the overlap integral of four LG modes
[12, 14, 15]. On the other hand, for the higher intensities
($\Omega=2\Gamma,4\Gamma$) there is a significant contribution from the mode
$u_{1}$, and a dark ring is present on the far-field profile. This transition
of the radial profile is due to the superposition of modes in the generated
beam that possess different mode numbers $N_{p}$, and thus acquire different
Gouy phases upon propagation [28, 29].

Sensible change of the longitudinal profile is achieved by increasing
$\Omega$. We highlight that the generated beam outside the medium clearly
indicates an intensity maximum at positions $z/z_{R}$ far from $0$, for all
incident beam amplitudes. In fact, for increasing $\Omega$, this position of
maximum intensity is shifted towards greater $z/z_{R}$, as seen from Fig.
2(b). This would make one expect the position of minimum $r_{\mathrm{rms}}(z)$
to be shifted as well. However, this is not the case, since the minimum of
$r_{\mathrm{rms}}(z)$ does not necessarily correspond to the position of
maximum local intensity. This is because $r_{\mathrm{rms}}(z)$ for an
arbitrary beam does not correspond to the radial position where the field
amplitude decreases by a factor of $1/e$ relative to the amplitude at the
center, as it does for a pure Gaussian beam.

####  III.1.2 Effect of detunings from resonance

Figure 3: (a) Behavior of $r_{\mathrm{rms}}(z)$ of the symmetric generated
beams with varying $\delta_{a}=\delta_{b}=\delta_{a,b}$ for $\Omega=\Gamma$
and $\Gamma/4$ (inset). (b) Longitudinal parameters for the symmetric
generated beams for the same values of $\delta_{a,b}$. Red and black curves
correspond to $\Omega=\Gamma$ and $\Omega=\Gamma/4$, respectively. Figure 4:
(a) Behavior of $r_{\mathrm{rms}}(z)$ of signals $S_{1}$ (top) and $S_{2}$
(bottom), and (b) longitudinal parameters of $S_{1}$ (red) and $S_{2}$ (black)
for $\delta_{a}=0$, and different values of $\delta_{b}/\Gamma$. Rabi
frequencies are $\Omega_{a,b}^{0}=\Gamma$.

To better understand the influence of the medium on the spatial
characteristics of the generated beams, we investigate the effect of the
detunings from resonance, $\delta_{a,b}$, on the free space propagation of the
generated beam. First, we consider a situation where both beams have equal
detunings, i.e., come from the same laser source. We see from Fig. 3 that
$\delta_{a,b}$ has an intuitive effect on the FWM beams. On resonance, the
beams are generated with the maximum radius and the minimum divergence. As we
move away from the resonance, the radius right outside the sample decreases
while the divergence angle increases. It is interesting to note that $z_{m}$
changes considerably for varying incident beam detunings. Above resonance,
$z_{m}$ is shifted to negative values, while below resonance, it is shifted to
positive values. These results suggest that it is possible to translate the
position where the minimum waist of the FWM beam occurs by controlling the
frequency of the incident beams. This translation comes with not much change
on the other beam parameters. Another parameter is the beam quality factor,
which is maximum at resonance and approaches unity as $\delta_{a,b}$ goes away
from resonance.

To evidence that the influence of the frequency degrees of freedom depends on
the incident beam amplitudes, we show in the inset of Fig. 3(a) the rms radii
of the symmetric generated beams when $\Omega_{a,b}^{0}=\Omega=\Gamma/4$. We
see that all five curves, at each $\delta_{a,b}$ value, are now closer to each
other. The longitudinal parameters for this low Rabi frequency are shown in
the black curves of Fig. 3(b), and we can see that $r_{m}$ varies much less
with the detunings and $z_{m}$ stays around $0$ for all $\delta_{a,b}$. This
is because, in the limit $\Omega^{0}_{a,b}\rightarrow 0$, the coupling
$\kappa$ becomes uniform and factors out of the integral in Eq. (24). As a
result, the frequency degrees of freedom would not affect the overall shape of
the generated beam, only the power conversion efficiency. Moreover, we can see
that the beam quality factor satisfies $\mathcal{M}^{2}\geq 1$ for all
detunings and approaches the value $\mathcal{M}^{2}=1$ for increasing
$|\delta_{a,b}|$.

We now turn to a situation where the two incident beams do not have the same
detuning but maintain the Rabi frequencies equal,
$\Omega_{a}^{0}=\Omega_{b}^{0}$. In this case, the symmetry under the exchange
$a\leftrightarrow b$ in Eq. (8) no longer holds, and the generated signals
$S_{1}$ and $S_{2}$ are shown to differ.

First, we set the frequency of $E_{a}$ on resonance, $\delta_{a}=0$, and make
$\delta_{b}$ vary around $\delta_{b}=0$. Figure 4(a) shows the radii of both
generated beams on free-space propagation in this situation. We see that by
varying only the detuning $\delta_{b}$, we also obtain changes in the focusing
region of both FWM beams. However, the positions of minimum radius, $z_{m}$,
of the two signals are translated to opposite directions (Figure 4(b)). A more
focused beam on direction $(2\mathbf{k}_{a}-\mathbf{k}_{b})$ is accompanied by
a more spread beam on direction $(2\mathbf{k}_{b}-\mathbf{k}_{a})$, and vice
versa. The values of the divergence angle $\theta_{\mathrm{rms}}$, minimum
radius $r_{m}$, and quality factor $\mathcal{M}^{2}$, also shown in Figure
4(b), are very similar for both signals. Further, their dependence on
$\delta_{b}$ is similar to that seen in the case with $\delta_{a}=\delta_{b}$.

In Figure 5(a) we show the longitudinal profile of $S_{1}$ and $S_{2}$ for
$\delta_{a}=0$ and $\delta_{b}=0.75\Gamma$, corresponding to the purple curves
of Figure 4(a). Referring to the discussion regarding the position of maximum
intensity, we see that in this case, the position of minimum
$r_{\mathrm{rms}}(z)$ seems to be closer to the position of maximum intensity
outside the nonlinear medium.

Figure 5: (a) Propagation of both generated beams outside the nonlinear medium
from $z/z_{R}=L/2z_{R}\approx 0$ to $z/z_{R}=1$, (b) mode components
${\eta}_{p}$ and (c) phases $\Phi_{p}$ of $S_{1}$ and $S_{2}$ for
$\delta_{a}=0$ and $\delta_{b}=0.75\Gamma$, corresponding to the purple curves
of Fig. 4(a).

The distributions of $\eta_{p}$, shown in Fig. 5(b) are similar for both
signals, with $u_{0}$ having the greatest contribution, as before, and
slightly different weights for the modes with $p\neq 0$. However, the
distribution of relative phases of the modes $u_{p}$ that are being
superimposed, shown in Fig. 5(c), is quite different between $S_{1}$ and
$S_{2}$. This is the dominant factor that leads to the differences in the
longitudinal profiles of Fig. 5(a).

Our calculations reveal an effect that resembles Kerr lensing [35], where the
total index of refraction in the medium can be written as $n=n_{0}+n_{2}I$,
where $n_{0}$ and $n_{2}$ are the linear and nonlinear refractive indices and
$I$ is the beam intensity. The total index of refraction seen by the strong
beam is thus modulated by its own spatial distribution. Evidently, the
generated beam is not strong enough to induce self-modulation. However, we can
explain qualitatively the observed focusing/defocusing effect by expressing
the total index of refraction seen by the generated field as

$n_{t}=n_{0}+\Delta n(I_{a},I_{b};\delta),$ (28)

where $\Delta n(I_{a},I_{b};\delta)\propto\mathrm{Re}\chi^{2a-b}(r;\delta)$ is
a nonlinear contribution modulated by the spatial distribution of both
incident beams.

Now we look at the first situation considered: incident fields with equal
detunings, $\delta_{a}=\delta_{b}$, and symmetric nonlinear signals. We show
in figures 6(a) and 6(b) plots of the real part of the nonlinear
susceptibility as a function of the radial coordinate for different values of
$\delta_{a}=\delta_{b}$. The Rabi frequencies are $\Omega_{a,b}^{0}=\Gamma$
(Fig. 6(a)) and $\Omega_{a,b}^{0}=\Gamma/4$ (Fig. 6(b)), corresponding to the
situations of Fig. 3(a) and its inset, respectively.

Figure 6: Radial distributions of the real part of the nonlinear
susceptibility for different cases. For (a)-(b) the detunings are equal
$\delta_{a}=\delta_{b}=\delta$, depicting the completely symmetric case,
$\chi^{2a-b}=\chi^{2b-a}=\chi$. Rabi frequencies are (a)
$\Omega^{0}_{a,b}=\Gamma$ and (b) $\Omega^{0}_{a,b}=\Gamma/4$. For (c)-(d)
$\delta_{a}$ is fixed on resonance, and we have the non-symmetric case. In (c)
and (d) we show respectively $\mathrm{Re}\chi^{2a-b}(r;\delta_{b})$ and
$\mathrm{Re}\chi^{2b-a}(r;\delta_{b})$ with $\Omega^{0}_{a,b}=\Gamma$.

We see that below resonance, $\delta_{a,b}<0$ (red and blue curves), the total
index of refraction is greater at the center, $r=0$, and decreases at greater
$r$ positions. Thus, the FWM beam is focused. On resonance, $\delta_{a,b}=0$,
the total index of refraction is unaffected by $\chi^{2a-b}=\chi^{2b-a}=\chi$
at all radial positions. The FWM beam is neither focused nor defocused. Above
resonance, $\delta_{a,b}>0$ (purple and green curves), $n_{t}$ is smaller at
the center, and increases as we move away from this position. As a result, the
FWM beam is defocused. For the lower intensity (Fig. 6(b)) the variation of
the refractive index is smaller. We then expect the (de)focusing effect to be
weaker. This agrees with the behavior of $z_{m}$ shown in Fig. 3 and its
inset.

The same analysis can be made regarding the second case, where $\delta_{a}$ is
fixed at resonance, and $\delta_{b}$ varies. As already discussed, the two FWM
signals in this case are non-symmetric. In figures 6(c) and 6(d) we present
plots of the real part of the nonlinear susceptibility associated with the
generation of FWM signals $S_{1}$ ($2\mathbf{k}_{a}-\mathbf{k}_{b}$) and
$S_{2}$ ($2\mathbf{k}_{b}-\mathbf{k}_{a}$) with the input Rabi frequency
$\Omega^{0}_{a,b}=\Gamma$. The nonlinear susceptibility related to $S_{2}$ has
a radial dependence similar to that shown in Fig. 6(a)-(b) for $\delta_{b}>0$
and $\delta_{b}<0$, while for $S_{1}$, the curves are switched. This indicates
that the (de)focusing effect is opposite between $S_{1}$ and $S_{2}$, in
agreement with the results for $z_{m}$ shown in Fig. 4(b). The opposite
behavior of the two generated beams is related to whether the beams that
participate with one and two photons in the nonlinear process have negative or
positive detunings.

When the incident fields carry topological charge, the radial profile of the
susceptibility is further affected by their non-Gaussian intensity profiles.
More interestingly, the phase distribution of the generated field is twisted
when the detunings are varied around resonance [23, 36]. The phase
discontinuities that arise from the azimuthal phase factors $e^{\pm
i\ell\phi}$, inherent to vortex beams, are the features that reveal this
twisting effect. We are conducting a separate study that deals with this
situation.

##  IV Conclusions

We have shown that the combined spatial and spectral degrees of freedom of the
nonlinear susceptibility in a thin medium lead to intuitive effects on the
free space propagation properties of the generated FWM beam. In particular, we
analyzed a situation where two nonlinear signals are induced by the same
incident beams, and studied their transverse and longitudinal characteristics.
The FWM fields were calculated as superpositions of paraxial modes, with the
expansion coefficients given by the overlap integral of the nonlinear
polarization. Effects of incident beam power and detunings from resonance were
investigated.

For increasing power, we showed that the distribution of coefficient
amplitudes is sensibly affected, leading to greater contributions from higher
radial order modes. The intensity profile at the medium exit was shown to
suffer notable change, with an intensity ring being formed as the pump power
increases. Physically this can be understood as a spatial saturation effect.
This intensity ring is not stable, and under free propagation to distances of
the order of the Rayleigh range outside the interaction medium, the ring-
shaped profile transitions into one that is Gaussian at the center.

When the detunings are varied around resonance, the distribution of phases of
the superimposed modes is mainly affected. In this case, the FWM beam is shown
to undergo focusing or defocusing, depending on whether the fields are above
or below resonance. We argued that it is possible to identify a nonlinear
contribution to the refractive index inside the atomic medium, similar to the
Kerr effect. The fundamental difference is that this nonlinear index is
affected by the intensity distribution of the strong pump beams, and not by
the generated beam itself. Moreover, in the configuration where two FWM
signals are generated simultaneously, we showed that the effect can be either
symmetric between the two beams, if the detunings are equal, or anti-
symmetric, if the detunings are different, e.g., one of the detunings fixed on
resonance.

The transfer of transverse structure of light beams in nonlinear processes is
strongly dictated by the overlap of incident beams, and this assumption leads
to remarkably accurate theoretical predictions both in third- and second-order
processes. However, by taking into account the spatial structure of the medium
susceptibility, effects due to pump power and frequency arise, and we show
that in this case the overall spatial shape of the generated signal can be
sensibly affected. The results presented in this work indicate that the FWM
beam structure, in the near-field, can reveal characteristics of the nonlinear
light-atom interaction, such as the resonances of the medium, saturation
effects and nonlinear variations of the refractive index; while after
propagation to distances of the order of a Rayleigh range, this information
fades away as the generated beam evolves into a stable form, which is strongly
dictated by the overlap of incident beam modes.

##  Appendix A Steady state OBEs

The modeling of the atomic quantities is based on that of Ref. [7], where the
usual case of the nonlinear process involving a strong pump and a weak probe
is considered. In our specific setup both fields possess similar intensities,
and thus we look for solutions to the relevant matrix elements of the density
operator that are influenced by both pumps. The set of steady state equations
for the Fourier components in Eqs. (6) and (7) is

$\displaystyle(\Delta\rho)^{a-b}$
$\displaystyle=\dfrac{2i\widetilde{\Omega}_{a}^{*}\sigma_{12}^{2a-b}+2i\widetilde{\Omega}_{b}^{*}\sigma_{12}^{a}-2i\widetilde{\Omega}_{a}\sigma_{12}^{b*}}{\left(i\delta_{a}-i\delta_{b}+\Gamma\right)},$
(29) $\displaystyle(\Delta\rho)^{\mathrm{dc}}$
$\displaystyle=(\Delta\rho)^{0}-\dfrac{4}{\Gamma}(\mathrm{Im}[\Omega_{a}^{*}\sigma_{12}^{a}]+\mathrm{Im}[\Omega_{b}^{*}\sigma_{12}^{b}]),$
(30) $\displaystyle\sigma_{12}^{a}$
$\displaystyle=\dfrac{i\widetilde{\Omega}_{a}(\Delta\rho)^{\mathrm{dc}}+i\widetilde{\Omega}_{b}(\Delta\rho)^{a-b}}{\left(i\delta_{a}+\Gamma/2\right)},$
(31) $\displaystyle\sigma_{12}^{b}$
$\displaystyle=\dfrac{i\widetilde{\Omega}_{b}(\Delta\rho)^{\mathrm{dc}}+i\widetilde{\Omega}_{a}\left[(\Delta\rho)^{a-b}\right]^{*}}{\left(i\delta_{b}+\Gamma/2\right)},$
(32) $\displaystyle\sigma_{12}^{2a-b}$
$\displaystyle=\dfrac{i\widetilde{\Omega}_{a}(\Delta\rho)^{a-b}}{\left(2i\delta_{a}-i\delta_{b}+\Gamma/2\right)}.$
(33)

With a direct substitution method, it is possible to arrive at the slow
coherence $\sigma_{12}^{2a-b}$ shown in the main text. The coherence related
to signal $S_{2}$, in direction $(2\mathbf{k}_{b}-\mathbf{k}_{a})$, can be
found by introducing in Eq. (6) a term $\rho_{12}^{2b-a}$ that oscillates at
$(2\omega_{b}-\omega_{a})$ and following the same procedure. It is evident
that this will result in a solution to the slow coherence $\sigma^{2b-a}_{12}$
that has the same form as $\sigma_{12}^{2a-b}$ with the exchange of labels
$a\leftrightarrow b$.

##  Appendix B Calculation of the rms parameters

Following Ref. [34] we can write the rms radius of field
$\mathcal{E}_{s}=\sum_{p}\mathcal{A}_{p}u_{p}$ as

$\begin{split}r_{\mathrm{rms}}(z)&=\frac{w(z)}{\sqrt{2}}\left[1+\braket{N}-\mathrm{Re}\left\\{\varphi
e^{2i\tan^{-1}(z/z_{R})}\right\\}\right]^{\frac{1}{2}},\end{split}$

where $\braket{N}=\sum_{p}{\eta}_{p}2p$ is the mean value of the mode order
$N_{p}=2p$ in the superposition and
$\varphi=\frac{1}{P}\sum_{p>0}2p\mathcal{A}_{p}\mathcal{A}_{p-1}^{*}$ is a
generally complex-valued factor. By straightforward manipulations we can
obtain the form shown in Eq. (27) of the main text,
$r_{\mathrm{rms}}(z)=\sqrt{r^{2}_{m}+\theta^{2}_{\mathrm{rms}}(z-z_{m})^{2}}$,
where we have explicitly

$\displaystyle r^{2}_{m}$
$\displaystyle=\frac{w^{2}_{o}}{2}\frac{[1+\braket{N}]^{2}-|\varphi|^{2}}{1+\braket{N}+\mathrm{Re}\\{\varphi\\}},$
(34) $\displaystyle z_{m}$
$\displaystyle=-z_{R}\frac{\mathrm{Im}\\{\varphi\\}}{1+\braket{N}+\mathrm{Re}\\{\varphi\\}},$
(35) $\displaystyle\theta_{\mathrm{rms}}$
$\displaystyle=\dfrac{w_{o}}{\sqrt{2}z_{R}}\left[1+\braket{N}+\mathrm{Re}\\{\varphi\\}\right]^{\frac{1}{2}}.$
(36)

###### Acknowledgements.

This work was supported by CAPES (PROEX 534/2018, No. 23038.003382/2018-39).
The authors would like to acknowledge the financial support from Brazilian
agencies CAPES and CNPq. M. R. L. da Motta acknowledges financial support from
CNPq (130306/2020-7) and CAPES (88887.623521/2021-00). A. A. C. de Almeida
acknowledges financial support by CNPq (141103/2019-1).

## References

  * Kogelnik and Li [1966] H. Kogelnik and T. Li, Appl. Opt. 5, 1550 (1966).
  * Allen  _et al._ [1992] L. Allen, M. W. Beijersbergen, R. J. C. Spreeuw, and J. P. Woerdman, Phys. Rev. A 45, 8185 (1992).
  * Kimel and Elias [1993] I. Kimel and L. R. Elias, IEEE J. Quantum Electron. 29, 2562 (1993).
  * Bandres and Gutiérrez-Vega [2004] M. A. Bandres and J. C. Gutiérrez-Vega, J. Opt. Soc. Am. A 21, 873 (2004).
  * Forbes  _et al._ [2021] A. Forbes, M. de Oliveira, and M. Dennis, Nat. Photonics 15, 253 (2021).
  * Harter and Boyd [1980] D. Harter and R. Boyd, IEEE J. Quantum Electron. 16, 1126 (1980).
  * Boyd  _et al._ [1981] R. W. Boyd, M. G. Raymer, P. Narum, and D. J. Harter, Phys. Rev. A 24, 411 (1981).
  * Yariv and Pepper [1977] A. Yariv and D. M. Pepper, Opt. Lett. 1, 16 (1977).
  * Abrams and Lind [1978] R. L. Abrams and R. C. Lind, Opt. Lett. 2, 94 (1978).
  * Fleischhauer  _et al._ [2005] M. Fleischhauer, A. Imamoglu, and J. P. Marangos, Rev. Mod. Phys. 77, 633 (2005).
  * Chopinaud  _et al._ [2018] A. Chopinaud, M. Jacquey, B. V. de Lesegno, and L. Pruvost, Phys. Rev. A 97, 063806 (2018).
  * Walker  _et al._ [2012] G. Walker, A. Arnold, and S. Franke-Arnold, Phys. Rev. Lett. 108, 243601 (2012).
  * Akulshin  _et al._ [2015] A. M. Akulshin, R. J. McLean, E. E. Mikhailov, and I. Novikova, Opt. Lett. 40, 1109 (2015).
  * Offer  _et al._ [2018] R. F. Offer, D. Stulga, E. Riis, S. Franke-Arnold, and A. S. Arnold, Commun. Phys. 1, 1 (2018).
  * Offer  _et al._ [2020] R. F. Offer, A. Daffurn, E. Riis, P. F. Griffin, A. S. Arnold, and S. Franke-Arnold, Phys. Rev. A 103.2 (2020).
  * Tabosa and Petrov [1999] J. W. R. Tabosa and D. Petrov, Phys. Rev. Lett. 83, 4967 (1999).
  * Barreiro and Tabosa [2003] S. Barreiro and J. W. R. Tabosa, Phys. Rev. Lett. 90, 133001 (2003).
  * Barreiro  _et al._ [2004] S. Barreiro, J. W. R. Tabosa, J. P. Torres, Y. Deyanova, and L. Torner, Opt. Lett. 29, 1515 (2004).
  * Moretti  _et al._ [2009] D. Moretti, D. Felinto, and J. W. R. Tabosa, Phys. Rev. A 79, 023825 (2009).
  * Ding  _et al._ [2013] D. S. Ding, Z. Y. Zhou, B. S. Shi, and G. C. Guo, Nat. Commun. 4, 1 (2013).
  * Lanning  _et al._ [2017] R. N. Lanning, Z. Xiao, M. Zhang, I. Novikova, E. E. Mikhailov, and J. P. Dowling, Phys. Rev. A 96 (2017).
  * Mallick and Dey [2020] N. S. Mallick and T. N. Dey, J. Opt. Soc. Am. B 37, 1857 (2020).
  * Hamedi  _et al._ [2018] H. R. Hamedi, J. Ruseckas, and G. Juzeliūnas, Phys. Rev. A 98, 013840 (2018).
  * Steck [2001] D. A. Steck, Rubidium 87 d line data, https://steck.us/alkalidata/rubidium87numbers.pdf (2001), accessed: 05-05-2022.
  * Jackson [1999] J. D. Jackson,  _Classical electrodynamics_ (Wiley, 1999).
  * Karimi  _et al._ [2014] E. Karimi, R. Boyd, P. De La Hoz, H. De Guise, J. Řeháček, Z. Hradil, A. Aiello, G. Leuchs, and L. L. Sánchez-Soto, Phys. Rev. A 89, 063813 (2014).
  * Plick and Krenn [2015] W. N. Plick and M. Krenn, Phys. Rev. A 92, 063841 (2015).
  * Pereira  _et al._ [2017] L. Pereira, W. Buono, D. Tasca, K. Dechoum, and A. Khoury, Phys. Rev A 96, 053856 (2017).
  * Wu  _et al._ [2020] H. J. Wu, L. W. Mao, Y. J. Yang, C. Rosales Guzmán, W. Gao, B. S. Shi, and Z. H. Zhu, Phys. Rev. A 101, 063805 (2020).
  * Schwob  _et al._ [1998] C. Schwob, P. F. Cohadon, C. Fabre, M. Marte, H. Ritsch, A. Gatti, and L. Lugiato, Appl. Phys. B 66, 685 (1998).
  * Alves  _et al._ [2018] G. Alves, R. Barros, D. Tasca, C. Souza, and A. Khoury, Phys. Rev. A 98, 063825 (2018).
  * Buono  _et al._ [2020] W. Buono, A. Santos, M. Maia, L. Pereira, D. Tasca, K. Dechoum, T. Ruchon, and A. Khoury, Phys. Rev. A 101, 043821 (2020).
  * Phillips and Andrews [1983] R. L. Phillips and L. C. Andrews, Appl. Opt. 22, 643 (1983).
  * Vallone  _et al._ [2016] G. Vallone, G. Parisi, F. Spinello, E. Mari, F. Tamburini, and P. Villoresi, Phys. Rev. A 94 (2016).
  * Boyd [2020] R. Boyd,  _Nonlinear Optics_ (Academic Press, 2020).
  * Yu and Wang [2021] C. Yu and Z. Wang, Phys. Rev. A 103, 013518 (2021).


[a]Andreas Ringwald

# Prospects to scrutinise or smash SM*A*S*H

###### Abstract

SM*A*S*H is an extension of the Standard Model of particle physics which has
just the minimal number of fields in order to solve six puzzles of particle
physics and cosmology in one smash: vacuum stability, inflation, baryon
asymmetry, neutrino masses, strong CP, and dark matter. The parameters of
SM*A*S*H are constrained by symmetries and requirements to solve these
puzzles. This provides various firm predictions for observables which can be
confronted with experiments. We discuss the prospects and timeline to
scrutinise or smash SM*A*S*H by cosmic microwave background polarisation
experiments, axion haloscopes, and future space-borne gravitational wave
detectors.

##  1 SM*A*S*H

The acronym “SM*A*S*H” stands for “Standard Model*Axion*Seesaw*Higgs-Portal
Inflation” – a minimal extension of the Standard Model (SM) which solves six
puzzles of particle physics and cosmology in one smash [1, 2]: vacuum
stability, inflation, baryon asymmetry, neutrino masses, strong CP, and dark
matter.111Similar models have been considered in Refs. [3, 4, 5, 6, 7].

$\begin{array}[]{|c|c|c|c|c|c|c|c|c|}\hline\cr
q&u&d&L&N&E&Q&\tilde{Q}&\sigma\\\ \hline\cr
1/2&-1/2&-1/2&1/2&-1/2&-1/2&-1/2&-1/2&1\\\\[2.15277pt] \hline\cr\end{array}$

Table 1: PQ-charge assignments of the fields in SM*A*S*H. The remaining SM
fields have no PQ charge.

A SM-singlet complex scalar field $\sigma$ (the Peccei-Quinn (PQ) field), a
vector-like quark $Q$ and three SM-singlet neutrinos $N_{i}$, with $i=1,2,3$,
are added to the SM. All the new fields, as well as the quarks and leptons of
the SM, are assumed to be charged under a global $U(1)_{\rm PQ}$ symmetry, cf.
Table 1. The scalar potential respecting this PQ symmetry has the general
form:
$V(H,\sigma)=\lambda_{H}\left(H^{\dagger}H-\frac{v^{2}}{2}\right)^{2}+\lambda_{\sigma}\left(|\sigma|^{2}-\frac{v_{\sigma}^{2}}{2}\right)^{2}+2\lambda_{H\sigma}\left(H^{\dagger}H-\frac{v^{2}}{2}\right)\left(|\sigma|^{2}-\frac{v_{\sigma}^{2}}{2}\right)$,
where $H$ is the SM Higgs doublet. In order to ensure that both the PQ and the
electroweak symmetry are broken by the vacuum expectation values $\langle
H^{\dagger}H\rangle=v^{2}/2$, $\langle|\sigma|^{2}\rangle=v_{\sigma}^{2}/2$,
where $v_{\sigma}\gg v\simeq 246$ GeV, the scalar couplings are required to
obey $\lambda_{H},\lambda_{\sigma}>0$,
$\lambda_{H\sigma}^{2}<\lambda_{H}\lambda_{\sigma}$. The $U(1)_{Y}$
hypercharge of the vector-like quark $Q$ is required to be $-1/3$, such that
the most general Yukawa interactions of $Q$ and $N_{i}$, allowed by SM gauge
and PQ symmetries, are ${\cal L}\supset-[F_{ij}\bar{N}_{j}P_{L}L_{i}\epsilon
H+\frac{1}{2}Y_{ij}\sigma\bar{N}_{i}P_{L}N_{j}+y\,\sigma\bar{Q}P_{L}Q+\,{y_{Q_{d}}}_{i}\sigma\bar{D}_{i}P_{L}Q+h.c.]$,
where $D_{i}$, $L_{i}$ denote the Dirac spinors associated with the down
quarks and leptons of the $i$th generation, while the $N_{i}$ are taken to be
Majorana spinors. Electroweak vacuum instability [8] – the instability of the
Higgs potential at large field values, present for the preferred value of the
top mass – can be avoided in SM*A*S*H by the stabilizing effect of the Higgs
portal coupling $\lambda_{H\sigma}$ [9, 10]. This requires
$\lambda^{2}_{H\sigma}/\lambda_{\sigma}$ to be between $\sim 10^{-2}$ and
$\sim 10^{-1}$ [2]. Inflation is realised in SM*A*S*H via non-minimal chaotic
inflation. The inflaton is comprised by a mixture of the real part of the PQ
field and the modulus of the Higgs. Reheating of the Universe after inflation
proceeds then efficiently via the Higgs portal. SM*A*S*H predicts a lower
bound on the ratio of the power in tensor to scalar fluctuations, $r\gtrsim
0.004$ [1, 2], a reheating temperature around $10^{12}$ GeV [11], and a second
order PQ phase transition at around $T_{c}\sim 10^{8}$ GeV [1, 2]. The strong
CP puzzle is solved in SM*A*S*H by the PQ mechanism [12]. The axion [13, 14] –
the pseudo Goldstone boson associated with the spontaneous breaking of the PQ
symmetry – is the prime candidate for cold dark matter in SM*A*S*H [15, 16,
17]. Requiring the latter to constitute 100% of dark matter, its decay
constant, $f_{a}=v_{\sigma}$, is predicted to be in the range $10^{10}\,{\rm
GeV}\lesssim f_{a}<2\times 10^{11}\,\mathrm{GeV}$. The PQ symmetry breaking
scale also gives rise to large Majorana masses for the heavy neutrinos. This
can explain the smallness of the masses of the active neutrions through the
seesaw mechanism [18, 19, 20, 21] and results in the generation of the baryon
asymmetry of the universe via thermal leptogenesis [22].

In summary: The parameters in SM*A*S*H are constrained by symmetries and
requirements to solve the puzzles of particle physics and cosmology,
respectively. This provides various firm predictions of observables which can
be confronted with experiments. In this contribution I will focus on some of
these predictions. The order of discussion follows the approximate timeline of
experiments being able to probe the predictions from SM*A*S*H.

Figure 1: Inflationary constraints/predictions on $\tilde{\lambda}_{\sigma}$
(left), Hubble scale at the beginning and end of inflation (center) and field
value at the end of inflation (right) as a function of $\xi_{\sigma}$ and for
$k_{*}=0.002$ Mpc-1 [11]. The blue regions are compatible within 95% C.L. with
the latest combination of Planck and BICEP/KECK data [30]. The red lines
correspond to the predictions when taking into account the SM*A*S*H prediction
of radiation domination immediately after inflation, and the red dots
correspond to the benchmark scenarios BP1 and BP2. The gray shaded regions to
the right of $\xi_{\sigma}>1$ indicate the region where the predictions may be
threatened by the breakdown of peturbative unitarity. Figures adapted from
Ref. [11].

##  2 Confronting SM*A*S*H with Cosmic Microwave Background Observations

Inflation is realised in SM*A*S*H by the dynamics of the PQ and Higgs fields
in the presence of their generically present non-minimal gravitational
couplings to the Ricci scalar $R$ [23, 24, 25, 26, 27], described, in the
Jordan frame, by the action $S\supset-\int
d^{4}x\sqrt{-g}\,\left[\frac{M^{2}}{2}+\xi_{H}\,H^{\dagger}H+\xi_{\sigma}\,\sigma^{*}\sigma\right]R$.
The mass scale $M$ is related to the reduced Planck mass ($M_{P}\simeq
2.435\times 10^{18}\,\rm GeV$) by
$M^{2}_{P}=M^{2}+\xi_{H}v^{2}+\xi_{\sigma}v^{2}_{\sigma}$. After a Weyl
transformation of the metric to the Einstein frame, the scalar potential
becomes flat for large field values. Eventual problems with perturbative
unitarity [28, 29] are avoided by requiring
$1\gtrsim\xi_{\sigma}\gg\xi_{H}\geq 0$. For positive $\lambda_{H\sigma}$,
inflation can take place along the direction of $\phi\equiv\sqrt{2}{\rm
Re}\,\sigma$, but in this case reheating can be shown to be problematic,
leading to an excess of dark radiation [1, 2]. In the case of a negative Higgs
portal coupling, $\lambda_{H\sigma}<0$, slow-roll inflation takes place along
a valley in the scalar potential that can be approximated by the line
$h/\phi=\sqrt{-\lambda_{H\sigma}/\lambda_{H}}$, where $h$ is the neutral
component of the Higgs doublet in the unitary gauge. The potential along the
inflaton valley,
$\tilde{V}(\chi)=\frac{1}{4}\tilde{\lambda}_{\sigma}\phi(\chi)^{4}\left(1+\xi_{\sigma}\frac{\phi(\chi)^{2}}{M_{P}^{2}}\right)^{-2}$,
is determined by two model parameters: an effective quartic coupling
$\tilde{\lambda}_{\sigma}=\lambda_{\sigma}-\lambda_{H\sigma}^{2}/\lambda_{H}$
and the non-minimal coupling $\xi_{\sigma}$. Here, $\phi$ and the canonically
normalized inflaton field $\chi$ are related by
$\Omega^{2}\,d\chi/d\phi\simeq(b\,\Omega^{2}+6\,\xi_{\sigma}^{2}\,\phi^{2}/M_{P}^{2})^{1/2}$,
with $\Omega^{2}=1+\xi_{\sigma}\frac{\phi(\chi)^{2}}{M_{P}^{2}}$ and
$b=1+|\lambda_{H\sigma}/\lambda_{H}|$. Quantum fluctuations during slow-roll
inflation along this potential produce power spectra of scalar metric
perturbations (density waves) and tensor metric perturbations (gravitational
waves (GWs)) which can be parametrized as
$\Delta^{2}_{s/t}(k)=A_{s/t}(k_{*})\left({k}/{k_{*}}\right)^{n_{s/t}(k_{*})-1+\cdots}$,
where $k_{*}$ is a given reference pivot scale. Fitting the amplitude of these
perturbations inferred from the observed cosmic microwave background (CMB)
temperature and polarisation maps imposes one relation between these two
parameters $\tilde{\lambda}_{\sigma}$, $\xi_{\sigma}$, which is approximately
given by $7\times 10^{-3}\lesssim\xi_{\sigma}\simeq 4\times
10^{4}\sqrt{\tilde{\lambda}_{\sigma}}\lesssim 1$, see Fig. 1 (left).
Accordingly, quantities during inflation can be characterized by a single
parameter, which can be chosen as $\xi_{\sigma}$, as illustrated in Fig. 1
(middle and right). The numerical predictions in the $r$ vs $n_{s}$ plane are
shown in Fig. 2 for $k_{*}=0.002$ Mpc-1. The thin solid lines indicate the
values of the non-minimal coupling, bounded by the quartic chaotic inflation
limit, $\xi_{\sigma}=0$, and the $\xi_{\sigma}\rightarrow\infty$ limit. The
dashed lines show the number of e-folds, $N=\Delta\log a$, of the cosmic scale
factor $a$ between the pivot scale’s crossing of the horizon and the end of
inflation, cf. Fig. 3. The figure also shows the 68% and 95% C.L. regions from
the combined data analysis of Planck and BICEP/KECK [30] at $k_{*}=0.002$
Mpc-1. We see, that $\xi_{\sigma}>10^{-2}$ is perfectly compatible with
current bounds on $n_{s}$ and $r$ and it also gives an adequate $N$.

Figure 2:  Inflationary predictions in SM*A*S*H in the $r$ vs $n_{s}$ plane
for $k_{*}=0.002$ Mpc-1. The green solid/dashed gray lines are contours of
constant $\xi_{\sigma}$/number of efolds, respectively. Accounting for the
post-inflationary expansion history gives the orange region, and the red dots
correspond to the benchmark scenarios BP1 and BP2. We also show the 68% and
95% C.L. contours arising from Planck and BICEP/KECK data [30], as well as the
95% projected sensitivities from BICEP Array [30], the Simons observatory
[31], LiteBird [32], and CMB-S4 [33]. Figure adapted from Ref. [11].

The latter can be determined by using the post-inflationary evolution of the
universe to match the scales of current CMB perturbations to their values at
horizon crossing during inflation. Inflation in SM*A*S*H ends when
$\phi\sim\mathcal{O}(M_{P})$, see Fig. 1 (right), after which the background
goes through Hubble-damped oscillations that have the equation of state of a
radiation bath. Hence radiation domination starts immediately after inflation,
which fixes $N$ to around $60$, as can be inferred from Fig. 3. This narrows
the SM*A*S*H prediction to the thick orange line in Fig. 2. This prediction
will be probed exhaustively in the next decade by CMB polarisation
experiments, such as BICEP Array [30], the Simons observatory [31], LiteBird
[32], and CMB-S4 [33], see Fig. 2. We conclude: SM*A*S*H is smashed if CMB-S4
or LiteBird do not discover B modes from inflationary tensor metric
perturbations (GWs)!

Figure 3: The expansion and thermal history of the Universe in SM*A*S*H,
emphasising the transition from inflation to radiation-domination-like
expansion $a{\cal H}\propto 1/a$ before standard matter and cosmological
constant domination epochs. Figure adapted from Ref. [2].

##  3 Scrutinising SM*A*S*H with Axion Dark Matter Experiments

For $\lambda_{H\sigma}<0$, the oscillations of the inflaton field after
inflation allow for an efficient reheating. They result in the copious
production of scalar field fluctuations via parametric resonance
(“preheating”), thereby restoring the PQ symmetry after a few oscillations.
The decay of the Higgs component of the inflaton into SM particles (mainly top
quarks) reheats the Universe efficiently. The reheating temperature has been
determined via lattice simulations of the evolution of the scalar fields in a
an expanding Friedmann-Robertson-Walker background to be in the range $T_{\rm
rh}\approx 10^{12-13}\,{\rm GeV}$ [11]. The predicted range is quite narrow,
since the parameter space in SM*AS*H for the bosonic couplings of the PQ-field
$\sigma$ is significantly constrained by requiring consistency of the
inflationary predictions with CMB observations and of vacuum stability in the
Higgs direction. Roughly, a given tensor-to-scalar ratio $r$ fixes
$\xi_{\sigma}$ (see Fig. 2) which determines the effective quartic coupling
$\tilde{\lambda}_{\sigma}$ (see Fig. 1 (left)). From the stability constraints
it follows that $\tilde{\lambda}_{\sigma}$ cannot be very different from
$\lambda_{\sigma}$. This means that choosing $r$ roughly specifies all the
bosonic couplings of $\sigma$, which then settles the scalar field dynamics
which determines preheating and reheating. For benchmark point BP1 in Fig. 1,
the simulations yielded a reheating temperature $T_{\rm rh}=9.7\times 10^{12}$
GeV, while for BP2 in the same figure, they resulted in $T_{\rm rh}=2.0\times
10^{12}$ GeV [11].

Figure 4: Electromagnetic coupling of the axion versus its mass. Haloscope
exclusion regions (filled blue) and projected sensitivities (dashed blue
lines) assume the axion to be 100% of the halo dark matter. The solid orange
band denotes the parameter region where the axion in SM*A*S*H is predicted to
constitute 100% of the cold dark matter in the Universe. Figure adapted from
Ref. [35].

The axion is “born” when the PQ symmetry, which was restored during
preheating, is broken again during the radiation dominated hot phase after
reheating. The PQ phase transition is predicted to be second order, with
critical temperature
$T_{c}\simeq{2\sqrt{6\lambda_{\sigma}}v_{\sigma}}/{\sqrt{8(\lambda_{\sigma}+\lambda_{H\sigma})+\sum_{i}Y_{ii}^{2}+6y^{2}}}\simeq\lambda_{\sigma}^{1/4}v_{\sigma}\simeq(2-8)\times
10^{-3}\,v_{\sigma}$ [2], where the last two estimates take into account the
parameter constraints arising from the requirements of vacuum stability and
consistency of the inflationary predictions with CMB observations,
respectively. After the PQ phase transition, the Universe consists of many
causally disconnected patches, with sizes of order the Hubble radius, in which
the axion field $\theta=a/f_{a}$ relaxes quickly to a uniform value between
$-\pi$ and $\pi$. At the boundaries of these patches, one-dimensional
topological defects - cosmic axion strings - are formed. The corresponding
string loop network evolves during the Universe’s expansion by tightening,
oscillating, reconnecting, and collapsing, thereby continuously producing a
population of effectively massless axions. This continues until the
temperature of the Universe drops to around a GeV, when the Hubble expansion
rate drops below the axion mass. Thereafter, axions are also produced from the
decay of two-dimensional topological defects - axion domain walls - and from
the realignment mechanism: within any Hubble size patch, the axion field
starts to oscillate around the minimum of its potential, behaving as a cold
dark matter fluid with equation of state $w_{a}=p_{a}/\rho_{p}\approx 0$.

A conservative upper bound on the PQ symmetry breaking scale,
$f_{a}=v_{\sigma}<2.3\times 10^{11}\,{\rm GeV}$, can be set by requiring that
the predicted contribution to the dark matter abundance from the realignment
mechanism alone does not exceed the observed dark matter abundance. This
translates directly into a lower bound on the axion mass, $m_{a}\simeq 5.9\
{\rm\mu eV}\left(\frac{10^{12}\,{\rm GeV}}{f_{a}}\right)>26\ \mu{\rm eV}$
[34]. Comparing with the sensitivity prospects of axion dark matter
experiments (axion haloscopes) shown in Fig. 4, we conclude: ADMX [36],
BabyIAXO-RADES [37], CAPP [38], FLASH [39], or DMRadio [40] have good
prospects to smash SM*A*S*H by discovering an axion with a mass below 26
micro-eV!

Unfortunately, the prediction of the axion dark matter contribution from
cosmic strings and domain walls suffers from significant theoretical
uncertainties which arise from the difficulty in quantifying precisely the
energy loss processes of topological defects and the generated axion spectra.
In fact, the results from present state-of-the-art first principle classical
field theory simulations [41, 42, 43, 44, 45, 46, 47] still allow for two
extreme possibilities: the contribution from the decay of topological defects
may either be negligible [42] or overwhelming [46] in comparison to the one
from the realignment mechanism. In the latter case, a conservative lower bound
on $f_{a}=v_{\sigma}\gtrsim 1.1\times 10^{10}\,{\rm GeV}$, corresponding to an
upper bound on $m_{a}\lesssim 0.5\ {\rm meV}$, is obtained from the
requirement that the produced axions constitute 100% of the cold dark matter
in the Universe. We infer from Fig. 4 that ALPHA [48] and MADMAX [49] have
good prospects to detect the axion from SM*A*S*H if the latter constitutes
100% of the dark matter in the galactic halo.222Large density variations in
the initial state of the axion field after the PQ phase transition may lead to
the formation of compact dark matter objects known as “miniclusters”. This
leads to an increased theoretical uncertainty in the local axion density for
dark matter detection in the laboratory.

##  4 Searching for SM*A*S*H Signatures in the Cosmic GW Background

In this section, we are looking into the prospects to scrutinise SM*A*S*H in
the far future, say in the second half of this century.

Suppose, the CMB polarisation experiments have discovered, in the present
decade, the B modes induced by the primordial tensor metric fluctuations
originating from quantum fluctuations during inflation. This will result in a
strong push towards the realisation of space-borne GW laser interferometers
which are designed to search directly for the GWs from inflation, such as the
Big Bang Observer (BBO) [50, 51, 52, 53] or the Deci-Hertz Interferometer
Gravitational Wave Observatory (DECIGO) [54, 55]. In fact, in a frequency
range around a Hz, these detectors have the prospected sensitivity to discover
the stochastic GW background predicted from inflation in SM*A*S*H [56], as can
be seen from Fig. 5.

Figure 5:  Characteristic amplitude of primordial GWs in SM*A*S*H (orange;
upper (lower) curve corresponding to the benchmark BP1 (BP2) in Fig. 2)
compared to present (shaded areas) and projected sensitivities (colored solid
lines) [57, 54, 58, 59, 50, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71,
72]. Indirect dark radiation constraints [73, 74, 75] are shown with dashed
lines. Figure adapted from Ref. [11].

After a successful detection, one may envisage an upgrade of the detectors,
such that their noise curves are determined by the quantum limit, for example
to an “ultimate” DECIGO [60]. The latter will enable a precise measurement of
the shape of the spectrum in this frequency range, see Fig. 5. Intriguingly,
SM*A*S*H predicts a step in the spectrum within the sensitivity frequency
window of ultimate DECIGO, see Fig. 6 (left). This step is an observational
signature of the second order PQ phase transition in SM*A*S*H, resulting in a
drastic change in the equation of state of the thermal plasma, notably in the
number of relativistic degrees of freedom, $g_{\ast\rho}(T)$, when the
temperature drops below the critical temperature, $2\times
10^{7}\,\mathrm{GeV}\lesssim
T_{c}\simeq\lambda_{\sigma}^{1/4}\,v_{\sigma}\lesssim 2\times
10^{9}\,\mathrm{GeV}$,333Here, we have used for the range of $v_{\sigma}$ the
one in which the axion can be 100% of the cold dark matter, see Fig. 6
(right). Entropy conservation implies that the growth rate of the Hubble
radius ${\mathcal{H}}^{-1}$ is diminished during the PQ phase transition.
Thus, the rate at which primordial tensor modes cross into the horizon is
changed during the transition, resulting in a step in the spectrum at
frequencies corresponding to the Hubble rate at the transition,
$f=\frac{a(T_{c}){\mathcal{H}}(Tc)}{2\pi
a_{0}}\approx\,\mathrm{Hz}\left(\frac{T_{c}}{10^{8}\,\mathrm{GeV}}\right)$.

Figure 6:  Left panel: The spectrum of primordial GWs predicted in SM*A*S*H
for BP2 and $f_{a}=v_{\sigma}=1.2\times 10^{11}\,\mathrm{GeV}$ (solid blue
line), corresponding to $T_{c}\simeq 5\times 10^{9}$ GeV, in comparison with
the results obtained by using the effective relativistic degrees of freedom in
the SM. Dark yellow bars correspond to the prospected sensitivity of ultimate
DECIGO. Right panel: Temperature dependence of $g_{*\rho}(T)$ for some
specific choice of the parameters in SM*A*S*H (similar to BP1 with
$f_{a}=v_{\sigma}=1.2\times 10^{11}\,\mathrm{GeV}$) that leads to $T_{c}\simeq
5\times 10^{8}$ GeV. Shaded regions represent the uncertainty due to the
choice of the renormalization scale and estimate of the QCD corrections. The
gray dashed line corresponds to the critical temperature of the PQ phase
transition. Both figures adapted from Ref. [56].

Clearly, the PQ scale $v_{\sigma}$ is the crucial parameter determining the
spectrum of the primordial GWs from inflation in SM*A*S*H, since it affects
$T_{c}\propto v_{\sigma}$ and therefore changes the frequency range at which
the step appears. However, in second half of this century, when ultimate
DECIGO may operate, and in the case that the axion constitutes the main
component of cold dark matter, axion haloscopes such as ALPHA and MADMAX
should have detected the axion already and determined its mass with high
precision, from which one can infer its decay constant, which is, in SM*A*S*H,
equal to the PQ scale, $f_{a}=v_{\sigma}$. Therefore, a detection of the step
in the GW spectrum will yield further information about SM*A*S*H, such as the
quartic PQ field coupling $\lambda_{\sigma}$ and the Higgs portal coupling
$\lambda_{H\sigma}$, which can be cross-correlated with the one obtained from
CMB polarisation observations. In this sense, future GW observations of the
GWs from inflation can be used to probe the details of the PQ sector that may
not be reached in other high-energy experiments.

But this is not the end of the story. SM*A*S*H has also definite predictions
for stochastic GW backgrounds from other epochs of its cosmological history
[11], notably from inflaton fragmentation during preheating and from the
thermal plasma at the beginning of the hot thermal radiation-dominated stage,
see Fig. 5. The contributions from the different epochs are not independent
and their features are correlated: each epoch determines the initial
conditions for the subsequent one. A hypothetical detection of the complete
spectrum in different frequency ranges would allow to cross-check for the
correlations predicted in SM*A*S*H, opening new possibilities for falsifying
the model. For example, a hypothetical measurement of the GW spectrum between
$0.1$ MHz and $100$ GHz could potentially determine the scale of inflaton
fragmentation after inflation, related to the peak frequency of the GWs from
preheating, and the maximum temperature and the number of relativistic degrees
of freedom of the hot Big Bang plasma, which fix the amplitude and peak of the
GWs from the thermal plasma [76, 75, 70]. This could provide an unprecedented
window into the physics of the very early universe. However, probing the GWs
generated by preheating and the thermal plasma requires much progress in the
detection of ultra high frequency GWs, as is obvious from Fig. 5. Fortunately,
a worldwide initiative towards this goal has already started [77].

##  5 I Have a Dream

Let us end this proceeding’s contribution with a pipe dream of mine concerning
the establishment of SM*A*S*H as the theory of particle physics and cosmology.

In the early 2030’s, the CMB polarisation experiments CMB-S4 and LiteBIRD
discover B modes induced by inflationary GWs consistent with the prediction
from SM*A*S*H in Fig. 2. Again in the 2030’s, the axion haloscopes ALPHA and
MADMAX discover axion dark matter in the mass region predicted by SM*A*S*H in
Fig. 4. In the 2060’s, a GW laser interferometer with sensitivity in the Hz
range similar to BBO or DECIGO in Fig. 5 is launched which directly detects
the stochastic GW background predicted from inflation in SM*A*S*H shown in the
same figure. In the 2080’s, an upgrade of the latter space-borne GW detector
similar to ultimate DECIGO discovers the step in the inflationary GW
background in Fig. 6 arising from the PQ phase transition in SM*A*S*H.
Moreover, the SM*A*S*H parameter ranges inferred from the position of the step
is consistent with the parameter ranges inferred from the measurements at CMB
polarisation experiments and axion haloscopes. Of course, this motivates then
to develop detectors also for ultra high frequency GWs in order to probe the
SM*A*S*H predictions in Fig. 5 from preheating and the thermal plasma.
However, I do not dare to give a time scale for this detection because this
experimental field is still in its infancy.

## Acknowledgments

Special thanks to Guillermo Ballesteros, Kenichi Saikawa, and Carlos Tamarit
for valuable comments on the draft. This work has been partially funded by the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy - EXC 2121 Quantum Universe - 390833306 and
under \- 491245950. This article/publication is based upon work from COST
Action COSMIC WISPers CA21106, supported by COST (European Cooperation in
Science and Technology).

## References

  * [1] G. Ballesteros, J. Redondo, A. Ringwald and C. Tamarit, Phys. Rev. Lett. 118 (2017) no.7, 071802 [arXiv:1608.05414 [hep-ph]].
  * [2] G. Ballesteros, J. Redondo, A. Ringwald and C. Tamarit, JCAP 08 (2017), 001 [arXiv:1610.01639 [hep-ph]].
  * [3] S. M. Boucenna, S. Morisi, Q. Shafi and J. W. F. Valle, Phys. Rev. D 90 (2014) no.5, 055023 [arXiv:1404.3198 [hep-ph]].
  * [4] G. Ballesteros, J. Redondo, A. Ringwald and C. Tamarit, Front. Astron. Space Sci. 6 (2019), 55 [arXiv:1904.05594 [hep-ph]].
  * [5] A. H. Sopov and R. R. Volkas, Phys. Dark Univ. 42 (2023), 101381 [arXiv:2206.11598 [hep-ph]].
  * [6] M. Berbig, JCAP 11 (2022), 042 [arXiv:2207.08142 [hep-ph]].
  * [7] J. Dutta, M. Matlis, G. Moortgat-Pick and A. Ringwald, [arXiv:2309.10857 [hep-ph]].
  * [8] G. Degrassi, S. Di Vita, J. Elias-Miro, J. R. Espinosa, G. F. Giudice, G. Isidori and A. Strumia, JHEP 08 (2012), 098 [arXiv:1205.6497 [hep-ph]].
  * [9] O. Lebedev, Eur. Phys. J. C 72 (2012), 2058 [arXiv:1203.0156 [hep-ph]].
  * [10] J. Elias-Miro, J. R. Espinosa, G. F. Giudice, H. M. Lee and A. Strumia, JHEP 06 (2012), 031 [arXiv:1203.0237 [hep-ph]].
  * [11] A. Ringwald and C. Tamarit, Phys. Rev. D 106 (2022) no.6, 063027 [arXiv:2203.00621 [hep-ph]].
  * [12] R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38 (1977), 1440-1443
  * [13] S. Weinberg, Phys. Rev. Lett. 40 (1978), 223-226
  * [14] F. Wilczek, Phys. Rev. Lett. 40 (1978), 279-282
  * [15] J. Preskill, M. B. Wise and F. Wilczek, Phys. Lett. B 120 (1983), 127-132
  * [16] L. F. Abbott and P. Sikivie, Phys. Lett. B 120 (1983), 133-136
  * [17] M. Dine and W. Fischler, Phys. Lett. B 120 (1983), 137-141
  * [18] P. Minkowski, Phys. Lett. B 67 (1977), 421-428
  * [19] M. Gell-Mann, P. Ramond and R. Slansky, Conf. Proc. C 790927 (1979), 315-321 [arXiv:1306.4669 [hep-th]].
  * [20] T. Yanagida, Conf. Proc. C 7902131 (1979), 95-99 KEK-79-18-95.
  * [21] R. N. Mohapatra and G. Senjanovic, Phys. Rev. Lett. 44 (1980), 912
  * [22] M. Fukugita and T. Yanagida, Phys. Lett. B 174 (1986), 45-47
  * [23] B. L. Spokoiny, Phys. Lett. B 147 (1984), 39-43
  * [24] T. Futamase and K. i. Maeda, Phys. Rev. D 39 (1989), 399-404
  * [25] D. S. Salopek, J. R. Bond and J. M. Bardeen, Phys. Rev. D 40 (1989), 1753
  * [26] R. Fakir and W. G. Unruh, Phys. Rev. D 41 (1990), 1783-1791
  * [27] F. L. Bezrukov and M. Shaposhnikov, Phys. Lett. B 659 (2008), 703-706 [arXiv:0710.3755 [hep-th]].
  * [28] J. L. F. Barbon and J. R. Espinosa, Phys. Rev. D 79 (2009), 081302 [arXiv:0903.0355 [hep-ph]].
  * [29] C. P. Burgess, H. M. Lee and M. Trott, JHEP 09 (2009), 103 [arXiv:0902.4465 [hep-ph]].
  * [30] P. A. R. Ade et al. [BICEP and Keck], Phys. Rev. Lett. 127 (2021) no.15, 151301 [arXiv:2110.00483 [astro-ph.CO]].
  * [31] P. Ade et al. [Simons Observatory], JCAP 02 (2019), 056 [arXiv:1808.07445 [astro-ph.CO]].
  * [32] E. Allys et al. [LiteBIRD], PTEP 2023 (2023) no.4, 042F01 [arXiv:2202.02773 [astro-ph.IM]].
  * [33] K. Abazajian et al., [arXiv:1907.04473 [astro-ph.IM]].
  * [34] S. Borsanyi et al., Nature 539 (2016) no.7627, 69-71 [arXiv:1606.07494 [hep-lat]].
  * [35] C. O’Hare, doi:10.5281/zenodo.3932430, https://github.com/cajohare/AxionLimits
  * [36] T. Braine et al. [ADMX], Phys. Rev. Lett. 124 (2020) 101303 [arXiv:1910.08638 [hep-ex]].
  * [37] S. Ahyoune et al., Annalen Phys. 535 (2023) no.12, 2300326 [arXiv:2306.17243 [physics.ins-det]].
  * [38] A. K. Yi et al., Phys. Rev. Lett. 130 (2023) no.7, 071002 [arXiv:2210.10961 [hep-ex]].
  * [39] D. Alesini et al., Phys. Dark Univ. 42 (2023), 101370 [arXiv:2309.00351 [physics.ins-det]].
  * [40] L. Brouwer et al. [DMRadio], Phys. Rev. D 106 (2022) no.10, 103008 [arXiv:2204.13781 [hep-ex]].
  * [41] M. Kawasaki, K. Saikawa and T. Sekiguchi, Phys. Rev. D 91 (2015) no.6, 065014 [arXiv:1412.0789 [hep-ph]].
  * [42] V. B. Klaer and G. D. Moore, JCAP 11 (2017), 049 [arXiv:1708.07521 [hep-ph]].
  * [43] M. Gorghetto, E. Hardy and G. Villadoro, JHEP 07 (2018), 151 [arXiv:1806.04677 [hep-ph]].
  * [44] M. Buschmann, J. W. Foster and B. R. Safdi, Phys. Rev. Lett. 124 (2020) no.16, 161103 [arXiv:1906.00967 [astro-ph.CO]].
  * [45] M. Hindmarsh, J. Lizarraga, A. Lopez-Eiguren and J. Urrestilla, Phys. Rev. Lett. 124 (2020) no.2, 021301 [arXiv:1908.03522 [astro-ph.CO]].
  * [46] M. Gorghetto, E. Hardy and G. Villadoro, SciPost Phys. 10 (2021) no.2, 050 [arXiv:2007.04990 [hep-ph]].
  * [47] M. Buschmann, J. W. Foster, A. Hook, A. Peterson, D. E. Willcox, W. Zhang and B. R. Safdi, Nature Commun. 13 (2022) no.1, 1049 [arXiv:2108.05368 [hep-ph]].
  * [48] M. Lawson, A. J. Millar, M. Pancaldi, E. Vitagliano and F. Wilczek, Phys. Rev. Lett. 123 (2019) no.14, 141802 [arXiv:1904.11872 [hep-ph]].
  * [49] S. Beurthey, N. Böhmer, P. Brun, A. Caldwell, L. Chevalier, C. Diaconu, G. Dvali, P. Freire, E. Garutti and C. Gooch, et al. [arXiv:2003.10894 [physics.ins-det]].
  * [50] S. Phinney et al., “The Big Bang Observer: Direct detection of gravitational waves from the birth of the Universe to the Present”, NASA Mission Concept Study (2004).
  * [51] J. Crowder and N. J. Cornish, Phys. Rev. D 72 (2005) 083005 [gr-qc/0506015].
  * [52] V. Corbin and N. J. Cornish, Class. Quant. Grav. 23 (2006) 2435 [gr-qc/0512039].
  * [53] G. M. Harry, P. Fritschel, D. A. Shaddock, W. Folkner and E. S. Phinney, Class. Quant. Grav. 23 (2006) 4887 Erratum: [Class. Quant. Grav. 23 (2006) 7361].
  * [54] N. Seto, S. Kawamura and T. Nakamura, Phys. Rev. Lett. 87 (2001), 221103 [arXiv:astro-ph/0108011 [astro-ph]].
  * [55] S. Kawamura et al., Class. Quant. Grav. 23 (2006) S125.
  * [56] A. Ringwald, K. Saikawa and C. Tamarit, JCAP 02 (2021), 046 [arXiv:2009.02050 [hep-ph]].
  * [57] B. P. Abbott et al. [LIGO Scientific], Class. Quant. Grav. 34 (2017) no.4, 044001 [arXiv:1607.08697 [astro-ph.IM]].
  * [58] M. Punturo et al., Class. Quant. Grav. 27 (2010), 194002
  * [59] P. Amaro-Seoane et al. [LISA], [arXiv:1702.00786 [astro-ph.IM]].
  * [60] S. Kuroyanagi, K. Nakayama and J. Yokoyama, PTEP 2015 (2015) no.1, 013E02 [arXiv:1410.6618 [astro-ph.CO]].
  * [61] N. Aggarwal, G. P. Winstone, M. Teo, M. Baryakhtar, S. L. Larson, V. Kalogera and A. A. Geraci, Phys. Rev. Lett. 128 (2022) no.11, 111101 [arXiv:2010.13157 [gr-qc]].
  * [62] A. S. Chou et al. [Holometer], Phys. Rev. D 95 (2017) no.6, 063002 [arXiv:1611.05560 [astro-ph.IM]].
  * [63] M. Goryachev and M. E. Tobar, Phys. Rev. D 90 (2014) no.10, 102005 [erratum: Phys. Rev. D 108 (2023) no.12, 129901] [arXiv:1410.2334 [gr-qc]].
  * [64] T. Akutsu et al., Phys. Rev. Lett. 101 (2008), 101101 [arXiv:0803.4094 [gr-qc]].
  * [65] V. Domcke and C. Garcia-Cely, Phys. Rev. Lett. 126 (2021) no.2, 021104 [arXiv:2006.01161 [astro-ph.CO]].
  * [66] A. Ito, T. Ikeda, K. Miuchi and J. Soda, Eur. Phys. J. C 80 (2020) no.3, 179 [arXiv:1903.04843 [gr-qc]].
  * [67] A. Ito and J. Soda, Eur. Phys. J. C 80 (2020) no.6, 545 [arXiv:2004.04646 [gr-qc]].
  * [68] A. Ejlli, D. Ejlli, A. M. Cruise, G. Pisano and H. Grote, Eur. Phys. J. C 79 (2019) no.12, 1032 [arXiv:1908.00232 [gr-qc]].
  * [69] V. Domcke, C. Garcia-Cely and N. L. Rodd, Phys. Rev. Lett. 129 (2022) no.4, 041101 [arXiv:2202.00695 [hep-ph]].
  * [70] A. Ringwald, J. Schütte-Engel and C. Tamarit, JCAP 03 (2021), 054 [arXiv:2011.04731 [hep-ph]].
  * [71] A. Berlin, D. Blas, R. Tito D’Agnolo, S. A. R. Ellis, R. Harnik, Y. Kahn and J. Schütte-Engel, Phys. Rev. D 105 (2022) no.11, 116011 [arXiv:2112.11465 [hep-ph]].
  * [72] K. Schmitz, “New Sensitivity Curves for Gravitational-Wave Experiments” [Data set], Zenodo, [http://doi.org/10.5281/zenodo.3689582].
  * [73] L. Pagano, L. Salvati and A. Melchiorri, Phys. Lett. B 760 (2016), 823-825 [arXiv:1508.02393 [astro-ph.CO]].
  * [74] T. J. Clarke, E. J. Copeland and A. Moss, JCAP 10 (2020), 002 [arXiv:2004.11396 [astro-ph.CO]].
  * [75] J. Ghiglieri, G. Jackson, M. Laine and Y. Zhu, JHEP 07 (2020), 092 [arXiv:2004.11392 [hep-ph]].
  * [76] J. Ghiglieri and M. Laine, JCAP 07 (2015), 022 [arXiv:1504.02569 [hep-ph]].
  * [77] N. Aggarwal et al., Living Rev. Rel. 24 (2021) no.1, 4 [arXiv:2011.12414 [gr-qc]].