# Cosmological phase transitions, gravitational waves and self-interacting
dark matter in the singlet extension of MSSM

Wenyu Wang<EMAIL_ADDRESS>Faculty of Science, Beijing University of
Technology, Beijing, P. R. China  Ke-Pan Xie<EMAIL_ADDRESS>Department of
Physics and Astronomy, University of Nebraska, Lincoln, NE 68588, USA  Wu-Long
Xu<EMAIL_ADDRESS>Faculty of Science, Beijing University of
Technology, Beijing, P. R. China  Jin Min Yang<EMAIL_ADDRESS>CAS Key
Laboratory of Theoretical Physics, Institute of Theoretical Physics, Chinese
Academy of Sciences, Beijing 100190, P. R. China  School of Physical Sciences,
University of Chinese Academy of Sciences, Beijing 100049, P. R. China

###### Abstract

In the minimal supersymmetric standard model (MSSM) extended by a singlet
superfield, when the coupling between the singlet sector and the MSSM sector
is tiny, the singlet sector can be a quasi dark sector with supersymmetry
(SUSY). We investigate the cosmological phenomena in this scenario and obtain
the following observations: (i) In the parameter space solving the small
cosmological scale anomalies via self-interacting singlino dark matter (SIDM),
a first-order phase transition (FOPT) can readily happen but requires rather
light dark matter below MeV; (ii) The corresponding parameter space indicated
by FOPT and SIDM can be partially covered by detecting the phase-transition
gravitational waves (GWs) at the near-future projects, such as LISA, TianQin
and Taiji. Therefore, the recently developed GW astronomy could be a novel
probe to such a SUSY scenario.

###### Contents

  1. I Introduction
  2. II Self-interacting singlino dark matter
  3. III Cosmic FOPT and gravitational waves
    1. III.1 First-order phase transition
    2. III.2 Gravitational waves
  4. IV Numerical results
  5. V Conclusion
  6. A Feynman rules in the singlet SUSY model

##  I Introduction

Dark matter (DM) and gravitational wave (GW) are hot topics in the frontier of
both theoretical and experimental physics studies. The DM is indicated from
lots of astrophysical observations, which is non-luminous and distributes all
over the Universe with a density approximately five times of the ordinary
matter Young (2017). Usually the DM is considered as a fundamental particle
which is weakly interacting with the Standard Model (SM) particles to give the
correct relic density, i.e., the so-called WIMP (Weakly Interacting Massive
Particle) miracle Dodelson (2003). As a mystery of the fundamental structure
of matters, the DM issue is on the list of the physics beyond the SM and needs
to be deeply studied. On the other hand, the GW, which was a simple prediction
of classical general relativity long time ago, has been measured very recently
due to the great improvements in the measurements of space and time.
Superficially these two topics are unrelated, which, however, become
correlated in case of a first-order cosmological phase transition (FOPT) Linde
(1983).

In the particle physics side, the spontaneous symmetry breaking of the Higgs
field, which is a fundamental field of phase transition, gives the masses of
particles. However, there is no DM candidate in the SM. And the SM suffers
from both aesthetical judgements and experimental observations. The hierarchy
problem, the strong CP problem and the Grand Unification Theory (GUT) require
new content of particles or new symmetries; in addition to DM, the neutrino
oscillation implies finite mass of neutrinos. Although finite neutrino mass
terms can be introduced conventionally, in the form of Higgs terms without a
significant modification of the theory, this makes the SM to appear rather
contrived. Thus the physics beyond the SM with an upgraded phase transition
usually has to solve the hierarchy problem, realize GUT, and give candidates
of the DM. In the GW side, if the potential of a scalar field has two
different energy minima, the FOPT can happen via quantum tunneling between the
two vacua separated by a barrier Coleman (1977); Callan and Coleman (1977). In
the early Universe, as temperature was dropping, this quantum tunneling may
also occur in the Higgs field. During the FOPT, the collision of vacuum
bubbles, sound waves in the plasma and the magneto-hydrodynamics turbulence
could radiate GWs that may be detected by the future space-based LISA Amaro-
Seoane et al. (2017), TianQin Luo et al. (2016) or Taiji Ruan et al. (2020)
experiments. Thus, for an unified description of the underlying fundamental
theory, the results of phase transition should be checked from both sides
since phase transition in fact gives a platform for both fundamental particles
and gravitational phenomena, or the detection of phase transition should give
constraints on fundamental particle models.

The Monte Carlo simulations demonstrate that the electroweak phase transition
in the SM is a crossover Kajantie et al. (1996); Fodor et al. (1995); Kajantie
et al. (1993). It shows that the extension of the SM is needed for a FOPT
Espinosa et al. (2012); Athron et al. (2019a); Apreda et al. (2002); Hasegawa
et al. (2019); Bian et al. (2019); Xie et al. (2020); Su et al. (2021); Ghosh
et al. (2021); Guo et al. (2021a); Gonçalves et al. (2021). Among various new
physics models, low energy supersymmetry (SUSY) is most widely studied since
it can elegantly solve both aesthetical and experimental problems of the SM
Wess and Zumino (1974); Sohnius (1985); Sakai (1981); Dimopoulos and Georgi
(1981) (for review reviews, see, e.g., Refs. Baer et al. (2020); Wang et al.
(2022)). However, the FOPT can be hardly realized in the minimal
supersymmetric standard model (MSSM) because the needed large barrier of the
SM Higgs potential requires light stops (superpartners of top quark) not
allowed by current search results at colliders. This problem can be easily
solved by adding a singlet superfield which can have tree-level barrier
freely. The FOPT of such singlet extension of MSSM, which is called NMSSM
(GNMSSM) with (without) $Z_{3}$ symmetry Ellwanger et al. (2010); Cao et al.
(2012), has been studied intensively in the literature Huang et al. (2015); Bi
et al. (2015); Bian et al. (2018); Athron et al. (2019a); Chatterjee et al.
(2022). In these studies, the phase transition is always considered at the
electroweak scale and the singlet couples with Higgs doublets to play the role
of spontaneous symmetry breaking. Generally speaking, the coupling between the
singlet and Higgs doublets, which is the source of electroweak phase
transition, is not the essential part in theoretical models. The phase
transition of singlet itself can have more advantages compared with the mixing
with the doublets. The most promising point is that the singlet superfield can
play a very subtle role in the scenario of light DM and self-interacting DM
(SIDM) favored by the observation of small structures of the Universe. The
singlino can be a light DM candidate, annihilating to the singlet scalars to
give correct relic abundance, interchanging with singlet to give appropriate
self-interaction of singlino DM and solve the small structure problem of the
WIMP. The most important thing is that the singlino DM with tiny coupling to
the doublets can easily escape the current stringent direct detection limits
of the DM. Thus, this scenario has charming features in DM physics, which can
realize SIDM scenario and solve the small cosmological scale anomalies Wang et
al. (2014); Zhu and Abdughani (2021).

With these advantages in fundamental particle physics, the scenario of a
singlet superfield with tiny couplings to the MSSM sector should be checked
from the cosmological phase transition and the detection of GWs. In details,
we should check: whether the FOPT can be maintained, whether the surviving
parameter space can be probed by the detection of the GWs, or whether the GW
detection can cover the parameter space of SIDM. In this work we examine these
issues. This work is organized as follows. In Section II we show the content
of the SUSY singlet model and show how the small structure problem is solved.
In Section III, the FOPT and GWs are discussed. In Section IV the numerical
results are shown. The conclusion is given in Section V.

##  II Self-interacting singlino dark matter

As addressed in the preceding section, the non-luminous DM is a necessary
component in the Universe. The standard cosmological model $\Lambda$CDM, which
assumes cold and collisionless DM, has achieved remarkable success in
describing the Universe at large scales greater than $\cal O$ (Mpc) today.
However, this model encounters lots of crises on the small-scale Universe at
which the structure formation becomes strongly non-linear and the $N$-body
simulation is the standard tool in this regime. The numerical results showed
that the weakly interacting DM cannot form the rich structures observed in the
Universe. Such discrepancies are called the small structure problems,
including the core-cusp problem Navarro et al. (1997); Moore et al. (1999),
the diversity problem Kuzio de Naray et al. (2010); Bullock et al. (2001);
Oman et al. (2015), the missing satellites problem Klypin et al. (1999);
Kauffmann et al. (1993); Zavala et al. (2009) and the too-big-to-fail problem
Boylan-Kolchin et al. (2012); Tollerud et al. (2014); Garrison-Kimmel et al.
(2014).

Though in the literature there are some debates about the robustness of these
problems, they may indeed indicate some shortcomings of CDM in describing the
small structures of the Universe. One way to solve these problems is to
propose warm DM, which, however, is not consistent with some other
observations Iršič et al. (2017); Viel et al. (2013); Menci et al. (2016) and
gives too small core sizes to solve the core-cusp problem Maccio et al.
(2012).

Another promising alternative approach to solve the small structure problems
is to introduce self-interacting dark matter (SIDM) Spergel and Steinhardt
(2000). For SIDM the elastic cross section is velocity-dependent and thus the
light mediator of SIDM can give sufficiently large cross sections to satisfy
the small structure requirements, while the annihilation can still maintain
weak to give a correct freezing-out relic density in the early Universe.

A detailed study (shown in the proceeding section) of SIDM shows that the
mediators or the force carriers prefer at about several MeV or even lighter.
Such a hierarchy between the mediator and electroweak scale makes it very
difficult to connect with the electroweak scale physics in the scheme of the
new physics beyond SM. Further more, as shown in Ref. Wang et al. (2014); Elor
et al. (2021), if the light mediator couples to the nucleons, the correlation
between the DM annihilation rate and DM-nucleon spin independent (SI) cross
section has almost excluded the scenario. Thus from the point view of SIDM,
the introduction of a singlet is not necessarily connected to the electroweak
sector.

With a tiny connection to electroweak sector, a singlet in the SUSY framework
can be a perfect building block for SIDM. The singlino composes the DM,
annihilating to the singlet bosons to give the correct relic abundance, and
interacting by exchanging singlet bosons to give large scattering elastic
cross sections in case of small velocities for the small structure
cosmological scales. In the following, we will show the details of such a
singlet SUSY model.

The superpotential of a singlet in the most general form is

$W=\eta\hat{S}+\frac{1}{2}\mu_{s}\hat{S}^{2}+\frac{1}{3}\kappa\hat{S}^{3},$
(1)

where $\hat{S}$ is the singlet superfield, $\eta$ and $\mu_{s}$ are parameters
with mass dimension $m^{2}$ and $m$ respectively, and $\kappa$ is a
dimensionless coupling strength of the field. The fermion part of $\hat{S}$
gives a Majorana fermion singlino, while the scalar part $S$ gives one scalar
and one pseudo-scalar. This is actually the modified Wess-Zumino model, with
simply degenerate mass spectrum. A physically allowed model needs soft
breaking terms

$\displaystyle V_{\rm soft}=m_{s}^{2}|S|^{2}+(C_{\eta}\eta
S+\frac{1}{2}B_{s}\mu_{s}S^{2}+\frac{1}{3}\kappa A_{k}S^{3}+{\rm h.c.}),$ (2)

and then the potential becomes

$V=V_{F}+V_{\rm soft}=|\kappa S^{2}|^{2}+m_{s}^{2}|S|^{2}+\left(C_{\eta}\eta
S+\frac{1}{2}B_{s}\mu_{s}S^{2}+\frac{1}{3}\kappa
A_{k}S^{3}+\kappa\mu_{s}S^{2}S^{*}+{\rm h.c.}\right).$ (3)

Note that some of the parameters are absorbed by the redefinition of the soft
parameters, giving a more economic potential. We can see that this potential
can easily spontaneously break to get the vacuum expectation value (VEV)
$v_{s}$. Thus we denote

$\displaystyle S=v_{s}+\frac{1}{\sqrt{2}}(H+iA).$ (4)

The tadpole vanishes in the condition of

$\displaystyle m_{s}^{2}=-2v_{s}^{2}\kappa^{2}-v_{s}\kappa
A_{\kappa}-\frac{\eta C_{\eta}}{v_{s}}-(3v_{s}\kappa+B_{s})\mu_{s}.$ (5)

After a short calculation, we can get the masses of the singlino $\chi$, CP-
even scalar $H$ and CP-odd pseudo-scalar $A$:

$\begin{split}m_{\chi}&=~{}\mu_{s}+2\kappa v_{s},\\\
m_{H}^{2}&=~{}6v_{s}^{2}\kappa^{2}+2v_{s}\kappa
A_{\kappa}+m_{s}^{2}+6v_{s}\kappa\mu_{s}+B_{s}\mu_{s},\\\
m_{A}^{2}&=~{}2v_{s}^{2}\kappa^{2}-2v_{s}\kappa
A_{\kappa}+m_{s}^{2}+2v_{s}\kappa\mu_{s}-B_{s}\mu_{s}.\end{split}$ (6)

The corresponding Feynman rules are shown in the appendix. Here we just
mention that singlino $\chi$ can have sizable self-interaction by exchanging a
light scalar, and the annihilation to the scalars can give a proper relic
abundance of singlino DM. The corresponding Feynman diagrams of the
annihilation are shown in Fig. 1.

Figure 1:  Feynman diagrams of singlino annihilation to the salars, with
$\phi=H$, $A$.

Note that, in this singlet SUSY model, the DM annihilation has more channels
than in the simple SIDM model with only one DM and one mediator, in which the
DM can only annihilate to the scalar via $t$\- and $u$-channels (the first two
diagrams in Fig. 1). In our model, the scalars have self-couplings, thus the
annihilation can have $s$-channel (the third diagram in Fig. 1). This will
make the parameter space different from the simple SIDM model. Nevertheless,
the self-interactions of scalars are in fact derived from the scalar potential
Eq. (3), which is the starting point of the FOPT for the detection of the GWs,
and thus detailed studies are needed to check the connection between
annihilation, self-interaction and GW detection.

After the derivation of the basic model, an important issue is the parameter
definition. Though we can use the soft parameters such as $\mu$ and $C_{\eta}$
to study the model, the physical implications of these parameters are obscure.
Thus we chose $m_{\chi}$, $m_{H}$, $m_{A}$, $\alpha_{\chi}$, $A_{\kappa}$,
$v_{s}$ as the parameters of the model, which are masses of the particles, the
coupling strength, the trilinear term and the scalar VEV, respectively. Note
that we chose $\alpha_{\chi}=4\kappa^{2}/\pi$ for a more convenient definition
of the coupling strength. Then, the scalar potential Eq. (3) can be obtained
from Eq. (5) and we have the following relations

$\begin{split}\kappa=&~{}\frac{\sqrt{\pi\alpha_{\chi}}}{2},~{}~{}~{}~{}\mu_{s}=m_{\chi}-2\kappa
v_{s},\\\ \eta
C_{\eta}=&~{}\frac{1}{2}v_{s}(-m_{H}^{2}+(A_{\kappa}+3m_{\chi})\sqrt{\pi\alpha_{\chi}}v_{s}-\pi\alpha_{\chi}v_{s}^{2}),\\\
B_{s}=&~{}\frac{m_{H}^{2}-m_{A}^{2}+v_{s}(-2m_{\chi}+\sqrt{\pi\alpha_{\chi}}(-4A_{\kappa}-3m_{\chi}+2v_{s})+\pi\alpha_{\chi}v_{s})}{4(m_{\chi}-\sqrt{\pi\alpha_{\chi}}v_{s})}.\end{split}$
(7)

These parameters are basic inputs for the scalar potential. Before going to
the detailed study of FOPT, we should note that our scalar potential is not a
simple construction from the scalar field with specific gauge symmetry,
instead it is just a general construction of a singlet superfield, or namely,
the soft breaking Wess-Zumino SUSY model. Thus the model has particular
meaning on the study of SUSY, a mystery symmetry of the space and time of our
Universe.

##  III Cosmic FOPT and gravitational waves

Cosmological FOPT is the decay of the Universe between the false and true
vacua, which are separated by a potential barrier Linde (1983). In our case,
the FOPT scalar field is the CP-even background scalar field $h$, whose
potential receives the quantum and thermal corrections from the CP-even scalar
itself, the CP-odd scalar and the singlino loops. A detailed study of FOPT and
detection of the GWs are discussed in the following.

###  III.1 First-order phase transition

We work under the CP conservation assumption such that only the real part of
the complex field $S=(H+iA)/\sqrt{2}$ gains a VEV, and hence we only consider
the CP-even background field $h$. At zero temperature, the effective potential
at the tree level is

$\displaystyle V_{0}(h)=\kappa^{2}h^{4}+m_{s}^{2}h^{2}+2\left(C_{\eta}\eta
h+\frac{1}{2}B_{s}\mu_{s}h^{2}+\frac{1}{3}\kappa
A_{k}h^{3}+\kappa\mu_{s}h^{3}\right).$ (8)

At one-loop level, the potential receives the Coleman-Weinberg corrections

$\displaystyle V_{1}^{T=0}(h)$ $\displaystyle=$
$\displaystyle\frac{1}{64\pi^{2}}m_{a}(h)^{4}\left[\log\frac{m_{a}^{2}(h)}{\mu^{2}}-\frac{3}{2}\right]+\frac{1}{64\pi^{2}}m_{h}^{4}(h)\left[\log\frac{m_{h}^{2}(h)}{\mu^{2}}-\frac{3}{2}\right],$
(9)
$\displaystyle\,-\frac{1}{32\pi^{2}}m_{\chi}^{4}(h)\left[\log\frac{m_{\chi}^{2}(h)}{\mu^{2}}-\frac{3}{2}\right],$

where we use the $\overline{\rm MS}$ subtraction scheme and take the
renormalization scale as $\mu=246~{}{\rm GeV}$. Note that there is a SUSY
limit in the potential. $V_{1}^{T=0}(h)$ will approach zero when the masses of
$\chi$, $H$, $A$ approach to degenerate, which is an interesting point in SUSY
theory. The field-dependent masses are obtained by replacing $v_{s}$ with $h$
in Eq. (6). To remain the tree level VEV and mass relations, we also add the
following counter terms

$\displaystyle\delta V_{1}^{T=0}(h)=\delta m_{s}^{2}h^{2}+2\left(\delta
C_{\eta}\eta h+\frac{1}{2}B_{s}\delta\mu_{s}h^{2}+\frac{1}{3}\delta\kappa
A_{k}h^{3}\right),$ (10)

and choose the renormalization condition as

$\displaystyle\frac{d}{dh}\left(V_{1}^{T=0}+\delta
V_{1}^{T=0}\right)\Big{|}_{h=v_{s}}=0,\quad\frac{d^{2}}{dh^{2}}\left(V_{1}^{T=0}+\delta
V_{1}^{T=0}\right)\Big{|}_{h=v_{s}}=0.$ (11)

At finite temperature, the potential receives one-loop thermal integral
corrections

$\displaystyle V_{1}^{T\neq
0}(h,T)=\frac{T^{4}}{2\pi^{2}}J_{B}\left(\frac{m_{h}^{2}(h)}{T^{2}}\right)+\frac{T^{4}}{2\pi^{2}}J_{B}\left(\frac{m_{a}^{2}(h)}{T^{2}}\right)-4\frac{T^{4}}{2\pi^{2}}J_{F}\left(\frac{m_{\chi}^{2}(h)}{T^{2}}\right),$
(12)

where

$J_{B/F}(y)=\int_{0}^{\infty}dxx^{2}\log\left(1\mp
e^{-\sqrt{x^{2}+y}}\right).$ (13)

We also include the daisy resummation contribution

$V_{\rm
daisy}(h)=-\frac{T}{12\pi}\left[\left(m_{h}^{2}(h)+\Pi_{h}T^{2}\right)^{3/2}+\left(m_{a}^{2}(h)+\Pi_{a}T^{2}\right)^{3/2}-m_{h}^{3}(h)-m_{a}^{3}(h)\right],$
(14)

where $\Pi_{h}=\kappa^{2}T^{2}$, $\Pi_{a}=\kappa^{2}T^{2}$.

The total effective potential up to one-loop can be written as

$V_{\rm eff}(h,T)=V_{0}(h)+V_{1}^{T=0}(h)+\delta V_{1}^{T=0}(h)+V_{1}^{T\neq
0}(h,T)+V_{\rm daisy}(h,T).$ (15)

The temperature that $V_{\rm eff}(h,T)$ has two degenerate vacua is defined as
the critical temperature $T_{c}$. Below $T_{c}$, the true vacuum is
energetically preferred and the Universe will decay from the false vacuum to
the true vacuum. The decay rate per unit volume is given by Linde (1983)

$\Gamma\sim T^{4}\left(\frac{S_{3}(T)}{2\pi T}\right)^{3/2}e^{-S_{3}(T)/T},$
(16)

where the Euclid action $S_{3}$ is

$S_{3}=\int_{0}^{\infty}4\pi
r^{2}dr\left[\frac{1}{2}\left(\frac{d\hat{h}}{dr}\right)^{2}+V_{\rm
eff}(\hat{h},T)\right],$ (17)

with $\hat{h}(r)$ being the $O(3)$-symmetric bounce solution determined by

$\frac{d^{2}\hat{h}}{dr^{2}}+\frac{2}{r}\frac{d\hat{h}}{dr}=\frac{\partial}{\partial\hat{h}}V_{\rm
eff}(\hat{h},T),\quad\lim_{r\to\infty}\hat{h}=0,\quad\frac{d\hat{h}}{dr}\Big{|}_{r=0}=0.$
(18)

The probability that a bubble nucleates inside a Hubble volume is

$N(T)=\int_{T}^{T_{c}}\frac{dT^{\prime}}{T^{\prime}}\frac{\Gamma(T^{\prime})}{H^{4}(T^{\prime})},$
(19)

where $H(T)$ is the Hubble constant. The temperature satisfying $N(T_{n})\sim
1$ is called the nucleation temperature $T_{n}$. This is resolved by

$\frac{S_{3}}{T_{n}}\simeq\ln\left[\frac{1}{4}\left(\frac{90}{8\pi^{3}g_{*}}\right)^{2}\right]+4\ln\left[\frac{M_{\rm
Pl}}{T_{n}}\right],$ (20)

where $M_{\rm Pl}=1.22\times 10^{19}$ GeV is the Planck scale, and $g_{*}$ is
the number of relativistic degrees of freedom. Here, we can take $T_{n}$ of
right hand side of Eq. (20) as $T_{c}$ approximately because of the
logarithmic function. In this article, we use the public package
cosmoTransition Wainwright (2012)111See Refs. Athron et al. (2020, 2019b) for
other packages for calculating the FOPT. to calculate the bounce solution and
determine $T_{n}$.

###  III.2 Gravitational waves

Next, we show the details of the radiation of the GWs. The volume fraction of
the false vacuum of the Universe is given by Guth and Weinberg (1981)

$p(T)=e^{-I(T)},\quad
I(T)=\frac{4\pi}{3}\int_{T}^{T_{c}}dT^{\prime}\frac{\Gamma(T^{\prime})}{T^{\prime
4}H(T^{\prime})}\left[\int_{T}^{T^{\prime}}d\tilde{T}\frac{v_{b}}{H(\tilde{T})}\right]^{3},$
(21)

which decreases to zero as the FOPT completes. Here $v_{b}$ is the bubble
expansion velocity relative to the plasma at finite distance. The temperature
that the true vacuum bubbles form an infinite connected cluster is called the
percolation temperature $T_{p}$. Numerical simulation shows that
$p(T_{p})=0.71$ Rintoul and Torquato (1997). $T_{p}$ is the characteristic
temperature for calculating the FOPT GWs Megevand and Ramirez (2017);
Kobakhidze et al. (2017); Ellis et al. (2019, 2020); Wang et al. (2020).

A FOPT generates stochastic GWs via bubble collision, sound waves and
turbulence of the magneto-hydrodynamics (MHD) in the plasma. Because of the
friction in the plasma-wall system, the expansion of bubble wall will
accelerate for only a short time and then reach the terminal velocity $v_{b}$.
Therefore, only a small fraction of FOPT energy is stored in the wall and the
bubble collision contribution is negligible. Instead, most energy is in the
fluid shells surrounding the wall, making the sound waves the dominant
contribution Ellis et al. (2019). The MHD turbulence contributes a sub-leading
source to the GWs. Defining the GW spectrum today as

$\Omega_{\rm GW}(f)=\frac{1}{\rho_{c}}\frac{\rho_{\rm GW}}{d\ln f},$ (22)

where $f$ is the frequency, $\rho_{\rm GW}$ is the GW energy density and
$\rho_{c}$ is the critical energy density of the present Universe, we
have222Here $h$ is $H_{0}/(100~{}{\rm km/s/Mpc})$, not to be confused with the
background field $h$.

$\Omega_{\rm GW}h^{2}\simeq\Omega_{\rm sw}h^{2}+\Omega_{\rm turb}h^{2},$ (23)

where the sound wave contribution $\Omega_{\rm sw}$ and turbulence
contribution $\Omega_{\rm turb}$ can be expressed as the numerical functions
of two parameters from the FOPT profile Grojean and Servant (2007); Caprini et
al. (2016, 2020): i) $\alpha$, defined as the ratio of the latent heat in the
FOPT to the radiation energy density of the Universe; and ii) $\beta/H_{p}$,
the inverse ratio of duration of FOPT and the Hubble time. Quantitatively,
their definitions are

$\alpha=\frac{1}{\rho_{R}(T_{p})}\left(T\frac{\partial\Delta V_{\rm
eff}}{\partial T}-\Delta V_{\rm
eff}\right)\Big{|}_{T_{p}};\quad\frac{\beta}{H_{p}}=T_{p}\frac{d(S_{3}/T)}{dT}\Big{|}_{T_{p}},$
(24)

where $\Delta V_{\rm eff}$ is the free energy difference between true and
false vacua and $\rho_{R}(T)=g_{*}\pi^{2}T^{4}/30$ is the radiation energy
density (with $g_{*}$ being the number of relativistic degrees of freedom). We
adopt $T_{p}\approx T_{n}$ for the GW calculation, as the supercooling effect
is not strong in our scenario, and the characteristic temperatures are very
close. Because the coupling between dark sector and SM is extremely weak,
after reheating stage, those two parts do not have any chance exchanging
energy and injecting entropy until the dark sector phase transition. Thus, the
total number of relativistic degrees of freedom is $g_{*}=g^{\rm
SM}_{*}+g^{d}_{*}(\frac{T}{T_{\rm SM}})$ when the dark sector temperature is
$T$ and the SM temperature is $T_{\rm SM}$. Here we take $g^{d}_{*}=3.75$
counting the degrees of freedom from $\chi$, $H$ and $A$. This latent heat
changes the dark sector temperature form $T_{p}$ to $T_{f}$. The result is
written as

$\frac{\pi^{2}g_{*}^{d}}{30}T_{f}^{4}=\frac{\pi^{2}g_{*}^{d}}{30}T_{p}^{4}+\rho_{\rm
dvac}.$ (25)

The parameter $\alpha$ and $\frac{\beta}{H_{p}}$ can be used to characterize
the strength of the FOPT.

Numerically, the sound wave contribution reads

$\displaystyle h^{2}\Omega_{\rm sw}(f)=2.65\times 10^{-6}(H_{p}\tau_{\rm
sw})\left(\frac{H_{p}}{\beta}\right)\left(\frac{\kappa_{v}\alpha}{1+\alpha}\right)^{2}\left(\frac{100}{g_{*}}\right)^{\frac{1}{3}}v_{b}S_{\rm
sw}(f),$ (26) $\displaystyle S_{\rm sw}(f)=\left(\frac{f}{f_{\rm
sw}}\right)^{3}\left(\frac{7}{4+3(\frac{f}{f_{\rm
sw}})^{2}}\right)^{\frac{7}{2}},$ (27) $\displaystyle f_{\rm sw}=1.9\times
10^{-2}~{}{\rm
mHz}\frac{1}{v_{b}}\left(\frac{\beta}{H_{p}}\right)\left(\frac{T_{p}}{100~{}{\rm
GeV}}\right)\left(\frac{g_{*}}{100}\right)^{\frac{1}{6}},$ (28)

where $\tau_{\rm sw}$ is the duration of the sound wave source Ellis et al.
(2019); Guo et al. (2021b), and $\kappa_{v}$ is the fraction of vacuum energy
distributed in the bubble kinetic energy and can be found in Ref. Espinosa et
al. (2010). The turbulence contribution is

$\displaystyle h^{2}\Omega_{{\rm turb}}(f)=3.35\times
10^{-4}\left(\frac{H_{p}}{\beta}\right)\left(\frac{\kappa_{\rm
turb}\alpha}{1+\alpha}\right)^{2}\left(\frac{100}{g_{*}}\right)^{\frac{1}{3}}v_{b}S_{\rm
sturb}(f),$ (29) $\displaystyle S_{\rm turb}(f)=\frac{(\frac{f}{f_{\rm
turb}})^{3}}{[1+(\frac{f}{f_{\rm turb}})]^{\frac{11}{3}}\left(1+\frac{8\pi
f}{h_{p}}\right)},$ (30) $\displaystyle f_{\rm sw}=2.7\times 10^{-2}~{}{\rm
mHz}\frac{1}{v_{b}}\left(\frac{\beta}{H_{p}}\right)\left(\frac{T_{p}}{100~{}{\rm
GeV}}\right)\left(\frac{g_{*}}{100}\right)^{\frac{1}{6}},$ (31)
$\displaystyle\kappa_{\rm turb}=\epsilon\kappa_{v},$ (32)

where $\kappa_{\rm turb}$ denotes the fraction of fraction of vacuum energy
distributed in the turbulence and $H_{p}$ is the Hubble parameter at $T_{p}$.

##  IV Numerical results

Before performing specific numerical calculation, we briefly review the
calculation of DM relic density and SIDM cross section. The Feynman diagrams
of the annihilation are shown in Fig. 1 and the annihilation cross section is
Drees and Nojiri (1993); Jungman et al. (1996)

$\sigma v=\frac{1}{4}\frac{\bar{\beta}_{f}}{8\pi
sS}\left[|A(^{1}S_{0})|^{2}+\frac{1}{3}\left(|A(^{3}P_{0})|^{2}+|A(^{3}P_{1})|^{2}\right)+|A(^{3}P_{2})|^{2}\right],$
(33)

where $S$ is the symmetry factor, amplitudes
$A(^{1}S_{0}),~{}A(^{3}P_{0}),~{}A(^{3}P_{1}),~{}A(^{3}P_{2})$ are the
contributions from different spin states of DM which are shown in the
appendix, and $\bar{\beta}_{f}$ is given by

$\displaystyle\bar{\beta}_{f}=\sqrt{1-2(m_{X}^{2}+m_{Y}^{2})/s+(m_{X}^{2}+m_{Y}^{2})^{2}/s^{2}}~{},$
(34)

with $X$, $Y$ being the final states.

For the SIDM, we use the cross sections defined in Ref. Krstić and Schultz
(1999),333We assume a instant freeze-out such that the annihilation cross
section can be adopted as the WIMP one Vanderheyden (2021). the transfer cross
section $\sigma_{T}$ and the viscosity (or conductivity) cross section
$\sigma_{V}$: Tulin et al. (2013a, b); Ko and Tang (2014)

$\sigma_{T}=\int
d\Omega\,(1-\cos\theta)\,\frac{d\sigma}{d\Omega}\,,\qquad\sigma_{V}=\int
d\Omega\,\sin^{2}\theta\,\frac{d\sigma}{d\Omega}\,.$ (35)

$\sigma_{T}$ is for the estimation of the Dirac DM, while $\sigma_{T}$ is for
the Majorana DM. Since the singlino $\chi$ is a Majorana fermion, the
viscosity cross section is defined with two variables:

$\displaystyle\frac{d\sigma_{VS}}{d\Omega}$ $\displaystyle=$
$\displaystyle\left|f(\theta)+f(\pi-\theta)\right|^{2}=\frac{1}{k^{2}}\left|\sum_{\ell(\mbox{\tiny\rm
EVEN\
number})}^{\infty}(2\ell+1)(\exp(2i\delta_{l})-1)P_{\ell}(\cos\theta)\right|^{2},$
(36) $\displaystyle\frac{d\sigma_{VA}}{d\Omega}$ $\displaystyle=$
$\displaystyle\left|f(\theta)-f(\pi-\theta)\right|^{2}=\frac{1}{k^{2}}\left|\sum_{\ell(\mbox{\tiny\rm
ODD\
number})}^{\infty}(2\ell+1)(\exp(2i\delta_{l})-1)P_{\ell}(\cos\theta)\right|^{2}.$
(37)

Using the orthogonality relation for the Legendre polynomials, we can get

$\displaystyle\frac{\sigma_{VS}k^{2}}{4\pi}$ $\displaystyle=$
$\displaystyle\sum_{\ell(\mbox{\tiny EVEN\
number})}^{\infty}4\sin^{2}(\delta_{\ell+2}-\delta_{\ell})(\ell+1)(\ell+2)/(2\ell+3),$
(38) $\displaystyle\frac{\sigma_{VA}k^{2}}{4\pi}$ $\displaystyle=$
$\displaystyle\sum_{\ell(\mbox{\tiny ODD\
number})}^{\infty}4\sin^{2}(\delta_{\ell+2}-\delta_{\ell})(\ell+1)(\ell+2)/(2\ell+3).$
(39)

The phase shift $\delta_{\ell}$ must be computed by solving the Schrödinger
equation

$\frac{1}{r^{2}}\frac{d}{dr}\Big{(}r^{2}\frac{dR_{\ell}}{dr}\Big{)}+\Big{(}k^{2}-\frac{\ell(\ell+1)}{r^{2}}-2m_{r}V(r)\Big{)}R_{\ell}=0,$
(40)

directly with partial wave expansion method. The total wave function of the
spin-1/2 fermionic DM must be antisymmetric with respect to the exchange of
two identical particles. Then the spatial wave function should be symmetric
when the total spin is 0 (singlet) while the spatial wave function should be
antisymmetric when the total spin is 1 (triplet). Then we consider that the DM
scatters with random orientations, thus the triplet is three times as likely
as the singlet and the average cross section is

$\displaystyle\sigma_{V}=\frac{1}{4}\sigma_{VS}+\frac{3}{4}\sigma_{VA}\;.$
(41)

We need to do a specific numerical calculation to check the existence of
parameter space which can satisfy the constraint of relic density, realize the
SIDM scenario and give detectable GWs.

After all preparations are completed, the parameters $m_{\chi}$, $m_{H}$,
$m_{A}$, $\kappa$, $A_{\kappa}$ are chosen in ranges of

$\begin{split}&10^{-6}~{}{\rm GeV}<m_{\chi}<10^{4}~{}{\rm GeV},\quad
10^{-6}~{}{\rm GeV}<m_{H}<10^{4}~{}{\rm GeV},\\\ &10^{-6}~{}{\rm
GeV}<m_{A}<10^{4}~{}{\rm GeV},\quad 10^{-6}<\kappa<1,\quad-10^{5}~{}{\rm
GeV}<A_{\kappa}<0.\end{split}$ (42)

According to the experimental constraints, the favourable relic density is
$(0.107,0.131)$. For the SIDM cross section $\sigma/m_{\chi}$, Ref. Tulin and
Yu (2018) shows the different value in the different small scale structure
problems. From a pure theoretical point of view, it is sufficient to choose
$\sigma/m_{\chi}$ in the range of $(0.1,10~{}\rm{cm^{2}/g})$ and the
characteristic velocity is $200~{}\rm{km/s}$. The parameter space allowed by
DM relic density and the small scale structure is shown in Fig. 2 as the blue
points. We can see that the allowed mass of SIDM for the small scale structure
can be from keV to tens of GeV, while the mass of dark scalar is below a few
MeV.

Figure 2: The blue points are allowed by DM relic density and the small scale
structure via SIDM; the red region can give a FOPT.

For the calculation of FOPT, there is one more free parameter, namely the VEV
$v_{s}$. By varying $v_{s}\in[10^{-6}~{}{\rm GeV},10^{4}~{}{\rm GeV}]$, we
obtain the parameter space giving a FOPT, shown as the red region in Fig. 2.
The intersection region between blue points and the red region in Fig. 2 is
the surviving parameter space which can give a correct relic density and
realize FOPT and SIDM. It should be noted that in recent studies the filtered
bubble can reflect the DM particles because of the energy defect between the
value of false vacuum and real vacuum of the scalar in which Fermi-ball or
Q-ball can be formed Krylov et al. (2013); Hong et al. (2020); Kawana and Xie
(2022). But in this SUSY model, the DM mass gap between the two vacua is
small. Thus, the DM particle can pass the bubble wall safely. Then this FOPT
can not affect the DM generation and the relic density of DM can be obtained
via the standard freeze-out mechanism. From Fig. 2, we can see that in the
survived region the SIDM is lighter than 1 MeV and the dark scalar is lighter
than 100 KeV.

Figure 3: The SNRs at different near-future detectors (left: LISA, middle:
TianQin, right: Taiji), with the blue points for SNR $<$ 10, the green points
for $10<$ SNR $<$ 50 and the red points for SNR $>$ 50\. The data-taking
duration is adopted as four years.

According to Refs. Marfatia and Tseng (2021); Di Bari et al. (2021), the
effective number of extra neutrino species contributed by the SUSY dark sector
$\Delta N_{\rm eff}\propto\frac{T_{f}}{T_{\rm SM}}$. Because the phase
transition belongs to low temperature FOPT, the constraints from the epochs of
Big Bang Nucleosynthesis (BBN), recombination, cosmic microwave
background(CMB) and baryon acoustic oscillations should be considered Izotov
et al. (2014); Aghanim et al. (2020); Riess et al. (2016); Mangano and Serpico
(2011). However, we do not know the temperature of the dark sector after
reheating stage. In our study the value of $\Delta N_{\rm eff}<0.5$ is given
and $\frac{T_{f}}{T_{\rm SM}}=0.47$ is chosen. Subsequently, these
intersection points in Fig. 2 are used to calculate the GW signal (with an
additional input parameter $h$, as mentioned before). The signal-to-noise
ratios (SNRs) are plotted in Fig. 3, which shows that in the future the GW
signal produced by the singlet SUSY model may be detected by LISA, TianQin and
Taiji detectors.

##  V Conclusion

We examined the cosmological phenomena in the minimal supersymmetric standard
model (MSSM) extended by a singlet superfield, with the coupling between the
singlet sector and the MSSM sector being tiny and hence the singlet sector
being a quasi dark sector with SUSY. We obtained the following observations:
(i) In the parameter space of SIDM which solves the small cosmological scale
anomalies, a FOPT can readily happen but requires rather light dark matter
below MeV; (ii) The parameter space required by FOPT and SIDM can be partially
covered by detecting the phase-transition GWs at the near-future projects,
such as LISA, TianQin and Taiji. Therefore, the recently developed GW
astronomy could be a good probe to such a SUSY scenario.

## Acknowledgements

We thank Bin Zhu and Yang Zhang for useful discussions. This work was
supported by the Natural Science Foundation of China under grant numbers
11775012, 12075300 and 11821505, by Peng-Huan-Wu Theoretical Physics
Innovation Center (12047503), by the CAS Center for Excellence in Particle
Physics (CCEPP), by Key R&D Program of Ministry of Science and Technology of
China under number 2017YFA0402204, and by Key Research Program of the Chinese
Academy of Sciences under grant No. XDPB15.

##  Appendix A Feynman rules in the singlet SUSY model

The interaction vertices are given by

$\begin{split}V_{HHH}=&~{}-2\sqrt{2}\kappa(6\kappa
v_{s}+A_{\kappa}+3\mu_{s}),\\\ V_{HAA}=&~{}-2\sqrt{2}\kappa(2\kappa
v_{s}-A_{\kappa}+\mu_{s}),\\\ V_{H\chi\chi}=&~{}-\sqrt{2}\kappa,\\\
V_{A\chi\chi}=&~{}\sqrt{2}i\kappa\gamma^{5}.\end{split}$ (43)

The amplitudes for the singlino DM annihilation processes are given by

  1. 1.

$\chi\chi\to HH$ :

$\displaystyle A(^{3}P_{0})$ $\displaystyle=$ $\displaystyle
2\sqrt{6}v\kappa^{2}\left[\frac{R(3m_{\chi}+A_{\kappa})}{4-R(m_{H})^{2}+iG_{h}}-2\frac{1+R(m_{\chi})}{P_{\chi}}+\frac{4}{3}\frac{\bar{\beta}_{f}^{2}}{P_{\chi}^{2}}\right],$
(44) $\displaystyle A(^{3}P_{2})$ $\displaystyle=$
$\displaystyle-(16/\sqrt{3})v\kappa^{2}\bar{\beta}_{f}^{2}/P_{\chi}^{2}~{}.$
(45)

  2. 2.

$\chi\chi\to AA$ :

$\displaystyle A(^{3}P_{0})$ $\displaystyle=$ $\displaystyle
2\sqrt{6}v\kappa^{2}\left[\frac{R(m_{\chi}-A_{\kappa})}{4-R(m_{H})^{2}+iG_{h}}-2\frac{1-R(m_{\chi})}{P_{\chi}}+\frac{4}{3}\frac{\bar{\beta}_{f}^{2}}{P_{\chi}^{2}}\right],$
(46) $\displaystyle A(^{3}P_{2})$ $\displaystyle=$
$\displaystyle-(16/\sqrt{3})v\kappa^{2}\bar{\beta}_{f}^{2}/P_{\chi}^{2}~{}.$
(47)

  3. 3.

$\chi\chi\to HA$ :

$\displaystyle A(^{1}S_{0})$ $\displaystyle=$
$\displaystyle-4\sqrt{2}\kappa^{2}\frac{R(m_{\chi}-A_{\kappa})}{4-R(m_{H})^{2}}\left(1+\frac{v^{2}}{8}\right)$
(48)
$\displaystyle+8\sqrt{2}\kappa^{2}\frac{R(m_{\chi})}{P_{\chi}}\left[1+v^{2}\left(\frac{1}{8}-\frac{1}{2P_{\chi}}+\frac{\bar{\beta}_{f}^{2}}{3P_{\chi}^{2}}\right)\right]$
$\displaystyle+2\sqrt{2}\kappa^{2}\left(R(m_{A})^{2}-R(m_{H})^{2}\right)\left[1+v^{2}\left(-\frac{1}{8}-\frac{1}{2P_{\chi}}+\frac{\bar{\beta}_{f}^{2}}{3P_{\chi}^{2}}\right)\right],$
$\displaystyle A(^{3}P_{1})$ $\displaystyle=$ $\displaystyle
8v\kappa^{2}\bar{\beta}_{f}^{2}/P_{\chi}^{2}~{},$ (49)

where

$\displaystyle R(m_{X})=\frac{m_{X}}{m_{\chi}},\quad
P_{j}=1+R(m_{j})^{2}-\frac{1}{2}(R(m_{X})^{2}+R(m_{Y})^{2}),\quad
G_{i}=\frac{\Gamma_{i}m_{i}}{m_{\chi}^{2}}.$ (50)

## References

  * Young (2017) B.-L. Young, Front. Phys. (Beijing) 12, 121201 (2017), [Erratum: Front.Phys.(Beijing) 12, 121202 (2017)].
  * Dodelson (2003) S. Dodelson, _Modern Cosmology_ (Academic Press, Amsterdam, 2003), ISBN 978-0-12-219141-1.
  * Linde (1983) A. D. Linde, Nucl. Phys. B216, 421 (1983), [Erratum: Nucl. Phys.B223,544(1983)].
  * Coleman (1977) S. R. Coleman, Phys. Rev. D15, 2929 (1977), [Erratum: Phys. Rev.D16,1248(1977)].
  * Callan and Coleman (1977) C. G. Callan, Jr. and S. R. Coleman, Phys. Rev. D16, 1762 (1977).
  * Amaro-Seoane et al. (2017) P. Amaro-Seoane et al. (LISA) (2017), eprint 1702.00786.
  * Luo et al. (2016) J. Luo et al. (TianQin), Class. Quant. Grav. 33, 035010 (2016), eprint 1512.02076.
  * Ruan et al. (2020) W.-H. Ruan, Z.-K. Guo, R.-G. Cai, and Y.-Z. Zhang, Int. J. Mod. Phys. A 35, 2050075 (2020), eprint 1807.09495.
  * Kajantie et al. (1996) K. Kajantie, M. Laine, K. Rummukainen, and M. E. Shaposhnikov, Phys. Rev. Lett. 77, 2887 (1996), eprint hep-ph/9605288.
  * Fodor et al. (1995) Z. Fodor, J. Hein, K. Jansen, A. Jaster, and I. Montvay, Nucl. Phys. B439, 147 (1995), eprint hep-lat/9409017.
  * Kajantie et al. (1993) K. Kajantie, K. Rummukainen, and M. E. Shaposhnikov, Nucl. Phys. B407, 356 (1993), eprint hep-ph/9305345.
  * Espinosa et al. (2012) J. R. Espinosa, T. Konstandin, and F. Riva, Nucl. Phys. B854, 592 (2012), eprint 1107.5441.
  * Athron et al. (2019a) P. Athron, C. Balazs, A. Fowlie, G. Pozzo, G. White, and Y. Zhang, JHEP 11, 151 (2019a), eprint 1908.11847.
  * Apreda et al. (2002) R. Apreda, M. Maggiore, A. Nicolis, and A. Riotto, Nucl. Phys. B631, 342 (2002), eprint gr-qc/0107033.
  * Hasegawa et al. (2019) T. Hasegawa, N. Okada, and O. Seto, Phys. Rev. D99, 095039 (2019), eprint 1904.03020.
  * Bian et al. (2019) L. Bian, Y. Wu, and K.-P. Xie, JHEP 12, 028 (2019), eprint 1909.02014.
  * Xie et al. (2020) K.-P. Xie, L. Bian, and Y. Wu, JHEP 12, 047 (2020), eprint 2005.13552.
  * Su et al. (2021) W. Su, A. G. Williams, and M. Zhang, JHEP 04, 219 (2021), eprint 2011.04540.
  * Ghosh et al. (2021) T. Ghosh, H.-K. Guo, T. Han, and H. Liu, JHEP 07, 045 (2021), eprint 2012.09758.
  * Guo et al. (2021a) H.-K. Guo, K. Sinha, D. Vagie, and G. White, JHEP 06, 164 (2021a), eprint 2103.06933.
  * Gonçalves et al. (2021) D. Gonçalves, A. Kaladharan, and Y. Wu (2021), eprint 2108.05356.
  * Wess and Zumino (1974) J. Wess and B. Zumino, Nucl. Phys. B 70, 39 (1974).
  * Sohnius (1985) M. F. Sohnius, Phys. Rept. 128, 39 (1985).
  * Sakai (1981) N. Sakai, Z. Phys. C 11, 153 (1981).
  * Dimopoulos and Georgi (1981) S. Dimopoulos and H. Georgi, Nucl. Phys. B 193, 150 (1981).
  * Baer et al. (2020) H. Baer, V. Barger, S. Salam, D. Sengupta, and K. Sinha, Eur. Phys. J. ST 229, 3085 (2020), eprint 2002.03013.
  * Wang et al. (2022) F. Wang, W. Wang, J. Yang, Y. Zhang, and B. Zhu, Universe 8, 178 (2022), eprint 2201.00156.
  * Ellwanger et al. (2010) U. Ellwanger, C. Hugonie, and A. M. Teixeira, Phys. Rept. 496, 1 (2010), eprint 0910.1785.
  * Cao et al. (2012) J.-J. Cao, Z.-X. Heng, J. M. Yang, Y.-M. Zhang, and J.-Y. Zhu, JHEP 03, 086 (2012), eprint 1202.5821.
  * Huang et al. (2015) W. Huang, Z. Kang, J. Shu, P. Wu, and J. M. Yang, Phys. Rev. D 91, 025006 (2015), eprint 1405.1152.
  * Bi et al. (2015) X.-J. Bi, L. Bian, W. Huang, J. Shu, and P.-F. Yin, Phys. Rev. D 92, 023507 (2015), eprint 1503.03749.
  * Bian et al. (2018) L. Bian, H.-K. Guo, and J. Shu, Chin. Phys. C 42, 093106 (2018), [Erratum: Chin.Phys.C 43, 129101 (2019)], eprint 1704.02488.
  * Chatterjee et al. (2022) A. Chatterjee, A. Datta, and S. Roy (2022), eprint 2202.12476.
  * Wang et al. (2014) F. Wang, W. Wang, J. M. Yang, and S. Zhou, Phys. Rev. D 90, 035028 (2014), eprint 1404.6705.
  * Zhu and Abdughani (2021) B. Zhu and M. Abdughani, JHEP 12, 059 (2021), eprint 2103.06050.
  * Navarro et al. (1997) J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophys. J. 490, 493 (1997), eprint astro-ph/9611107.
  * Moore et al. (1999) B. Moore, T. R. Quinn, F. Governato, J. Stadel, and G. Lake, Mon. Not. Roy. Astron. Soc. 310, 1147 (1999), eprint astro-ph/9903164.
  * Kuzio de Naray et al. (2010) R. Kuzio de Naray, G. D. Martinez, J. S. Bullock, and M. Kaplinghat, Astrophys. J. Lett. 710, L161 (2010), eprint 0912.3518.
  * Bullock et al. (2001) J. S. Bullock, T. S. Kolatt, Y. Sigad, R. S. Somerville, A. V. Kravtsov, A. A. Klypin, J. R. Primack, and A. Dekel, Mon. Not. Roy. Astron. Soc. 321, 559 (2001), eprint astro-ph/9908159.
  * Oman et al. (2015) K. A. Oman et al., Mon. Not. Roy. Astron. Soc. 452, 3650 (2015), eprint 1504.01437.
  * Klypin et al. (1999) A. A. Klypin, A. V. Kravtsov, O. Valenzuela, and F. Prada, Astrophys. J. 522, 82 (1999), eprint astro-ph/9901240.
  * Kauffmann et al. (1993) G. Kauffmann, S. D. M. White, and B. Guiderdoni, Mon. Not. Roy. Astron. Soc. 264, 201 (1993).
  * Zavala et al. (2009) J. Zavala, Y. P. Jing, A. Faltenbacher, G. Yepes, Y. Hoffman, S. Gottlober, and B. Catinella, Astrophys. J. 700, 1779 (2009), eprint 0906.0585.
  * Boylan-Kolchin et al. (2012) M. Boylan-Kolchin, J. S. Bullock, and M. Kaplinghat, Mon. Not. Roy. Astron. Soc. 422, 1203 (2012), eprint 1111.2048.
  * Tollerud et al. (2014) E. J. Tollerud, M. Boylan-Kolchin, and J. S. Bullock, Mon. Not. Roy. Astron. Soc. 440, 3511 (2014), eprint 1403.6469.
  * Garrison-Kimmel et al. (2014) S. Garrison-Kimmel, M. Boylan-Kolchin, J. S. Bullock, and E. N. Kirby, Mon. Not. Roy. Astron. Soc. 444, 222 (2014), eprint 1404.5313.
  * Iršič et al. (2017) V. Iršič et al., Phys. Rev. D 96, 023522 (2017), eprint 1702.01764.
  * Viel et al. (2013) M. Viel, G. D. Becker, J. S. Bolton, and M. G. Haehnelt, Phys. Rev. D 88, 043502 (2013), eprint 1306.2314.
  * Menci et al. (2016) N. Menci, A. Grazian, M. Castellano, and N. G. Sanchez, Astrophys. J. Lett. 825, L1 (2016), eprint 1606.02530.
  * Maccio et al. (2012) A. V. Maccio, S. Paduroiu, D. Anderhalden, A. Schneider, and B. Moore, Mon. Not. Roy. Astron. Soc. 424, 1105 (2012), eprint 1202.1282.
  * Spergel and Steinhardt (2000) D. N. Spergel and P. J. Steinhardt, Phys. Rev. Lett. 84, 3760 (2000), eprint astro-ph/9909386.
  * Elor et al. (2021) G. Elor, R. McGehee, and A. Pierce (2021), eprint 2112.03920.
  * Wainwright (2012) C. L. Wainwright, Comput. Phys. Commun. 183, 2006 (2012), eprint 1109.4189.
  * Athron et al. (2020) P. Athron, C. Balázs, A. Fowlie, and Y. Zhang, Eur. Phys. J. C80, 567 (2020), eprint 2003.02859.
  * Athron et al. (2019b) P. Athron, C. Balázs, M. Bardsley, A. Fowlie, D. Harries, and G. White, Comput. Phys. Commun. 244, 448 (2019b), eprint 1901.03714.
  * Guth and Weinberg (1981) A. H. Guth and E. J. Weinberg, Phys. Rev. D 23, 876 (1981).
  * Rintoul and Torquato (1997) M. D. Rintoul and S. Torquato, Journal of physics a: mathematical and general 30, L585 (1997).
  * Megevand and Ramirez (2017) A. Megevand and S. Ramirez, Nucl. Phys. B 919, 74 (2017), eprint 1611.05853.
  * Kobakhidze et al. (2017) A. Kobakhidze, C. Lagger, A. Manning, and J. Yue, Eur. Phys. J. C 77, 570 (2017), eprint 1703.06552.
  * Ellis et al. (2019) J. Ellis, M. Lewicki, and J. M. No, JCAP 04, 003 (2019), eprint 1809.08242.
  * Ellis et al. (2020) J. Ellis, M. Lewicki, and J. M. No, JCAP 07, 050 (2020), eprint 2003.07360.
  * Wang et al. (2020) X. Wang, F. P. Huang, and X. Zhang, JCAP 05, 045 (2020), eprint 2003.08892.
  * Grojean and Servant (2007) C. Grojean and G. Servant, Phys. Rev. D 75, 043507 (2007), eprint hep-ph/0607107.
  * Caprini et al. (2016) C. Caprini et al., JCAP 1604, 001 (2016), eprint 1512.06239.
  * Caprini et al. (2020) C. Caprini et al., JCAP 03, 024 (2020), eprint 1910.13125.
  * Guo et al. (2021b) H.-K. Guo, K. Sinha, D. Vagie, and G. White, JCAP 01, 001 (2021b), eprint 2007.08537.
  * Espinosa et al. (2010) J. R. Espinosa, T. Konstandin, J. M. No, and G. Servant, JCAP 1006, 028 (2010), eprint 1004.4187.
  * Drees and Nojiri (1993) M. Drees and M. M. Nojiri, Phys. Rev. D 47, 376 (1993), eprint hep-ph/9207234.
  * Jungman et al. (1996) G. Jungman, M. Kamionkowski, and K. Griest, Phys. Rept. 267, 195 (1996), eprint hep-ph/9506380.
  * Krstić and Schultz (1999) P. S. Krstić and D. R. Schultz, Phys. Rev. A 60, 2118 (1999).
  * Vanderheyden (2021) L. Vanderheyden (2021), eprint 2107.13845.
  * Tulin et al. (2013a) S. Tulin, H.-B. Yu, and K. M. Zurek, Phys. Rev. Lett. 110, 111301 (2013a), eprint 1210.0900.
  * Tulin et al. (2013b) S. Tulin, H.-B. Yu, and K. M. Zurek, Phys. Rev. D 87, 115007 (2013b), eprint 1302.3898.
  * Ko and Tang (2014) P. Ko and Y. Tang, JCAP 05, 047 (2014), eprint 1402.6449.
  * Tulin and Yu (2018) S. Tulin and H.-B. Yu, Phys. Rept. 730, 1 (2018), eprint 1705.02358.
  * Krylov et al. (2013) E. Krylov, A. Levin, and V. Rubakov, Phys. Rev. D 87, 083528 (2013), eprint 1301.0354.
  * Hong et al. (2020) J.-P. Hong, S. Jung, and K.-P. Xie, Phys. Rev. D 102, 075028 (2020), eprint 2008.04430.
  * Kawana and Xie (2022) K. Kawana and K.-P. Xie, Phys. Lett. B 824, 136791 (2022), eprint 2106.00111.
  * Marfatia and Tseng (2021) D. Marfatia and P.-Y. Tseng, JHEP 11, 068 (2021), eprint 2107.00859.
  * Di Bari et al. (2021) P. Di Bari, D. Marfatia, and Y.-L. Zhou, JHEP 10, 193 (2021), eprint 2106.00025.
  * Izotov et al. (2014) Y. I. Izotov, T. X. Thuan, and N. G. Guseva, Mon. Not. Roy. Astron. Soc. 445, 778 (2014), eprint 1408.6953.
  * Aghanim et al. (2020) N. Aghanim et al. (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], eprint 1807.06209.
  * Riess et al. (2016) A. G. Riess et al., Astrophys. J. 826, 56 (2016), eprint 1604.01424.
  * Mangano and Serpico (2011) G. Mangano and P. D. Serpico, Phys. Lett. B 701, 296 (2011), eprint 1103.1261.


# Tropical cryptography III: digital signatures

Jiale Chen  Department of Mathematics, The City College of New York, New York,
NY 10031<EMAIL_ADDRESS>,  Dima Grigoriev  CNRS, Mathématiques,
Université de Lille, 59655, Villeneuve d’Ascq, France
<EMAIL_ADDRESS>and  Vladimir Shpilrain  Department of
Mathematics, The City College of New York, New York, NY 10031
<EMAIL_ADDRESS>
###### Abstract.

We use tropical algebras as platforms for a very efficient digital signature
protocol. Security relies on computational hardness of factoring one-variable
tropical polynomials; this problem is known to be NP-hard.

Keywords: tropical algebra, digital signature, factoring polynomials

##  1\. Introduction

In our earlier papers [2], [3], we employed tropical algebras as platforms for
cryptographic schemes by mimicking some well-known classical schemes, as well
as newer schemes like [4], [5], in the “tropical” setting. What it means is
that we replaced the usual operations of addition and multiplication by the
operations $\min(x,y)$ and $x+y$, respectively.

An obvious advantage of using tropical algebras as platforms is unparalleled
efficiency because in tropical schemes, one does not have to perform any
multiplications of numbers since tropical multiplication is the usual
addition, see Section 2. On the other hand, “tropical powers” of an element
exhibit some patterns, even if such an element is a matrix over a tropical
algebra. This weakness was exploited in [7] to arrange a fairly successful
attack on one of the schemes in [2].

In this paper, we offer a digital signature scheme that uses tropical algebra
of one-variable polynomials. Security of this scheme is based on computational
hardness of factoring one-variable tropical polynomials. This problem is known
to be NP-hard, see [6].

##  2\. Preliminaries

We start by giving some necessary information on tropical algebras here; for
more details, we refer the reader to the monograph [1].

Consider a tropical semiring ${\bf A}$, also known as the min-plus algebra due
to the following definition. This semiring is defined as a linearly ordered
set (e.g., a subset of reals) that contains 0 and is closed under addition,
with two operations as follows:

$x\oplus y=\min(x,y)$

$x\otimes y=x+y$.

It is straightforward to see that these operations satisfy the following
properties:

associativity:
$x\oplus(y\oplus z)=(x\oplus y)\oplus z$
$x\otimes(y\otimes z)=(x\otimes y)\otimes z$.

commutativity:
$x\oplus y=y\oplus x$
$x\otimes y=y\otimes x$.

distributivity:
$(x\oplus y)\otimes z=(x\otimes z)\oplus(y\otimes z)$.

There are some “counterintuitive” properties as well:

$x\oplus x=x$

$x\otimes 0=x$

$x\oplus 0$ could be either 0 or $x$.

There is also a special “$\epsilon$-element” $\epsilon=\infty$ such that, for
any $x\in S$,

$\epsilon\oplus x=x$

$\epsilon\otimes x=\epsilon$.

###  2.1. Tropical polynomials

A (tropical) monomial in $S$ looks like a usual linear function, and a
tropical polynomial is the minimum of a finite number of such functions, and
therefore a concave, piecewise linear function. The rules for the order in
which tropical operations are performed are the same as in the classical case,
see the example below. Still, we often use parenthesis to make a tropical
polynomial easier to read.

######  Example 1.

Here is an example of a tropical monomial: $x\otimes x\otimes y\otimes
z\otimes z$. The (tropical) degree of this monomial is 5. We note that
sometimes, people use the alternative notation $x^{\otimes 2}$ for $x\otimes
x$, etc.

An example of a tropical polynomial is: $p(x,y,z)=5\otimes x\otimes y\otimes
z\oplus x\otimes x\oplus 2\otimes z\oplus 17=(5\otimes x\otimes y\otimes
z)\oplus(x\otimes x)\oplus(2\otimes z)\oplus 17.$ This polynomial has
(tropical) degree 3, by the highest degree of its monomials.

We note that, just as in the classical case, a tropical polynomial is
canonically represented by an ordered set of tropical monomials (together with
finite coefficients), where the order that we use here is deglex.

We also note that some tropical polynomials may look “weird”:

######  Example 2.

Consider the polynomial $p(x)=(0\otimes x)\oplus(0\otimes x\otimes x)$. All
coefficients in this polynomial are 0, and yet it is not the same as the
polynomial $q(x)=0$.

Indeed, $q(x)\otimes r(x)=r(x)$ for any polynomial $r(x)$. On the other hand,
if, say, $r(x)=2\otimes x$, then $p(x)\otimes r(x)=(2\otimes x\otimes
x)\oplus(2\otimes x\otimes x\otimes x)\neq r(x)$.

In the following example, we show in detail how two tropical polynomials are
multiplied and how similar terms are collected.

######  Example 3.

Let $p(x)=(2\otimes x)\oplus(3\otimes x\otimes x)$ and $q(x)=5\oplus(1\otimes
x)$. Then $p(x)\otimes q(x)=[(2\otimes x)\otimes 5]\oplus[(2\otimes
x)\otimes(1\otimes x)]\oplus[(3\otimes x\otimes x)\otimes 5]\oplus[(3\otimes
x\otimes x)\otimes(1\otimes x)]=(7\otimes x)\oplus(3\otimes x\otimes
x)\oplus(8\otimes x\otimes x)\oplus(4\otimes x\otimes x\otimes x)=(7\otimes
x)\oplus(3\otimes x\otimes x)\oplus(4\otimes x\otimes x\otimes x).$

In this paper, our focus is on one-variable tropical polynomials, although one
can use multivariate tropical polynomials instead.

##  3\. Digital signature scheme description

Let ${\bf T}$ be the tropical algebra of one-variable polynomials over ${\bf
Z}$, the ring of integers.

The signature scheme is as follows.

Private: two polynomials $X,Y\in{\bf T}$ of the same degree $d$, with all
coefficients in the range $[0,r]$, where $r$ is one of the parameters of the
scheme.

Public:

– polynomial $M=X\otimes Y$

– a hash function $H$ (e.g., SHA-512) and a (deterministic) procedure for
converting values of $H$ to one-variable polynomials from the algebra ${\bf
T}$ (see Section 4.1).

Signing a message $m$:

  1. (1)

Apply a hash function $H$ to $m$. Convert $H(m)$ to a polynomial $P$ of degree
$d$ from the algebra ${\bf T}$ using a deterministic public procedure.

  2. (2)

Select two random private polynomials $U,V\in{\bf T}$ of degree $d$, with all
coefficients in the range $[0,r]$. Denote $N=U\otimes V$.

  3. (3)

The signature is the triple of polynomials $(P\otimes X\otimes U,~{}P\otimes
Y\otimes V,~{}N)$.

Verification:

1. The verifier computes the hash $H(m)$ and converts $H(m)$ to a polynomial $P$ of degree $d$ from the algebra ${\bf T}$ using a deterministic public procedure.

2. (a) The verifier checks that the degrees of the polynomials $P\otimes X\otimes U$ and $P\otimes Y\otimes V$ (the first two polynomials in the signature) are both equal to $3d$, and the degree of the polynomial $N$ is equal to $2d$. If not, then the signature is not accepted.

(b) The verifier checks that neither $P\otimes X\otimes U$ nor $P\otimes
Y\otimes V$ is a constant multiple (in the tropical sense) of $P\otimes M$ or
$P\otimes N$. If it is, then the signature is not accepted.

(c) The verifier checks that all coefficients in the polynomials $P\otimes
X\otimes U$ and $P\otimes Y\otimes V$ are in the range $[0,3r]$, and all
coefficients in the polynomial $N$ are in the range $[0,2r]$. If not, then the
signature is not accepted.

3. The verifier computes $W=(P\otimes X\otimes U)\otimes(P\otimes Y\otimes V)$. The signature is accepted if and only if $W=P\otimes P\otimes M\otimes N$.

Correctness is obvious since $W=(P\otimes X\otimes U)\otimes(P\otimes Y\otimes
V)=P\otimes P\otimes(X\otimes Y)\otimes(U\otimes V)=P\otimes P\otimes M\otimes
N$.

######  Remark 1.

Step 2 in the verification algorithm is needed to prevent trivial forgery,
e.g. signing by the triple of polynomials $(P\otimes M,~{}P\otimes N,~{}N)$,
in which case $(P\otimes M)\otimes(P\otimes N)=P\otimes P\otimes M\otimes N$.

######  Remark 2.

Here is how one can check whether or not one given tropical polynomial, call
it $R(x)$, is a constant multiple (in the tropical sense) of another given
tropical polynomial (of the same degree), call it $S(x)$.

Let $r_{i}\in{\mathbf{Z}}$ denote the coefficient at the monomial $x^{\otimes
i}$ in $R(x)$, and $s_{i}\in{\mathbf{Z}}$ denote the coefficient at the
monomial $x^{\otimes i}$ in $S(x)$. If $R(x)=c\otimes S(x)$ for some
$c\in{\mathbf{Z}}$, then $r_{i}=s_{i}+c$ for every $i$. Here “+” means the
“classical” addition in ${\mathbf{Z}}$.

Therefore, to check if $R(x)$ is a constant multiple of $S(x)$, one checks if
$(r_{i}-s_{i})$ is the same integer for every $i$.

##  4\. Key generation and suggested parameters

The suggested degree $d$ of each polynomial $X,Y,U,V$ is 150.

All coefficients of monomials in the polynomials $X,Y,U,V$ are selected
uniformly at random from integers in the range $[0,r]$, where $r=127$. We
emphasize that, in contrast with the “classical” case, if the coefficient at a
monomial is 0, this does not mean that this monomial is “absent” from the
polynomial.

###  4.1. Converting $H(m)$ to a tropical polynomial over ${\mathbf{Z}}$

We suggest using a hash functions from the SHA-2 family, specifically SHA-512.
We assume the security properties of SHA-512, including collision resistance
and preimage resistance. We also assume that there is a standard way to
convert $H(m)$ to a bit string of length 512. Then a bit string can be
converted to a tropical polynomial $P=P(x)$ over ${\mathbf{Z}}$ using the
following ad hoc procedure.

Let $B=H(m)$ be a bit string of length 512. We will convert $B$ to a one-
variable tropical polynomial $P$ of degree $d=150$ over ${\mathbf{Z}}$. We
therefore have to select 151 coefficients for monomials in $P$, and we want to
have these coefficients in the range [0, 127]. With 7 bits for each
coefficient, we need $151\cdot 7=1057$ bits in total.

(1) Concatenate 3 copies of the bit string $B$ to get a bit string of length
1536.

(2) Going left to right, convert 7-bit block $\\#j$ to an integer and use it
as the coefficient at the monomial $x^{\otimes j}$.

(3) After we use $7\cdot 151=1057$ bits, all monomials in the polynomial
$P=P(x)$ will get a coefficient.

###  4.2. Multiplying two tropical polynomials

Let $R(x)$ and $S(x)$ be two one-variable tropical polynomials of degree $d$
and $g$, respectively. We want to compute $R(x)\otimes S(x)$.

Note that a one-variable tropical monomial, together with a coefficient, can
be represented by a pair of integers $(k,l)$, where $k$ is the coefficient and
$l$ is the degree. Our goal is therefore to compute the coefficient at every
monomial of degree from 0 to $d+g$ in the product $R(x)\otimes S(x)$.

Suppose we want to compute the coefficient at the monomial of degree
$m,~{}0\leq m\leq d+g$. Then we go over all coefficients $r_{i}$ at the
monomials of degrees $i\leq m$ in the polynomial $R(x)$ and add (in the
“classical” sense) $r_{i}$ to $s_{j}$, where $s_{j}$ is the coefficient at the
monomial of degree $j=m-i$ in the polynomial $S(x)$.

Having computed all such sums $r_{i}+s_{j}$, we find the minimum among them,
and this is the coefficient at the monomial of degree $m$ in the polynomial
$R(x)\otimes S(x)$.

##  5\. What is the hard problem here?

The (computationally) hard problem that we employ in our construction is
factoring one-variable tropical polynomials. This problem is known to be NP-
hard, see [6].

Since recovering the private tropical polynomials $X$ and $Y$ from the public
polynomial $M=X\otimes Y$ is exactly the factoring problem, we see that
inverting our candidate one-way function $f(X,Y)=X\otimes Y$ is NP-hard.

However, the private tropical polynomials $X$ and $Y$ are involved also in the
signature. For example, from the polynomial $W=P\otimes X\otimes U$ the
adversary can recover $X\otimes U$ because the polynomial $P$ is public. The
polynomial $U$ is private, so it looks like the adversary is still facing the
factoring problem. However, the adversary now knows two products involving the
polynomial $X$, namely $X\otimes Y$ and $X\otimes U$. Therefore, we have a
somewhat different problem here: finding a common divisor of two given
polynomials.

This problem is easy for “classical” one-variable polynomials over
${\mathbf{Z}}$. In particular, any classical one-variable polynomial over
${\mathbf{Z}}$ has a unique factorization (up to constant multiples) as a
product of irreducible polynomials. In contrast, a one-variable tropical
polynomial can have an exponential number of incomparable factorizations [6].
Furthermore, it was shown in [6] that two one-variable tropical polynomials
may not have a unique g.c.d. All this makes it appear likely that the problem
of finding the g.c.d. of two (or more) given one-variable tropical polynomials
is computationally hard. No polynomial-time algorithm for solving this problem
is known. More about this in Section 6.

##  6\. Possible attacks

The most straightforward attack is trying to factor the tropical polynomial
$M=X\otimes Y$ as a product of two tropical polynomials $X$ and $Y$. As we
have pointed out before, this problem is known to be NP-hard [6]. In our
situation, there is an additional restriction on $X$ and $Y$ to be of the same
degree $d$, to pass Step 2 of the verification procedure.

If one reduces the equation $M=X\otimes Y$ to a system of equations in the
coefficients of $X$ and $Y$, then one gets a system of $2d+1$ quadratic
(tropical) equations in $2(d+1)$ unknowns. With $d$ large enough, such a
system is unapproachable; in fact, solving a system of quadratic tropical
equations is known to be NP-hard [9]. The size of the key space for $X$ (or
$Y$) with suggested parameters is $128^{151}=2^{1057}$, so the brute force
search is infeasible.

It is unclear whether accumulating (from different signatures) many tropical
polynomials of the form $M_{i}=X\otimes U_{i}$, with different (still unknown)
polynomials $U_{i}$ can help recover $X$. With each new $M_{i}$, the attacker
gets $d+1$ new unknowns (these are coefficients of $U_{i}$) and $2d+1$ new
equations. There is a well-known trick of reducing a system of quadratic
equations to a system of linear equations by replacing each product of two
unknowns by a new unknown. However, the number of pairs of unknowns increases
by $d^{2}$ with each new $U_{i}$. Therefore, a system of linear equations like
that will be grossly underdetermined, resulting in a huge number of solutions
for the new unknowns, thus making solving the original system (in the old
unknowns) hard, especially given the restrictions on the old unknowns tacitly
imposed by Step 2(c) of the verification procedure.

##  7\. Performance and signature size

For our computer simulations, we used Apple MacBook Pro, M1 CPU (8 Cores), 16
GB RAM computer. Python code is available, see [8].

We note that a one-variable tropical monomial, together with a coefficient,
can be represented by a pair of integers $(k,l)$, where $k$ is the coefficient
and $l$ is the degree of the monomial. Then a one-variable tropical polynomial
of degree 150 is represented by 151 such pairs of integers, by the number of
monomials. If $k$ is selected uniformly at random from integers in the range
$[0,127]$, then the size of such a representation is about 2000 bits on
average. Indeed, 151 coefficient of the average size of 6 bits give about 900
bits. Then, the degrees of the monomials are integers from 0 to 150. These
take up $(\sum_{k=1}^{7}k\cdot(2^{k}-2^{k-1}+1))+8\cdot(150-127)\approx 1000$
bits. Thus, it takes about 2000 bits to represent a single tropical polynomial
with suggested parameters.

Since a private key is comprised of two such polynomials, this means that the
size of the private key in our scheme is about 4000 bits (or 500 bytes) on
average.

The public key is a product of two polynomials of degree 150 and is therefore
a polynomial of degree 300. Coefficients in this polynomial are in the range
$[0,254]$. Using the same argument as in the previous paragraph, we estimate
the size of such polynomial to be about 4500 bits (or 562 bytes) on average.

The signature is a triple of polynomials, two of them are products of 3
tropical polynomials of degree 150 and one of them is a product of 2 tropical
polynomials of degree 150. Therefore, the signature size is about 16,000 bits
(or 2000 bytes) on average.

In the table below, we have summarized performance data for several parameter
sets. Most columns are self-explanatory; the last two columns show memory
usage during verification and during the whole process of signing and
verification.

Performance metrics for various parameter values
---
degree of private polynomials |  range for coefficients in private polynomials |  verification time (sec) |  signature size (Kbytes) |  public key size (Kbytes) |  private key size (Kbytes) |  memory usage, verification (Mbytes) |  memory usage, whole process (Mbytes)
100 |  [0,127] |  ¡0.1 |  1.3 |  0.37 |  0.33 |  0.4 |  0.4
150 |  [0,127] |  0.15 |  2 |  0.56 |  0.5 |  0.37 |  0.5
200 |  [0,127] |  0.25 |  2.6 |  0.74 |  0.67 |  0.47 |  0.6

## References

  * [1] P. Butkovic, Max-linear systems: theory and algorithms, Springer-Verlag London, 2010.
  * [2] D. Grigoriev, V. Shpilrain, Tropical cryptography, Comm. Algebra. 42 (2014), 2624–2632.
  * [3] D. Grigoriev and V. Shpilrain, Tropical cryptography II: extensions by homomorphisms, Comm. Algebra. 47 (2019), 4224–4229.
  * [4] M. Habeeb, D. Kahrobaei, C. Koupparis, V. Shpilrain, Public key exchange using semidirect product of (semi)groups, in: ACNS 2013, Lecture Notes Comp. Sc. 7954 (2013), 475–486.
  * [5] D. Kahrobaei, V. Shpilrain, Using semidirect product of (semi)groups in public key cryptography, in: CiE 2016, Lecture Notes Comp. Sc. 9709 (2016), 132–141.
  * [6] K. H. Kim, F. W. Roush, Factorization of polynomials in one variable over the tropical semiring, https://arxiv.org/pdf/math/0501167.pdf
  * [7] M. Kotov, A. Ushakov, Analysis of a key exchange protocol based on tropical matrix algebra, J. Math. Cryptology 12 (2018), 137–141.
  * [8] Python code for the tropical digital signature scheme, https://shpilrain.ccny.cuny.edu/tropicalDS.txt
  * [9] T. Theobald, On the frontiers of polynomial computations in tropical geometry, J. Symbolic Comput. 41 (2006), 1360–1375.