# Why Topological Data Analysis Detects Financial Bubbles?

Samuel W. Akingbade  Department of Mathematical Sciences, Yeshiva University,
USA<EMAIL_ADDRESS>,  Marian Gidea  Department of Mathematical
Sciences, Yeshiva University, USA<EMAIL_ADDRESS>,  Matteo Manzi
CrunchDAO, France<EMAIL_ADDRESS>and  Vahid Nateghi  Aerospace
Sciences and Technology Department, Politecnico di Milano, Italy
<EMAIL_ADDRESS>
###### Abstract.

We present a heuristic argument for the propensity of Topological Data
Analysis (TDA) to detect early warning signals of critical transitions in
financial time series. Our argument is based on the Log-Periodic Power Law
Singularity (LPPLS) model, which characterizes financial bubbles as super-
exponential growth (or decay) of an asset price superimposed with oscillations
increasing in frequency and decreasing in amplitude when approaching a
critical transition (tipping point). We show that whenever the LPPLS model is
fitting with the data, TDA generates early warning signals. As an application,
we illustrate this approach on a sample of positive and negative bubbles in
the Bitcoin historical price.

##  1\. Introduction

Topological Data Analysis (TDA) has emerged in recent years as a powerful
methodology in time-series analysis and signal processing. A prevalent
approach relies on time-delay coordinate embedding, which is used to
construct, from the time-series, a point-cloud in some high-dimensional
Euclidean space [Tak81]. The dynamics on the point-cloud allows one to recover
the evolution laws driving the time-series. The shape of the point-cloud can
be characterized via Persistent Homology – a fundamental tool in TDA – which
measures the number of holes of a geometric approximation of the point-cloud,
as well as their ‘visibility’ at any given ‘resolution’ level. Applying a
sliding window to a time series allows one to measure how the shape of the
point-cloud associated to each window changes over time. The robustness of
persistent homology of the point-cloud under small perturbations makes this
approach particularly suitable for analyzing noisy time series. The method can
be used, for example, for detection of periodic behavior, attractors, and
bifurcations. More broadly, it can be used for detection of critical
transitions (tipping points). These represent abrupt shifts of a system from
one steady state to another, due to small changes in external conditions.
Critical transitions are relevant to many disciplines, including finance,
climate, ecology, and biomedicine; see, e.g., [SBB+09, Len11, TS11, DCB+12,
BSP+21]. Some applications of TDA to detect critical transitions have been
developed in [BGVJ14, BG14, SSGC+18, IYG20, DUC20].

An important class of applications of TDA concerns critical transitions in
financial time series, particularly market crashes and financial bubbles. A
positive bubble represents a rapid increase of the price of an asset – which
does not reflect the asset’s intrinsic value – followed by a crash. A negative
bubble represents a rapid decay followed by a quick rebound. In this context,
the asset price peak in a positive bubble prior to a crash, and the trough in
a negative bubble prior to a rebound, represent tipping points of the system.
There is a practical interest in the detection of early warning signals of
critical transitions, which would allow minimizing the losses from a bubble
crash, or maximize the profit from a rebound. Some related references include:
[PY11, YWS12, FC16, Gid17, GK18, GGK+20, KB21, RQD21, DR21].

It is important to distinguish between crashes of an endogenous origin,
essentially caused by speculative, unsustainable, accelerating bubbles, and
crashes of an exogenous origin, caused by external factors (e.g., natural
cataclysms, pandemics, fraud, political changes). Exogenous crashes are less
likely to yield early warning signals in comparison to the endogenous ones
[SSZ22, SSZ21].

Most of the evidence, so far, on the ability of TDA to detect financial
bubbles has been empirical. A key question remains: what features of a time
series yield to significant changes in the shape of the associated point
cloud, more precisely, in terms of the holes that appear in the point-cloud?
Some previous works attempted to explain this by drawing analogies between a
market undergoing a crash and bifurcations in a complex dynamical system, or
in terms of changes in drift and volatility of the strong non-stationary
process described by a market near a crash [Gid17, GK18, GGK+20, AKV21, GK23]

In this paper, we propose a heuristic argument for why TDA is able to detect
financial bubbles. As we are interested in early warning signals, we focus on
the time period of a bubble before the tipping point. The ansatz is that the
time series during such a period follows the Log-Periodic Power Law
Singularity (LPPLS) model, which characterizes positive bubbles by super-
exponential growth superimposed with oscillations of increasing frequency and
decreasing amplitude when approaching the tipping point. A similar
characterization holds for negative bubbles. The LPPLS model is based on the
theory of rational expectation, and provides a signature for endogenous
bubbles. See, e.g., [Bla79, BW82, FF96, JSW+96, JLS00, YRWS12, LS13, SS16,
Sor17, GDS18].

From a topological view-point, the time-delay reconstruction of LPPLS
oscillatory signals yields loops (surrounding holes) in the point-clouds, and
changes in the frequency and/or amplitude of the oscillations yield changes in
the structure of holes in the point-cloud, which can be measured via
persistent homology. When approaching the tipping point of a bubble, there is
a significant change in the features of the oscillatory signal and
consequently in the TDA output.

The upshot of our work is that whenever the LPPLS model for the detection of a
positive (negative) bubble applies to some data set, the TDA method on that
data set provides early warning signals of critical transition. As far as we
know, this is the first time when the TDA method is justified in terms of a
deterministic model for the expected dynamics of financial bubbles.

As an illustration, we apply both the LPPLS model and the TDA to Bitcoin – the
largest cryptocurrency by market cap. Bitcoin has been experiencing numerous
phases of extreme price growth and massive crashes, hence it provides an
excellent test bed for detection of financial bubbles. Some related references
include [GDS18, ADGK18, RCPB19, GGK+20, ALGK20, SZ20b, SNC20].

Our experiments show a good agreement between the TDA method and the LPPLS
model. Moreover, even when the LPPLS model fits poorly with some of the data,
the TDA method can show relatively strong signal prior to the tipping point.
While the TDA method is robust, the LPPLS appears to be very sensitive to
noise [BCP13].

##  2\. Background

###  2.1. Log-Periodic Power Law Singularity (LPPLS) Model

Below we summarize the LPPLS model following [JSW+96, JLS00, YRWS12, SS16,
Sor17, GDS18]. This is a model for both positive bubbles as well as for
negative bubbles. In either case the market reaches a point when it no longer
reflects the real underlying value of the asset (or it violates the market
efficiency with respect to its time-scale [SS17]), which triggers a critical
transition.

The LPPLS model can be derived from a network description of the market, which
assumes that the trading of the asset is driven by two types of agents: a
group consisting of traders with rational expectations, and another group
formed by noise traders. Noise traders are likely to deviate from fundamental
valuation in an accelerating mode, due to positive feedback and herding
behaviour as a group. The collective behavior of the agents eventually leads
to market instability by pushing the market away from an equilibrium between
supply and demand. See [JLS00, SH05, SS16, Sor17].

It is important to keep in mind that the LPPLS model is concerned only with
bubble generated endogenously by the dynamics of the market, and does not
apply to bubbles of exogenous origins produced by external shocks.

According to the model, the asset price $p(t)$ during a (positive or negative)
bubble evolves according to the stochastic differential equation

(2.1)
$\begin{split}\frac{\textrm{d}p}{p}=&B^{\prime}(t_{c}-t)^{m-1}+C^{\prime}(t_{c}-t)^{m-1}\cos(\omega\ln(t_{c}-t)-\phi^{\prime})+\sigma(t)\textrm{d}W(t)\\\
\end{split}$

where $t_{c}$ represents the critical time,
$B^{\prime},C^{\prime},\phi^{\prime}\in\mathbb{R}$, $0<m<1$ are parameters,
$\sigma(t)$ is the volatility, and $dW(t)$ is the increment of a standard
Wiener process with zero mean and unit variance. Hence the expectation of the
logarithm of the asset price satisfies

(2.2) $\begin{split}E[\ln
p(t)]=&A+B(t_{c}-t)^{m}+C(t_{c}-t)^{m}\cos(\omega\ln(t_{c}-t)-\phi)\\\
=&A+B(t_{c}-t)^{m}+C_{1}(t_{c}-t)^{m}\cos(\omega\ln(t_{c}-t))\\\
&\qquad\qquad\qquad\,\,\,+C_{2}(t_{c}-t)^{m}\sin(\omega\ln(t_{c}-t)).\end{split}$

Here $A>0$, $0<m<1$, $B,C$ and $\phi$ are parameters, whose role is explained
below. Also, $C_{1}=C\cos(\phi)$, $C_{2}=-C\sin(\phi)$, and
$C=\sqrt{C_{1}^{2}+C_{2}^{2}}$.

The function in (2.2) has a singularity at $t=t_{c}$. Since

$\lim_{t\to t_{c}}(t_{c}-t)^{m}\cos(\omega\ln(t_{c}-t))=0,$

the function can be extended by continuity at $t=t_{c}$ by $E[\ln
p(t_{c})]=A$. The critical time $t_{c}$ marks the termination of the bubble
regime and the transition time to a different regime, which could be a be a
major crash or rebound. The parameter $A$ represents the expected value of the
log price $\ln p(t)$ at the critical time $t_{c}$.

The condition $m<1$ ensures that there is a singularity at the critical time
$t_{c}$, while $m>0$ ensures that the expected price remains finite at
$t_{c}$. The parameter $m$ is associated to the non-linearity of the trend,
that is, $m\approx 0$ corresponds to a strong super-linear trend of $\ln
p(t)$, while $m\approx 1$ corresponds to a nearly linear trend. The parameter
$m$ is also affecting, jointly with the parameters $C_{1}$ and $C_{2}$, the
changing amplitude of the oscillations.

The case when $B<0$ corresponds to a super-exponential growth of the price
$p(t)$ as the time $t$ approaches $t_{c}$, while $B>0$ corresponds to a super-
exponential decay.

Thus, the main features of the model (2.1) are the following

  * •

the price of the asset oscillates about a trend-line that is characterized by
super-exponential growth or decay;

  * •

the amplitude of the oscillations decreases when approaching the critical time
$t_{c}$;

  * •

the frequency of the oscillations increases when approaching the critical time
$t_{c}$;

  * •

the oscillations are superimposed with Gaussian noise.

$\begin{array}[]{cc}\includegraphics[width=138.76157pt]{Figures/LPPLS_synthetic}&\includegraphics[width=138.76157pt]{Figures/LPPLS_synthetic_noise}\\\
\end{array}$

Figure 1. Synthetic LPPLS signal: without noise (left), with noise (right).

The LPPLS model can be fitted to data and allows for predictions of the time
of the crash [YRWS12, GDS18]. However, we note that the LPPLS model is fairly
sensitive to the system parameters, to where to start and end the fit of a
bubble, and to the length of the data window (see [BCP13] and the references
listed there).

###  2.2. Persistent homology

In this section we summarize the main topological tools that we will use:
persistent homology and persistence landscapes. See [EH22, B+15, BD17].
Persistent homology converts a given point cloud into a geometric shape and
characterizes it in terms of its homology generators. Persistence landscapes
convert the information derived from persistent homology into elements of a
functional Banach space. Through this process, a point cloud can be
characterized by a single numerical value, which is the norm of the
corresponding persistence landscape.

####  2.2.1. Point-cloud

The input of persistent homology is a point-cloud
$Z=\\{z_{0},\ldots,z_{m-1}\\}$ embedded in some Euclidean space
$\mathbb{R}^{N}$. The computation of persistent homology involves a sequence
of steps.

####  2.2.2. Vietoris-Rips simplicial complex.

We associate to the point cloud $Z$ a topological space, thought of as a
‘shape’, depending on some parameter $\varepsilon>0$, where $\varepsilon$
represents the resolution at which the data is analyzed. The Vietoris-Rips
simplicial complex $R(Z,\varepsilon)$, or, simply Rips complex, is defined as
follows:

  * •

For each $k=0,1,\ldots$, a $k$-simplex of vertices
$\\{z_{i_{0}},\ldots,z_{i_{k}}\\}$ is part of the simplicial complex
$R(Z,\varepsilon)$ if and only if the mutual distance between any pair of its
vertices is less than $\varepsilon$, that is

$d(z_{i_{j}},z_{i_{j^{\prime}}})<\varepsilon,\textrm{ for all
}z_{i_{j}},z_{i_{i_{j^{\prime}}}}\in\\{z_{i_{0}},\ldots,z_{i_{k}}\\}.$

In other words, a $k$-simplex is included in $R(Z,\varepsilon)$ whenever the
vertices of that simplex are ‘indistinguishable from one another’ at a
resolution level $\varepsilon$.

####  2.2.3. Filtration of simplicial complexes.

The Rips simplicial complexes $R(Z,\varepsilon)$ form a filtration with
respect to the resolution parameter $\varepsilon$, that is,

$R(Z,\varepsilon)\subseteq R(Z,\varepsilon^{\prime})\textrm{ provided
}\varepsilon<\varepsilon^{\prime}.$

####  2.2.4. Homology groups of simplicial complexes.

At each resolution level $\varepsilon$, the corresponding shape can be
characterized in terms of its homology.

  * •

For each complex $R(Z,\varepsilon)$, we compute its $n$-dimensional homology
$H_{n}(R(Z,\varepsilon))$ with coefficients in some field, e.g.,
$\mathbb{Z}_{2}$.

The generators of the $0$-dimensional homology module
$H_{0}(R(Z,\varepsilon))$ correspond to the connected components of
$R(Z,\varepsilon)$, the generators of the $1$-dimensional homology module
$H_{1}(R(Z,\varepsilon))$ correspond to ‘holes’ in $R(Z,\varepsilon)$, the
generators of the $2$-dimensional homology module $H_{2}(R(Z,\varepsilon))$
correspond to ‘voids’ in $R(Z,\varepsilon)$, etc. For the rest of the paper we
will use only $1$-dimensional homology.

####  2.2.5. Filtration of homology modules.

The filtration property of the Rips complexes induces a filtration on the
corresponding homologies, that is

$H_{n}(R(Z,\varepsilon))\subseteq H_{n}(R(Z,\varepsilon^{\prime}))\textrm{
provided }\varepsilon<\varepsilon^{\prime}\,\textrm{ for each }n.$

These inclusions determine canonical homomorphisms

$H_{n}(R(Z,\varepsilon))\hookrightarrow
H_{n}(R(Z,\varepsilon^{\prime})),\textrm{ for
}\varepsilon<\varepsilon^{\prime}.$

####  2.2.6. Persistent homology.

For each non-zero $n$-dimensional homology generator $\alpha$ there exists a
pair of values $\varepsilon_{1}<\varepsilon_{2}$, such that:

  * •

$\alpha$ is in $H_{n}(R(Z,\varepsilon_{1}))$ but is not in the image of any
element in $H_{n}(R(Z,\varepsilon_{1}-\delta))$ under the corresponding
homomorphism, for $\delta>0$;

  * •

the image of $\alpha$ in $H_{n}(R(Z,\varepsilon^{\prime}))$ is non-zero for
all $\varepsilon_{1}<\varepsilon^{\prime}<\varepsilon_{2}$, but the image of
$\alpha$ in $H_{n}(R(Z,\varepsilon_{2}))$ is zero.

In this case, one says that the generator $\alpha$ is ‘born’ at the parameter
value $b_{\alpha}:=\varepsilon_{1}$, and ‘dies’ at the parameter value
$d_{\alpha}:=\varepsilon_{2}$. The pair $(b_{\alpha},d_{\alpha})$ represents
the ‘birth’ and ‘death’ indices of $\alpha$. The multiplicity
$\mu_{\alpha}(b_{\alpha},d_{\alpha})$ of the point $(b_{\alpha},d_{\alpha})$
equals the number of classes $\alpha$ that are born at ${b_{\alpha}}$ and die
at ${d_{\alpha}}$. This multiplicity is finite since the simplicial complex is
finite.

####  2.2.7. Persistence diagram.

The information on the $n$-dimensional homology generators at all resolution
levels can be encoded in a _persistence diagram_ $P_{n}$. Such a diagram
consists of:

  * •

for each $n$-dimensional homology class $\alpha$ one assigns a point
$p_{\alpha}=p_{\alpha}(b_{\alpha},d_{\alpha})\in\mathbb{R}^{2}$ together with
its multiplicity $\mu_{\alpha}=\mu_{\alpha}(b_{\alpha},d_{\alpha})$;

  * •

in addition, $P_{n}$ contains all points on the positive diagonal of
$\mathbb{R}^{2}$; these points represent all trivial homology generators that
are born and instantly die at every level; each point on the diagonal has
infinite multiplicity.

The axes of a persistence diagram are birth values on the horizontal axis and
death values on the vertical axis.

####  2.2.8. Persistence landscapes.

The space of persistence diagrams can be embedded into a Banach space. One
such embedding is based on _persistence landscapes_ , consisting of sequences
of functions in the Banach space $L^{p}(\mathbb{N}\times\mathbb{R})$, with
$p\geq 1$.

  * •

For each birth-death point $(b_{\alpha},d_{\alpha})\in P_{n}$, we first define
a piecewise linear function

(2.3)
$f_{(b_{\alpha},d_{\alpha})}(x)=\left\\{\begin{array}[]{ll}x-b_{\alpha},&\hbox{if
$x\in\left(b_{\alpha},\frac{b_{\alpha}+d_{\alpha}}{2}\right]$;}\\\
-x+d_{\alpha},&\hbox{if
$x\in\left(\frac{b_{\alpha}+d_{\alpha}}{2},d_{\alpha}\right)$;}\\\ 0,&\hbox{if
$x\not\in(b_{\alpha},d_{\alpha})$.}\end{array}\right.$

  * •

To a persistence diagram $P_{n}$ consisting of a finite number of off-diagonal
points, we associate a sequence of functions
$\lambda_{n}=(\lambda_{i})_{i\in\mathbb{N}}$, where
$\lambda_{n}(i):\mathbb{R}\to[0,+\infty]$ is given by

(2.4)
$\lambda_{n}(i)(x)=i\textrm{-max}\\{f_{(b_{\alpha},d_{\alpha})}(x)\,|\,(b_{\alpha},d_{\alpha})\in
P_{n}\\}$

where $i\textrm{-max}$ denotes the $i$-th largest value of a function. We set
$\lambda_{n}(i)(x)=0$ if the $i$-th largest value does not exist. Thus, the
persistence landscapes form a subset of the Banach space
$L^{p}(\mathbb{N}\times\mathbb{R})$.

####  2.2.9. Norm of a persistence landscape.

For each persistence landscape we can compute its norm, and so we can compare
persistence landscapes by computing the mutual distances, using the metric
derived from the norm.

  * •

The norm of $\lambda_{n}\in L^{p}(\mathbb{N}\times\mathbb{R})$ is given by

(2.5)
$\|\lambda_{n}\|_{p}=\left(\sum_{i=1}^{\infty}\|\lambda_{n}(i)\|_{p}^{p}\right)^{1/p}.$

Above, $\|\cdot\|_{p}$ denotes the $L^{p}$-norm, $p\geq 1$, i.e.,
$\|f\|_{p}=\left(\int_{\mathbb{R}}|f|^{p}\right)^{1/p}$, where the integration
is with respect to the Lebesgue measure on $\mathbb{R}$.

We recall from [B+15, BD17] that the norms $\|\cdot\|_{1}$ and
$\|\cdot\|_{\infty}$ can be computed via direct formulas

(2.6)
$\begin{split}\|\lambda_{n}\|_{1}=&\frac{1}{4}\sum_{\alpha}(d_{\alpha}-b_{\alpha})^{2},\\\
\|\lambda_{n}\|_{\infty}=&\frac{1}{2}\sup_{\alpha}(d_{\alpha}-b_{\alpha}).\end{split}$

Throughout the paper we will only refer to the $L^{1}$-norms of persistence
diagrams.

####  2.2.10. Stability

Persistence diagrams and persistence landscapes are ‘stable’ under small
perturbations of a point-cloud [CDSO14, B+15].

###  2.3. Topological data analysis of time-series

As input, consider a time series of numerical values of length $L$:

(2.7) $X=\\{x_{0},x_{1},\ldots,x_{L-1}\\}.$

####  2.3.1. Time-series of time-delay coordinate vectors.

We choose and fix an embedded dimension $N$ and a time-delay $d\geq 1$.

  * •

We transform the time series $X$ into a sequence of $N$-dimensional time-delay
coordinate vectors

(2.8) $\begin{split}z_{0}&=(x_{0},x_{d},\ldots,x_{(N-1)d}),\\\
z_{1}&=(x_{1},x_{1+d},\ldots,x_{1+(N-1)d}),\\\ &\cdots\\\
z_{t}&=(x_{t},x_{t+d},\ldots,x_{t+(N-1)d}),\\\ &\cdots\\\
z_{L-1-(N-1)d}&=(x_{L-1-(N-1)d},x_{L-1-(N-2)d},\ldots,x_{L-1}).\end{split}$

In our applications below, we will only be interested in detecting
$D=1$-dimensional objects (loops), in which case it will be sufficient to
choose $N\geq 2D+1=3$. The theoretical considerations behind the choice of the
embedding dimension $N$ rely on Takens’ embedding theorem [Tak81] and its
generalization by Sauer, Yorke and Casdagli [TJM91]. Also, in our applications
we choose the delay $d$ empirically. If we choose $d$ too small, then the
coordinates $x_{t+id}$ and $x_{t+(i+1)d}$ can be too close to each other, and
so they do not represent two ‘independent’ coordinates in a statistical sense.
Similarly, if $d$ is too large, then $x_{t+id}$ and $x_{t+(i+1)d}$ are
completely independent of each other, and so geometric information gets lost.
Some empirical rules on choosing a suitable delay $d$ are discussed in
[ABST93].

####  2.3.2. Sliding windows.

To detect qualitative changes along a time series, we apply a sliding window
and assess how its features change along the sliding window. In our case:

  * •

We apply a sliding window of size $w$ to the sequence (2.8), for $w$
sufficiently large with $N\ll w\ll L$, obtaining a time-varying point-cloud
embedded in $\mathbb{R}^{N}$:

(2.9) $Z^{t}=\\{z_{t},z_{t+1},\ldots,z_{t+w-1}\\},\textrm{ for
}t\in\\{0,\ldots,K-1\\},$

where we denote $K-1=L-1-(N-1)d-(w-1)$.

####  2.3.3. Time series of norms of persistent landscapes.

Fix $n\geq 1$ and $p\geq 1$.

  * •

For each point-cloud $Z^{t}$ we compute the norm of the persistent landscape
$\|\lambda^{t}_{n}\|_{p}$ , thus obtaining a time-series of norms

$\|\lambda^{t}_{n}\|_{p},\textrm{ for }t\in\\{0,\ldots,K-1\\}.$

##  3\. TDA applied to oscillatory time-series

The main question is to understand what features of a time-series yield
$1$-dimensional ‘holes’ in the persistence diagram and result in ‘peaks’ of
the norms of the persistence landscapes.

Time-series generated by periodic and quasi-periodic signals are some basic
type of signals that result in ‘holes’ in the persistence diagram. We will
also consider oscillatory time series whose frequency changes over time.

###  3.1. Periodic time-series

First, consider time-series generated by a periodic signal of the form

$X(t)=A_{1}\sin(\omega_{1}t+\phi_{1})+\sigma G(t)$

where $A_{1},\omega_{1},\phi_{1}$ are the parameters of the deterministic
signal, $\sigma$ is the noise intensity, and $G(t)$ is Gaussian noise with
zero mean and unit variance. When $\sigma=0$, the embedding of $X(t)$ in
$\mathbb{R}^{N}$ with $N\geq 2$ fills densely a topological circle, which is a
homeomorphic copy of $\mathbb{T}^{1}$. It has a single hole, which is
$1$-dimensional, thus the $1$-dimensional homology is
$H_{1}(\mathbb{T}^{1})=\mathbb{Z}_{2}$ whose generator is a $1$-dimensional
loop. When we add noise with intensity $\sigma$ small, we obtain a ‘noisy’
circle.

To apply the TDA procedure, we first discretize the time $t$, i.e., we let
$t_{i+1}=t_{i}+\Delta t$, with a step size of $\Delta t>0$ small. Then we
choose the window size $w$ and the delay $d$. It is important that the size of
the sliding window $w$ is larger than the period $T_{1}=2\pi/\omega_{1}$ of
the signal, and that the delay $d$ is less than $2\pi/(\omega_{1}\Delta T)$
but not too small. Otherwise, if $T_{1}>w$, then the embedding of a window
does not yield a closed curve. If $d$ is too small, the embedding of the
window may yield a point cloud with a narrow ‘hole’. Indeed, if $\Delta t$ is
small and $d$ is small then the components of
$z_{t}=(x_{t},x_{t+d},\ldots,x_{t+(N-1)d})$ are very close to one another, and
thus the points $z_{t}$ are all located within a small neighborhood of the
diagonal set
$\\{\textbf{x}=(x_{1},\ldots,x_{N})\in\mathbb{R}^{N}\,:\,x_{1}=\ldots=x_{N}\\}$.
In such a case, the persistence of the homology generator for the topological
circle may be relatively small, and hence the corresponding norm of the
persistence landscape may be dominated by the effect of the noise. Some
precise formulations on how to choose the window size $w$ and the delay $d$
can be found in [PH15]. Related works on applying TDA to periodic signals can
be found in [KM18, DQR19].

An example of a periodic time-series without noise is shown in Fig. 2. We note
that delay $d=5$ yields a stronger signal than $d=1$ in both the persistence
diagram and the norm of the persistence landscape.

In either case, the norms of the persistence landscapes stay nearly constant,
since the reconstructed set in $\mathbb{R}^{N}$ – a topological circle – is
the same for each sliding window. The deviations from a constant value are due
to the discretization of the signal.

An example with a periodic time-series with noise is shown in Fig. 3; in this
case we add noise with $\sigma=0.1$. We note that delay $d=5$ still yields a
stronger signal than $d=1$. We also notice that adding the noise weakens the
TDA signal from the no-noise case. This is because the noise narrows the
‘hole’ in the point cloud.

$\begin{array}[]{cccc}\includegraphics[width=95.39693pt]{Figures/circle2d_delay1}&\includegraphics[width=108.405pt]{Figures/circle3d_delay1}&\includegraphics[width=95.39693pt]{Figures/PD_periodic_delay1}&\includegraphics[width=104.07117pt]{Figures/PL_periodic_delay1}\\\
\includegraphics[width=95.39693pt]{Figures/circle2d_delay5}&\includegraphics[width=108.405pt]{Figures/circle3d_delay5}&\includegraphics[width=95.39693pt]{Figures/PD_periodic_delay5}&\includegraphics[width=104.07117pt]{Figures/PL_periodic_delay5}\\\
\end{array}$

Figure 2. Periodic signal without noise reconstructed in 2D and 3D,
persistence diagram and norms of persistence landscapes: $d=1$ (top), and
$d=5$ (bottom)

$\begin{array}[]{cccc}\includegraphics[width=95.39693pt]{Figures/circle2d_noisy_delay1}&\includegraphics[width=108.405pt]{Figures/circle3d_noisy_delay1}&\includegraphics[width=95.39693pt]{Figures/PD_periodic_noisy_delay1}&\includegraphics[width=95.39693pt]{Figures/PL_periodic_noisy_delay1}\\\
\includegraphics[width=95.39693pt]{Figures/circle2d_noisy_delay5}&\includegraphics[width=108.405pt]{Figures/circle3d_noisy_delay5}&\includegraphics[width=95.39693pt]{Figures/PD_periodic_noisy_delay5}&\includegraphics[width=104.07117pt]{Figures/PL_periodic_noisy_delay5}\end{array}$

Figure 3. Periodic signal with noise reconstructed in 2D and 3D, persistence
diagram and norms of persistence landscapes: $d=1$ (top), and $d=5$ (bottom)

###  3.2. Quasi-periodic time-series

Second, consider time-series generated by a quasi-periodic signal with two
incommensurate frequencies

$X(t)=A_{1}\sin(\omega_{1}t+\phi_{1})+A_{2}\sin(\omega_{2}t+\phi_{2})+\sigma
G(t)$

where $A_{1},\omega_{1},\phi_{1},A_{2},\omega_{2},\phi_{2}$ are the parameters
of the deterministic signal, and $\sigma$, $G(t)$ are as before; we assume
that $\omega_{1}/\omega_{2}$ is an irrational number. When $\sigma=0$, for
suitable window size $w$ and delay $d$, the embedding of $X(t)$ in
$\mathbb{R}^{N}$ with $N\geq 3$ yields a topological $2$-dimensional torus
$\mathbb{T}^{2}$. This has two $1$-dimensional homology generators, since
$H_{1}(\mathbb{T}^{2})=\mathbb{Z}^{2}_{2}$. In addition, there is a single
$2$-dimensional homology generator, i.e., $H_{2}(\mathbb{T}^{2})=\mathbb{Z}$.
See Fig. 4.

The norms of the persistence landscapes stay nearly constant, since the
reconstructed set in $\mathbb{R}^{N}$ – a topological torus – is the same for
each sliding window. The deviations from a constant value are due to the
discretization of the signal.

$\begin{array}[]{cccc}\includegraphics[width=95.39693pt]{Figures/torus_2d_delay_5}&\includegraphics[width=108.405pt]{Figures/torus_3d_delay_5}&\includegraphics[width=95.39693pt]{Figures/PD_quasi_freq_2_delay_5}&\includegraphics[width=104.07117pt]{Figures/PL_quasi_freq_2_delay_5}\\\
\end{array}$

Figure 4. Quasi-periodic signal without noise reconstructed in 2D and 3D,
persistence diagram and norms of persistence landscapes: $d=1$ (top), and
$d=5$ (bottom)

Similarly, we can consider general quasi-periodic signals with $m\geq 1$
incommensurate frequencies

$X(t)=\sum_{i=1}^{m}A_{i}\sin(\omega_{i}t+\phi_{i})+\sigma G(t).$

There are specific conditions on how to choose the windows size $w$, the delay
$d$ and the embedding dimension $N$; see [PH15, Per16]. The reconstructed set
is an $m$-dimensional torus $\mathbb{T}^{m}$, whose homology is given by
$H_{k}(\mathbb{T}^{m})=\mathbb{Z}_{2}^{m\choose k}$.

###  3.3. Oscillatory time-series with changing frequency

Another model that produces ‘holes’ in point-cloud is an oscillatory signal of
the form

$X(t)=A_{1}\sin(\omega_{1}\ln(t_{c}-t)+\phi_{1})+\sigma G(t),$

where $t_{c}$ represents a critical time and $t\in[0,t_{c})$. Note that when
$t\to t_{c}$ the frequency of the oscillations
$\omega_{1}\ln(t_{c}-t)\to\infty$ at a logarithmic rate, while the amplitude
of the oscillations stays constant. This type of signal is related to the
LPPLS model (2.2).

There is a strong dependence of the TDA signal on the range $[0,t_{c})$ of the
time-series, on the parameter $\omega_{1}$, and on size of the sliding window
$w$. Note that at any moment of time $t$ the approximate period of the signal
is $T=\frac{2\pi}{\omega_{1}\ln(t_{c}-t)}$. Each full oscillation in the time-
series produces a hole in point-cloud. As $t$ increases, more holes appear in
point-cloud for the corresponding sliding window. When $w$ is larger than the
initial period $T=\frac{2\pi}{\omega_{1}\ln(t_{c})}$, then the TDA signal will
capture the first oscillation in the time-series as well as the subsequent
ones. However, when $w<\frac{2\pi}{\omega_{1}\ln(t_{c})}$, the TDA signal will
only capture oscillations in the time-series whose period
$\frac{2\pi}{\omega_{1}\ln(t_{c}-t)}<w$.

In Fig. 5 we show an oscillatory signal as above with no noise, and its TDA
signal for various sizes of the sliding window.

$\begin{array}[]{ccc}\includegraphics[width=138.76157pt]{Figures/oscillatory}&\includegraphics[width=138.76157pt]{Figures/oscillatory_2D}&\includegraphics[width=138.76157pt]{Figures/oscillatory_3D}\\\
\includegraphics[width=138.76157pt]{Figures/TDA_oscillatory_w100}&\includegraphics[width=138.76157pt]{Figures/TDA_oscillatory_w150}&\includegraphics[width=138.76157pt]{Figures/TDA_oscillatory_w230}\\\
\end{array}$

Figure 5. Oscillatory time-series with changing frequency and its
reconstructions in 2D and 3D (top), norms of persistence landscapes for
$w=100,150,230$ (bottom).

##  4\. TDA applied to the LPPLS model

We recall that for the LPPLS model, the price of the asset oscillates with
increasing frequency and decreasing amplitude about a super-exponential trend
line. When we apply the TDA procedure with a fixed window size we observe the
following:

  * •

At the beginning of the time-series, the behavior is nearly periodic, and the
TDA output is similar to that of a periodic signal, producing a single,
persistent $1$-loop dimensional loop in the point cloud;

  * •

As the time increases, the time-series exhibits more than one frequency, and
TDA output produces a few persistent $1$-dimensional loops in the point cloud;

  * •

As the time approaches the critical time, the time-series exhibits
oscillations of growing frequencies, and the TDA output yields more and more
persistent $1$-dimensional loops, which decay in size. In addition, the super-
exponential trend makes the points in the point-cloud more spread out, so the
growing-in-size ‘gaps’ in the point-cloud contribute with additional
$1$-dimensional loops. Consequently, the norms of the persistent landscapes
spike in the proximity of the critical time.

To capture such behavior, the size of the sliding window $w$ should be
suitably chosen; this depends on the length of the data set as well as on the
parameters of the LPPLS model. For instance, using a larger window size $w$
will result in the larger loops in the earlier part of the time series to
carry more weight towards the TDA signal, which will result in a spike earlier
in the range as opposed to the end of the range. On the other hand, a smaller
window size $w$ will not capture the larger loops but only the smaller ones in
the latter part of the time series, which will result in a spike towards the
end of the range.

In practice, one needs to explore a range of window size $w$ and observe the
presence of a spike in the TDA signal that appears earlier in the range for
large $w$, and shifts towards the end of the range for smaller $w$.

As an example, we consider a noisy LPPLS model for a positive bubble, with the
following parameters

(4.1) $\begin{split}t_{c}=637,\,m=3.003e-01,\,\omega=6.889e+00,\\\
A=1.111e+01,\,B=-2.937e-04\\\
C=5.515e-05,C_{1}=4.372e-05,\,C_{2}=-3.362e-05.\end{split}$

These parameters corresponds to fitting the LPPLS model to the Bitcoin price
between the dates 2021/03/22/–2021/04/16, sampled hourly.

We use the following TDA sliding-window size $w$, delay $d$, and embedding
dimension $N$:

(4.2) $w=72,\,d=5,\,N=4.$

The TDA output for various levels of added noise is shown in Fig. 6. We
observe a peak in the TDA signal prior to the critical time, and note that the
growing trend towards the peak appears earlier when the noise intensity is
larger.

$\begin{array}[]{ccc}\includegraphics[width=138.76157pt]{Figures/Model3_03_22_2021}&\includegraphics[width=138.76157pt]{Figures/LPPLS_model_3_noisy_neg_B}&\includegraphics[width=138.76157pt]{Figures/LPPLS_model_3_more_noisy_neg_B}\\\
\includegraphics[width=138.76157pt]{Figures/TDA_Model_03_22_2021}&\includegraphics[width=138.76157pt]{Figures/LPPLS_model_3_TDA_win_72_delay_5_noisy_neg_B}&\includegraphics[width=138.76157pt]{Figures/LPPLS_model_3_TDA_win_72_delay_5_more_noisy_neg_B}\end{array}$

Figure 6. LPPLS model (top) and norm of persistence landscapes (bottom) for a
positive bubble and increasing noise intensity

###  4.1. Dependence of the TDA signal on parameters

We consider synthetic time series generated by (2.2), of $200$ data points and
with $t_{c}=200$. We start with the following model parameters $m=0.3$,
$\omega=6.7$, $A=11$, $B=-0.0003$, $C_{1}=0.000044$, $C_{2}=-0.000034$, and
with a sliding widow of size $w=48$. These parameters are similar to the ones
in (4.1), except for the length of the time series, which is chosen shorter
for computational efficiency. In the experiments below we explore the changes
in the TDA signal when we change the parameters $w,\omega,m,B,C_{1},C_{2}$,
one at a time. A more detailed explorations would be needed to study the
interdependence among these parameters, which appears to be quite intricate.

####  4.1.1. Dependence on sliding window size $w$.

The experiments displayed in Fig. 7 show that as the window size is decreased,
the peak of the TDA signal gets shifted towards the end of the range of the
time series. This is because windows of larger size $w$ capture more of the
larger loops in the point-cloud, while windows of smaller size capture more of
the smaller loops in the point-cloud. This is consistent with the remarks in
Section 3.

$\begin{array}[]{cc}\includegraphics[width=138.76157pt]{Figures/synthetic_1}&\includegraphics[width=138.76157pt]{Figures/TDA_synthetic_1_w72}\\\
\includegraphics[width=138.76157pt]{Figures/TDA_synthetic_1_w60}&\includegraphics[width=138.76157pt]{Figures/TDA_synthetic_1_w48}\end{array}$

Figure 7. Synthetic LPPLS signal and corresponding TDA signal for window sizes
$w=72,60,48$.

####  4.1.2. Dependence on frequency $\omega$

The experiments displayed in Fig. 8 show that as the frequency $\omega$ is
increased, the TDA signal increases in strength.

$\begin{array}[]{ccc}\includegraphics[width=138.76157pt]{Figures/synthetic_o5}&\includegraphics[width=138.76157pt]{Figures/synthetic_o7}&\includegraphics[width=138.76157pt]{Figures/synthetic_o9}\\\
\includegraphics[width=138.76157pt]{Figures/TDA_synthetic_o5}&\includegraphics[width=138.76157pt]{Figures/TDA_synthetic_o7}&\includegraphics[width=138.76157pt]{Figures/TDA_synthetic_o9}\end{array}$

Figure 8. Synthetic LPPLS signal (top) and corresponding TDA signal (bottom)
for $\omega=5,7,9$.

####  4.1.3. Dependence on non-linearity parameter $m$

We note that when $m\approx 1$, the main trend line of the log price is nearly
linear. The experiments displayed in Fig. 9 show that as the frequency $m$ is
decreased, the TDA signal gets shifted towards the end of the range.

$\begin{array}[]{ccc}\includegraphics[width=138.76157pt]{Figures/synthetic_m07}&\includegraphics[width=138.76157pt]{Figures/synthetic_m05}&\includegraphics[width=138.76157pt]{Figures/synthetic_m01}\\\
\includegraphics[width=138.76157pt]{Figures/TDA_synthetic_m07}&\includegraphics[width=138.76157pt]{Figures/TDA_synthetic_m05}&\includegraphics[width=138.76157pt]{Figures/TDA_synthetic_m01}\end{array}$

Figure 9. Synthetic LPPLS signal (top) and corresponding TDA signal (bottom)
for $m=0.7,0.5,0.1$

####  4.1.4. Dependence of positive/negative bubble parameter $B$

The sign of the parameter $b$ is associated with the bubble-type: $B<0$
corresponds to positive bubble and $b>0$ corresponds to negative bubble. In
Figure 10 we show a synthetic LPPLS model for a negative bubble. It is
remarkable that for both positive and negative bubbles the TDA signal spikes
when approaching the critical time.

$\begin{array}[]{cc}\includegraphics[width=138.76157pt]{Figures/synthetic_b0_0003}&\includegraphics[width=138.76157pt]{Figures/TDA_synthetic_b0_0003}\\\
\end{array}$

Figure 10. Synthetic LPPLS model and norm of persistence landscapes for a
negative bubble $B>0$.

####  4.1.5. Dependence on amplitude parameters $C_{1},C_{2}$

We explore the dependence of the TDA signal on the parameters $C_{1},C_{2}$
which are responsible of the amplitude of the oscillations. When the amplitude
of the signal is increased, the TDA signal grows in strength and becomes more
prominent towards the end of the range. See Fig. 11.

$\begin{array}[]{ccc}\includegraphics[width=138.76157pt]{Figures/synthetic_c1c20_00002}&\includegraphics[width=138.76157pt]{Figures/synthetic_c1c20_00005}&\includegraphics[width=138.76157pt]{Figures/synthetic_c1c20_00007}\\\
\includegraphics[width=138.76157pt]{Figures/TDA_synthetic_c1c20_00002}&\includegraphics[width=138.76157pt]{Figures/TDA_synthetic_c1c20_00005}&\includegraphics[width=138.76157pt]{Figures/TDA_synthetic_c1c20_00007}\end{array}$

Figure 11. Synthetic LPPLS signal (top) and corresponding TDA signal (bottom)
for $C_{1}=C_{2}=0.00002,0.00005,0.00007$.

##  5\. Data segmentation into positive and negative bubbles

In this section we describe a method to segment a given financial time series
into positive and negative bubbles. For a positive bubble the main idea is
there will be a succession of positive returns, which may be interrupted by
negative returns no larger in amplitude than some pre-specified tolerance
level. An analogous description holds for negative bubbles. We recall the data
segmentation algorithm from [GDS18], and some variants [SZ20a, SZ20b,
YRdSTT20].

Let $x(t_{i})$ be the time-series of the price of an asset, where
$t_{i}=t_{i-1}+\Delta t$. The log-return of the asset at time $t_{i}$ is

$p_{i}=\ln x(t_{i})-\ln x(t_{i-1}).$

Assume that some initial time $t_{i_{0}}$ is the beginning of a up- (down-)
trend. The cumulative return at time $i$ is

$p_{i_{0},i}=\sum_{i=i_{0}+1}^{i}(\ln x(t_{i})-\ln x(t_{i-1}))=\ln
x(t_{i})-\ln x(t_{i_{0}}).$

The largest deviation $\delta_{i_{0},i}$ of the price trajectory from a
previous minimum (maximum) is given by

(5.1) $\delta_{i_{0},i}=\left\\{\begin{array}[]{ll}\max_{i_{0}\leq k\leq
i}p_{i_{0},k}-p_{i_{0},i},&\hbox{{ for an upward-trend}},\\\
p_{i_{0},i}-\min_{i_{0}\leq k\leq i}p_{i_{0},k},&\hbox{{ for a downward-
trend}.}\end{array}\right.$

To find the end of the upward-trend, we first compute $i_{\textrm{cross}_{1}}$
the first time $i>i_{0}$ when $\delta_{i_{0},i}$ crosses a certain tolerance
level $\varepsilon_{i}$,

(5.2) $\delta_{i_{0},i}-\varepsilon_{i}>0$

and then define the end of the upward- (downward-) trend as the time

(5.3) $t_{i_{1}}=\left\\{\begin{array}[]{ll}\underset{i_{0}\leq k\leq
i_{\textrm{cross}_{1}}}{\textrm{arg max}}p_{i_{0},k},&\hbox{{ for an upward-
trend}},\\\\[5.0pt] \underset{{i_{0}\leq k\leq
i_{\textrm{cross}_{1}}}}{\textrm{arg min}}p_{i_{0},k},&\hbox{{ for a downward-
trend}.}\end{array}\right.$

The time $t_{i_{1}}$ is considered as the beginning of the next downward
(upward) trend. Through this procedure, the data is successively segmented in
upward and downward trends, which alternate one after the other. The procedure
is carried out for the whole time series and yields a sequence of peaks

$\\{t_{u,1},t_{u,2},\ldots\\},$

and a sequence of troughs

$\\{t_{d,1},t_{d,2},\ldots\\},$

such that peak is always followed by a trough and vice-versa.

At the end of this procedure, each upward-trend corresponds to a positive
bubble, and each downward-trend corresponds to a negative bubble.

The tolerance level $\varepsilon_{i}$ can be chosen in multiple ways. One
option is to choose

$\varepsilon_{i}=\varepsilon,$

where $\varepsilon>0$ is some preset constant. Another option is to choose a
time-dependent tolerance level $\varepsilon_{i}$ proportional to the
volatility (standard deviation) $\sigma_{i}=\sigma_{i}(w_{0})$ at time $i$
estimated over a preceding time window of some fixed size $w_{0}$, i.e.,

(5.4) $\varepsilon_{i}=\varepsilon_{0}\sigma_{i},$

where $\varepsilon_{0}$ is some scaling parameter. Thus the tolerance is more
permissive when the market presents high volatility and is stricter during
calmer periods. There are other ways to choose the tolerance level (see, e.g.,
[YRdSTT20]).

The segmentation procedure in [GDS18] does not fix a single scaling parameter
$\varepsilon_{0}$ and a single window-size $w_{0}$ for the computation of the
standard deviation, but rather considers some range of scaling parameters
$\varepsilon_{0}\in[\varepsilon_{1},\varepsilon_{2}]$ and of window-sizes
$w_{0}\in[w_{1},w_{2}]$, each sampled with some frequency, and computes a
collection of peaks and a collection of troughs for each $\varepsilon_{0}$ and
$w_{0}$. Finally, it selects the most frequent peaks and troughs that appear
over the whole range of $\varepsilon_{0},w_{0}$.

In Fig. 12 we show the log-price of Bitcoin over the period
10/01/2021-07/01/2022 sampled at 1 hour intervals (resulting in 6552 data
points), and the standard deviation measured over a preceding time window of
size $w_{0}=240$.

$\begin{array}[]{cc}\includegraphics[width=138.76157pt]{Figures/Logprice_bitcoin}&\includegraphics[width=138.76157pt]{Figures/Standard_deviation_bitcoin}\end{array}$

Figure 12. Log-price and standard deviation of Bitcoin over
10/01/2021-07/01/2022

In the examples in Section 6.3, as we are interesting in ‘large’ bubbles, we
will fix the scaling parameter at some large value of $\varepsilon_{0}$.

##  6\. Application to detection of Bitcoin bubbles

In this section we present a series of numerical experiments on hourly Bitcoin
price during the period October 1st 2021 and July 1st 2022. In Section 6.1 we
segment the data into positive bubbles and negative bubbles following the
methodology described in Section 5. In Section 6.2 we fit the LPPLS model to
each of the segments. In Section 6.3 we apply the TDA procedure to each
segment, following the methodology described in Sections 2.2 and 2.3.

###  6.1. Data segmentation

The hourly Bitcoin price data set between October 1st 2021 and July 1st 2022
consists of 6552 data points. We segment the data following the procedure
described in Section 5. To compute the standard deviation (volatility)
$\sigma_{i}$ we use a fixed window size of $w_{0}=240$ data points. We also
use a fixed scaling parameter $\varepsilon_{0}=15.0$. Thus, the tolerance
level is taken as $\varepsilon_{i}=\varepsilon_{0}\sigma_{i}$, where
$\sigma_{i}$ is computed over a sliding window of length $w_{0}$.

In the table below we identify the peak-time (P) and the crossing-time (C) for
each upward-trend.

Table 2. Upward-trend peak and crossing times:
P 471 977 2105 3184 3662 3831 4290 4861 5447 5824 6443 C 538 1104 2130 3355
3718 3934 4370 4887 5729 5894 6536

In the table below we identify the trough-time (T) and the crossing-time (C)
for each downward-trend.

Table 3. Downward-trend troughs and crossing times:
T 648 1873 2772 3509 3787 3934 4780 5358 5729 6260 C 783 2010 2821 3531 3831
4075 4835 5601 5804 6415

The combined peaks and troughs from Tables 1 and 2 are shown in Fig. 13.

Figure 13. Peaks and troughs

Then we define upward-trend (U) and downward-trend (D) segments to be fitted
with LPPLS models as well as analyzed with TDA.

First, each segment is determined by the following criterion. For an upward-
(downward-) trend, the end-point is a peak (trough) from the list and the
starting-point of the trend is the nearest crossing time in the past from the
previous downward- (upward-) trend.

Table 4. Upward-trend and downward-trend segments
U 783-977 2010-2105 2821-3184 3531-3831 4075-4290 4835-4861 U 4835-5447
5804-5824 6415-6443 D 538-648 1104-1873 2130-2772 3355-3509 3718-3787
3831-3934 D 4370-4780 4370-5358 4887-5729 5894-6260

Second, we perform some additional adjustments. We examine each of the
resulting segments and if we observe significant peaks or troughs inside a
segment we do an ad-hoc adjustment of the segment: we cut the portion from the
beginning of the segment to the last observed a peak/through. This adjustment
is occasionally needed since using the segmentation procedure with fixed
$w_{0}$ and fixed $\varepsilon_{0}$ can miss some peaks/troughs.

Third, from the remaining segments we dismiss those that are very short, since
applying the LPPLS fitting and the TDA procedure to them is unreliable. Thus,
the segments $4835-4861$, $5804-5824$ and $6415-6443$, all corresponding to
upward-trends, were deemed as too short to perform the analysis. We also
disregard the first segment $1-471$ which was used to initialize the process.

We retain the following:

Table 5. Adjusted upward-trend and downward-trend segments
U 783-977 2010-2105 3020-3184 3531-3831 4175-4290 U 5358-5447 D 538-648
1104-1873 2500-2772 3355-3509 3718-3787 D 3831-3934 4370-4780 5150-5358
4887-5729 5894-6260

###  6.2. Fitting the LPPLS model to Bitcoin time-series

Once a specific segment is defined, we address the optimization problem, in
which we aim at best fitting the LPPLS model to the time series. More
formally, using (6.1) to obtain a vector of estimates of the expected log-
price $\hat{y}=LPPLS(y)$ from the segment $x$, the goal is to solve the
unconstrained, non-convex (Figure 14) optimization problem:

(6.1) $\min_{tc,m,\omega,A,B,C_{1},C_{2}}\quad||\hat{y}-y||^{2}_{2}$

Figure 14. Fitness function in the domain $[m,\omega]$, for constant values of
the remaining parameters. The objective function displays strong non-
convexity.

Given the assumption for the time series to be characterized by Gaussian
noise, the choice of such fitness function is equivalent to maximum likelihood
estimation as depicted in Figure 15.

Figure 15. Distribution of residuals for a fitted model.

For this procedure, together with Gaussianity, we are assuming time-
independent variance; in the case in which heteroscedasticity is taken into
account, a weighting scheme can be introduced in the definition of the fitness
function.

Given the geometry of the problem, a global optimization scheme based on
differential evolution [SP97] has been used. This method, as opposed to
gradient-based local optimizers, is based on performing function evaluation on
a set of sample points drawn from a large candidate space followed by a local
optimizer at the end. While similar works propose the use of grid search in
the space of initial conditions, each followed by local optimization, also
[SZ20b] makes use of evolutionary algorithms (specifically, covariance matrix
adaptation evolution strategy).

In order to speed up convergence, the problem is reformulated as

(6.2)
$\min_{p=[tc,m,\omega]}\quad||\hat{y}(p,A(p),B(p),C_{1}(p),C_{2}(p))-y||^{2}_{2}$

In which $q=[A(p),B(p),C_{1}(p),C_{2}(p)]$ is determined, for each $p$,
solving a linear least square problem:

(6.3) $q=(\textbf{X}^{T}\textbf{X})^{-1}\textbf{X}^{T}\textbf{p}.$

The software developed in this context is built on top of the open-source
package [Bou]. The results of the LPPLS fitting are depicted in Appendix A.

###  6.3. TDA applied to Bitcoin time-series and comparison to fitted LPPLS
models

We apply TDA to each adjusted upward-trend and downward-trend segment from
Table 5. For each segments we have chosen a suitable set of TDA parameters
from the following range:

(6.4) $w=48,60,72,\,d=3,5,\,N=4.$

As noted in Section 4.1, adjustments of the TDA parameter are needed because
the segments have different lengths, and the trends of the corresponding time-
series (reflected by the fitted LPPLS models) are different.

In Appendix A we show, side-by-side, the fitted LPPLS model and TDA signal for
each upward-/downward-trend, in chronological order (i.e., alternating upward-
and downward-trends).

We observe that most of the segments show a primary peak in the TDA signal
towards the end of the range. Some of the segments show some secondary peaks
in the TDA signal (e.g., Appendix A – Figs.17(a), 17(b), 17(d), 18(b), 19(c)),
which seem to capture some intermediate price jumps or drops.

In a few instances (e.g., Appendix A – Fig. 18(e), 19(b))) the primary peak in
the TDA signal occurs earlier than the end of the range. We note that in these
instance that an apparent tipping point seems to occur not at the end of the
range but closer to the middle. The fact that the segmentation procedure used
in Section 5 does not identify this apparent tipping point is likely due to
the dependence of the threshold (5.4) on volatility; if the preceding time
period is highly volatile, the threshold for the crossing time (5.2) is
higher, and if the volatility becomes lower at a later time, then the
threshold will also become lower, and hence the crossing time will be
registered at a later time.

We also note that some of the fitted LPPLS models do not seem to display the
desired characteristics, that is, super-exponential growth (or decay)
superimposed with oscillations increasing in frequency and decreasing in
amplitude. For instance, Appendix A – Figs. 17(c), 18(c), 18(e)) show very few
oscillations, or weak nonlinearity. Such LPPLS models should, in fact, be
dismissed; only LPPLS models that pass certain filters should be selected as
viable solutions (see [GDS18] and the references listed there). Nevertheless,
the TDA signal is clearly distinguishable even in cases when the LPPLS fit is
flawed.

##  7\. Conclusions

We provide an empirical argument on why TDA is able to provide early warning
signals for critical transitions in financial time series, for both positive
and negative financial bubbles. The premise of our investigation is that
financial bubbles are reflected by the LPPLS model, asserting that the price
of an asset exhibits super-exponential growth (or decay) superimposed with
oscillations that increase in frequency and decrease in amplitude when
approaching a critical transition. When the time series is transformed into a
point cloud via time-delay coordinate embedding, these oscillations give rise
to persistent $1$-dimensional homology generators (loops), which can be
quantified via persistence landscapes. By applying a sliding window to the
time series, one can capture the changes in the oscillations through the
growth in the norms of the corresponding persistence landscapes.

Our experiments show that, when the TDA procedure is applied to synthetic
LPPLS time-series, the norms of the persistence landscapes grow as predicted.
They also show an intricate dependence of the TDA signal on the parameters of
the LPPLS model and of the TDA procedure. Further work is needed to fully
analyze this dependence and automatize the parameter selection.

As a practical application, we run the TDA procedure on a collection of
positive and negative bubbles of the Bitcoin price between October 1st 2021
and July 1st 2022. In most cases, the TDA gives early warnings for tipping
points, even when the fitted LPPLS appears to be flawed.

## Appendix A

$\begin{array}[]{cc}\includegraphics[width=225.48424pt]{Figures/graph_538_648}&\includegraphics[width=177.78578pt]{Figures/TDA_538_648}\end{array}$

(a) Data set 538-648

$\begin{array}[]{cc}\includegraphics[width=225.48424pt]{Figures/graph_783_977}&\includegraphics[width=177.78578pt]{Figures/TDA_783_977}\end{array}$

(a) Data set 783-977

$\begin{array}[]{cc}\includegraphics[width=225.48424pt]{Figures/graph_1104_1873}&\includegraphics[width=177.78578pt]{Figures/TDA_1104_1873}\end{array}$

(b) Data set 1104-1873

$\begin{array}[]{cc}\includegraphics[width=225.48424pt]{Figures/graph_2010_2105}&\includegraphics[width=177.78578pt]{Figures/TDA_2010_2105}\end{array}$

(c) Data set 2010-2105

$\begin{array}[]{cc}\includegraphics[width=225.48424pt]{Figures/graph_2500_2772}&\includegraphics[width=177.78578pt]{Figures/TDA_2500_2772}\end{array}$

(d) Data set 2500-2772

$\begin{array}[]{cc}\includegraphics[width=225.48424pt]{Figures/graph_3020_3184}&\includegraphics[width=177.78578pt]{Figures/TDA_3020_3184}\end{array}$

(e) Data set 3020-3184

$\begin{array}[]{cc}\includegraphics[width=225.48424pt]{Figures/graph_3355_3509}&\includegraphics[width=177.78578pt]{Figures/TDA_3355_3509}\end{array}$

(f) Data set 3355-3509

$\begin{array}[]{cc}\includegraphics[width=225.48424pt]{Figures/graph_3718_3787}&\includegraphics[width=177.78578pt]{Figures/TDA_3718_3787}\end{array}$

(a) Data set 3718-3787

$\begin{array}[]{cc}\includegraphics[width=225.48424pt]{Figures/graph_3531_3831}&\includegraphics[width=177.78578pt]{Figures/TDA_3531_3831}\end{array}$

(b) Data set 3531-3831

$\begin{array}[]{cc}\includegraphics[width=225.48424pt]{Figures/graph_3831_3934}&\includegraphics[width=177.78578pt]{Figures/TDA_3831_3934}\end{array}$

(c) Data set 3831-3934

$\begin{array}[]{cc}\includegraphics[width=225.48424pt]{Figures/graph_4175_4290}&\includegraphics[width=177.78578pt]{Figures/TDA_4175_4290}\end{array}$

(d) Data set 4175-4290

$\begin{array}[]{cc}\includegraphics[width=225.48424pt]{Figures/graph_4370_4780}&\includegraphics[width=177.78578pt]{Figures/TDA_4370_4780}\end{array}$

(e) Data set 4370-4780

$\begin{array}[]{cc}\includegraphics[width=225.48424pt]{Figures/graph_5150_5358}&\includegraphics[width=177.78578pt]{Figures/TDA_5150_5358}\end{array}$

(f) Data set 5150-5358

$\begin{array}[]{cc}\includegraphics[width=225.48424pt]{Figures/graph_5358_5447}&\includegraphics[width=177.78578pt]{Figures/TDA_5358_5447}\end{array}$

(a) Data set 5358-5447

$\begin{array}[]{cc}\includegraphics[width=225.48424pt]{Figures/graph_4887_5729}&\includegraphics[width=177.78578pt]{Figures/TDA_4887_5729}\end{array}$

(b) Data set 4887-5729

$\begin{array}[]{cc}\includegraphics[width=225.48424pt]{Figures/graph_5894_6260}&\includegraphics[width=177.78578pt]{Figures/TDA_5894_6260}\end{array}$

(c) Data set 5894-6260

Figure 19. Fitted LPPLS models and norms of persistence landscapes for the
segments from Table 5

## References

  * [ABST93] Henry DI Abarbanel, Reggie Brown, John J Sidorowich, and Lev Sh Tsimring.  The analysis of observed chaotic data in physical systems.  Reviews of modern physics, 65(4):1331, 1993.
  * [ADGK18] Cuneyt Gurcan Akcora, Matthew F Dixon, Yulia R Gel, and Murat Kantarcioglu.  Bitcoin risk modeling with blockchain graphs.  Economics Letters, 173:138–142, 2018.
  * [AKV21] Lloyd L Aromi, Yuri A Katz, and Josep Vives.  Topological features of multivariate distributions: Dependency on the covariance matrix.  Communications in Nonlinear Science and Numerical Simulation, 103:105996, 2021.
  * [ALGK20] Cuneyt G Akcora, Yitao Li, Yulia R Gel, and Murat Kantarcioglu.  Bitcoinheist: Topological data analysis for ransomware prediction on the bitcoin blockchain.  In Proceedings of the Twenty-Ninth International Joint Conference on Artificial Intelligence, 2020.
  * [B+15] Peter Bubenik et al.  Statistical topological data analysis using persistence landscapes.  J. Mach. Learn. Res., 16(1):77–102, 2015.
  * [BCP13] David S Brée, Damien Challet, and Pier Paolo Peirano.  Prediction accuracy and sloppiness of log-periodic functions.  Quantitative Finance, 13(2):275–280, 2013.
  * [BD17] Peter Bubenik and Paweł Dłotko.  A persistence landscapes toolbox for topological statistics.  Journal of Symbolic Computation, 78:91–114, 2017.
  * [BG14] Jesse Berwald and Marian Gidea.  Critical transitions in a model of a genetic regulatory system.  Mathematical Biosciences and Engineering, 11(4):723–740, 2014.
  * [BGVJ14] JJ Berwald, M Gidea, and Mikael Vejdemo-Johansson.  Automatic recognition and tagging of topologically different regimes in dynamical systems.  Discontinuity, Nonlinearity, and Complexity, 3(4):413–426, 2014\.
  * [Bla79] Olivier Jean Blanchard.  Speculative bubbles, crashes and rational expectations.  Economics letters, 3(4):387–389, 1979.
  * [Bou] Boulder-Investment-Technologies, howpublished = https://github.com/boulder-investment-technologies/lppls, note = Accessed: 2023-03-15.
  * [BSP+21] Thomas M Bury, RI Sujith, Induja Pavithran, Marten Scheffer, Timothy M Lenton, Madhur Anand, and Chris T Bauch.  Deep learning for early warning signals of tipping points.  Proceedings of the National Academy of Sciences, 118(39):e2106140118, 2021.
  * [BW82] Olivier J Blanchard and Mark W Watson.  Bubbles, rational expectations and financial markets, 1982.
  * [CDSO14] Frédéric Chazal, Vin De Silva, and Steve Oudot.  Persistence stability for geometric complexes.  Geometriae Dedicata, 173(1):193–214, 2014.
  * [DCB+12] Vasilis Dakos, Stephen R Carpenter, William A Brock, Aaron M Ellison, Vishwesha Guttal, Anthony R Ives, Sonia Kéfi, Valerie Livina, David A Seekell, Egbert H van Nes, et al.  Methods for detecting early warnings of critical transitions in time series illustrated using simulated ecological data.  PloS one, 7(7):e41010, 2012.
  * [DQR19] Paweł Dłotko, Wanling Qiu, and Simon Rudkin.  Cyclicality, periodicity and the topology of time series.  arXiv preprint arXiv:1905.12118, 2019.
  * [DR21] P Dłotko and ST Rudkin.  Financial ratios and stock returns reappraised through a topological data analysis lens.  The European Journal of Finance, pages 1–25, 2021.
  * [DUC20] Meryll Dindin, Yuhei Umeda, and Frederic Chazal.  Topological data analysis for arrhythmia detection through modular neural networks.  In Advances in Artificial Intelligence: 33rd Canadian Conference on Artificial Intelligence, Canadian AI 2020, Ottawa, ON, Canada, May 13–15, 2020, Proceedings 33, pages 177–188. Springer, 2020.
  * [EH22] Herbert Edelsbrunner and John L Harer.  Computational topology: an introduction.  American Mathematical Society, 2022.
  * [FC16] John Fry and Eng-Tuck Cheah.  Negative bubbles and shocks in cryptocurrency markets.  International Review of Financial Analysis, 47:343–352, 2016.
  * [FF96] James A. Feigenbaum and Peter G.O. Freund.  Discrete scale invariance in stock markets before crashes.  International Journal of Modern Physics B, 10(27):3737–3745, 1996\.
  * [GDS18] Jan-Christian Gerlach, Guilherme Demos, and Didier Sornette.  Dissection of Bitcoin’s multiscale bubble history.  Swiss Finance Institute Research Paper, (18-30), 2018.
  * [GGK+20] Marian Gidea, Daniel Goldsmith, Yuri Katz, Pablo Roldan, and Yonah Shmalo.  Topological recognition of critical transitions in time series of cryptocurrencies.  Physica A: Statistical Mechanics and its Applications, 548:123843, 2020.
  * [Gid17] Marian Gidea.  Topological data analysis of critical transitions in financial networks.  In International conference and school on network science, pages 47–59. Springer, 2017.
  * [GK18] Marian Gidea and Yuri Katz.  Topological data analysis of financial time series: Landscapes of crashes.  Physica A: Statistical Mechanics and its Applications, 491:820–834, 2018.
  * [GK23] Victor E Gluzberg and Yuri A Katz.  Topological data analysis of noise: Uniform unimodal distributions.  Communications in Nonlinear Science and Numerical Simulation, page 107216, 2023.
  * [IYG20] Umar Islambekov, Monisha Yuvaraj, and Yulia R Gel.  Harnessing the power of topological data analysis to detect change points.  Environmetrics, 31(1):e2612, 2020.
  * [JLS00] Anders Johansen, Olivier Ledoit, and Didier Sornette.  Crashes as critical points.  International Journal of Theoretical and Applied Finance, 3(02):219–255, 2000.
  * [JSW+96] Anders Johansen, Didier Sornette, Hiroshi Wakita, Urumu Tsunogai, William I Newman, and Hubert Saleur.  Discrete scaling in earthquake precursory phenomena: evidence in the Kobe earthquake, Japan.  Journal de Physique I, 6(10):1391–1402, 1996.
  * [KB21] Yuri A Katz and Alain Biem.  Time-resolved topological data analysis of market instabilities.  Physica A: Statistical Mechanics and its Applications, 571:125816, 2021.
  * [KM18] Firas A Khasawneh and Elizabeth Munch.  Topological data analysis for true step detection in periodic piecewise constant signals.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 474(2218):20180027, 2018.
  * [Len11] Timothy M Lenton.  Early warning of climate tipping points.  Nature climate change, 1(4):201–209, 2011.
  * [LS13] Li Lin and Didier Sornette.  Diagnostics of rational expectation financial bubbles with stochastic mean-reverting termination times.  The European Journal of Finance, 19(5):344–365, 2013.
  * [Per16] Jose A Perea.  Persistent homology of toroidal sliding window embeddings.  In 2016 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 6435–6439. IEEE, 2016.
  * [PH15] Jose A Perea and John Harer.  Sliding windows and persistence: An application of topological methods to signal analysis.  Foundations of Computational Mathematics, 15(3):799–838, 2015.
  * [PY11] Peter CB Phillips and Jun Yu.  Dating the timeline of financial bubbles during the subprime crisis.  Quantitative Economics, 2(3):455–491, 2011.
  * [RCPB19] Rodrigo Rivera-Castro, Polina Pilyugina, and Evgeny Burnaev.  Topological data analysis for portfolio management of cryptocurrencies.  In 2019 International Conference on Data Mining Workshops (ICDMW), pages 238–243. IEEE, 2019.
  * [RQD21] Simon Rudkin, Wanling Qiu, and Pawel Dlotko.  Uncertainty, volatility and the persistence norms of financial time series.  arXiv preprint arXiv:2110.00098, 2021.
  * [SBB+09] Marten Scheffer, Jordi Bascompte, William A Brock, Victor Brovkin, Stephen R Carpenter, Vasilis Dakos, Hermann Held, Egbert H Van Nes, Max Rietkerk, and George Sugihara.  Early-warning signals for critical transitions.  Nature, 461(7260):53–59, 2009.
  * [SH05] K Suchecki and JA Holyst.  Log-periodic oscillations in degree distributions of hierarchical scale-free networks.  Acta Physica Polonica. Series B, 35, 2005.
  * [SNC20] P. Saengduean, S. Noisagool, and F. Chamchod.  Topological data analysis for identifying critical transitions in cryptocurrency time series.  In 2020 IEEE International Conference on Industrial Engineering and Engineering Management (IEEM), pages 933–938, 2020.
  * [Sor17] Didier Sornette.  Why stock markets crash.  In Why Stock Markets Crash. Princeton university press, 2017.
  * [SP97] Rainer Storn and Kenneth Price.  Differential evolution-a simple and efficient heuristic for global optimization over continuous spaces.  Journal of global optimization, 11(4):341, 1997.
  * [SS16] Maximilian Seyrich and Didier Sornette.  Micro-foundation using percolation theory of the finite time singular behavior of the crash hazard rate in a class of rational expectation bubbles.  International Journal of Modern Physics C, 27(10):1650113, 2016\.
  * [SS17] H Sohn and Didier Sornette.  Bubbles as violations of efficient time-scales.  Institute for New Economic Thinking Working Paper Series, (65), 2017\.
  * [SSGC+18] Manish Saggar, Olaf Sporns, Javier Gonzalez-Castillo, Peter A Bandettini, Gunnar Carlsson, Gary Glover, and Allan L Reiss.  Towards a new approach to reveal dynamical organization of the brain using topological data analysis.  Nature communications, 9(1):1399, 2018.
  * [SSZ21] Min Shu, Ruiqiang Song, and Wei Zhu.  The ‘covid’crash of the 2020 us stock market.  The North American Journal of Economics and Finance, 58:101497, 2021\.
  * [SSZ22] Ruiqiang Song, Min Shu, and Wei Zhu.  The 2020 global stock market crash: Endogenous or exogenous?  Physica A: Statistical Mechanics and Its Applications, 585:126425, 2022.
  * [SZ20a] Min Shu and Wei Zhu.  Detection of Chinese stock market bubbles with LPPLS confidence indicator.  Physica A: Statistical Mechanics and Its Applications, 557:124892, 2020.
  * [SZ20b] Min Shu and Wei Zhu.  Real-time prediction of Bitcoin bubble crashes.  Physica A: Statistical Mechanics and its Applications, 548:124477, 2020.
  * [Tak81] Floris Takens.  Detecting strange attractors in turbulence.  In David Rand and Lai-Sang Young, editors, Dynamical Systems and Turbulence, Warwick 1980, pages 366–381, Berlin, Heidelberg, 1981. Springer Berlin Heidelberg.
  * [TJM91] Sauer Tim, A Yorke James, and Casdagli Martin.  Embedology.  Journal of Statistical Physics, 65(3-4):579–616, 1991.
  * [TS11] J Michael T Thompson and Jan Sieber.  Climate tipping as a noisy bifurcation: a predictive technique.  IMA Journal of Applied Mathematics, 76(1):27–46, 2011.
  * [YRdSTT20] Arthur Matsuo Yamashita Rios de Sousa, Hideki Takayasu, and Misako Takayasu.  Segmentation of time series in up-and down-trends using the epsilon-tau procedure with application to USD/JPY foreign exchange market data.  PloS one, 15(9):e0239494, 2020.
  * [YRWS12] Wanfeng Yan, Reda Rebib, Ryan Woodard, and Didier Sornette.  Detection of crashes and rebounds in major equity markets.  International Journal of Portfolio Analysis and Management, 1(1):59–79, 2012.
  * [YWS12] Wanfeng Yan, Ryan Woodard, and Didier Sornette.  Diagnosis and prediction of rebounds in financial markets.  Physica A: Statistical Mechanics and its Applications, 391(4):1361–1380, 2012.


# s- and p-superfluidity of Fermi atoms in Bose-Fermi mixtures

E.V. Gorbar1,2, Y.O. Nikolaieva1, I.V. Oleinikova3, S.I. Vilchinskii1, A.I.
Yakimenko1 1 Department of Physics, Taras Shevchenko National University of
Kyiv, 64/13, Volodymyrska Street, Kyiv 01601, Ukraine
2 Bogolyubov Institute for Theoretical Physics, 14-b, Metrologichna, Kyiv
03143, Ukraine
3 National University of Technologies and Design, 2, Nemirovich-Danchenko
Street, Kyiv 01011, Ukraine

###### Abstract

The $p$-wave superfluid is characterized by nontrivial topological
characteristics essential for fault-tolerant quantum state manipulation.
However, the practical realization of the $p$-wave state remains a challenging
problem. We study the $s$\- and $p$-wave superfluidity in mixtures of
fermionic and spinor bosonic gases and derive a general set of the gap
equations for these superfluid states. Numerically solving the gap equations
for the $s$-wave state, we quantify the physical conditions for the
realization of the pure $p$-wave state in a well-controlled environment of
atomic physics in the absence of an admixture of the $s$-wave state.

##  I Introduction

Chiral $p$-wave superfluids and superconductors have attracted much attention
because of their rich physics and remarkable potential applications. Vortex
excitations in a two-dimensional version of such a condensate are anyonic and
obey non-Abelian statistics. Topologically degenerate quasiparticle states of
chiral superconductors can be used for the creation of fault-tolerant
topological quantum computers Nayak et al. (2008). Cold atomic gases provide a
promising physical platform for the realization of topologically non-trivial
states. The $p$-wave fermionic superfluidity in cold gases has become the
topic of extensive theoretical research Inotani et al. (2016); Wu and Bruun
(2016); Cooper and Shlyapnikov (2009); Zhang et al. (2008); Juliá-Díaz et al.
(2013); Han et al. (2009); Bardyn et al. (2012); Thompson et al. (2020); Sato
et al. (2009). Although there are various $p$-wave Fermi superfluids such as
superfluid phases of ${}^{3}\mbox{He}$ Leggett (1975) and heavy-fermion
superconductors Sigrist and Ueda (1991); Stewart (1984); MacKenzie and Maeno
(2003), the realization of $p$-wave superfluid state in a Fermi gas with
tunable pairing interaction would strongly advance the physics of
unconventional Fermi superfluids both theoretically and experimentally.

The properties of fermionic pairing are determined by the agent, which
mediates the effective attraction between fermions. In quantum gases, a Bose-
Einstein condensate (BEC) can serve as such a mediator in Bose-Fermi mixtures.
Recent theoretical studies reveal rich physical properties of the fermi
superfluidity in Bose-Fermi mixtures Wu et al. (2020); Bazak and Petrov
(2018); Wang et al. (2018). In Ref. Efremov and Viverit (2002), it was shown
that the $p$-wave Cooper pairing of spin-polarized fermions (this makes
impossible the $s$-wave pairing) takes place in a binary fermion-boson mixture
due to the exchange of density fluctuations of the bosonic medium. The main
problem is that the $s$-wave pairing (with zero orbital angular momentum of
the pair) has typically a higher critical temperature than the $p$-wave
pairing and, therefore, is more favorable. As was demonstrated experimentally
Regal et al. (2003) in a single-component ultracold Fermi gas of 40K atoms, by
using the Feshbach resonance, it is possible to essentially enhance the
$p$-wave collision cross-section to values larger than the background $s$-wave
cross-section between atoms in different spin states. However, the Feshbach
resonance also significantly enhances three-body inelastic collisional loss
Levinsen et al. (2007, 2008); Jona-Lasinio et al. (2008).

A spinor Bose-Einstein condensate (BEC) supports various types of excitations
not available in single-component BECs, which offers considerable scope for
the realization of unconventional superfluid states of quantum gases. In our
recent paper Matsyshyn et al. (2018), we analyzed the $p$-wave polar phase in
mixtures of ultracold Fermi and $F=1$ Bose gases. It was found that the spin-
induced interactions in a mixture of spinor atomic BEC and fermionic atoms
naturally suppress the $s$-wave pairing. At the same time, the $p$-wave
pairing is driven both by density-induced and spin-induced interactions. Thus,
in contrast to the $s$-wave state, the $p$-wave one emerges even for repulsive
interaction between fermions. Although we demonstrated in Ref. Matsyshyn et
al. (2018) that the polar $p$-wave phase can be realized, its competition with
the $s$-wave state was left as an open question.

In the present work, we study the competition between the $s$\- and $p$-wave
superfluidity of the fermi atoms. We found numerical solutions to the gap
equation in $s-$wave state and quantify the physical conditions for the
realization of the pure $p$-wave state in the absence of an admixture of the
$s$-wave state.

The rest of the paper is organized as follows. The model is introduced in Sec.
II. In Sec. III, we derive the general gap equation and present the polar
$p$-wave and $s$-wave gap equations. In Sec. IV, we analyze the $s$-wave state
and present our numerical results of the corresponding gap equation. The paper
is concluded in Sec. V.

##  II Model

Our model of a mixture of electrically neutral $F=1$ bosonic and
spin-$\frac{1}{2}$ fermionic atoms is defined by the following second
quantized Hamiltonian of dilute Fermi and spinor Bose gases:

$\hat{H}=\int\left[\mathcal{\hat{H}}_{B}(\textbf{r})+\mathcal{\hat{H}}_{F}(\textbf{r})+\mathcal{\hat{H}}_{BF}(\textbf{r})\right]d\textbf{r}.$
(1)

The first term $\mathcal{\hat{H}}_{B}(\textbf{r})$ in the Hamilton density
describes interacting $F=1$ spinor bosonic atoms with no external trapping
potential Ho (1998); Kawaguchi and Ueda (2012)

$\displaystyle\mathcal{\hat{H}}_{B}(\textbf{r})=\sum_{m}\hat{\Psi}_{m}^{+}(\textbf{r})\left(-\frac{\hbar^{2}}{2M_{B}}\nabla^{2}\right)\hat{\Psi}_{m}(\textbf{r})+\frac{1}{2}c_{0}:\hat{n}_{B}^{2}(\textbf{r}):+\frac{1}{2}c_{2}:\textbf{F}^{2}:,$

where $:\hat{A}:$ denotes normally ordered operator $\hat{A}$ and the field
operators $\hat{\Psi}(\textbf{r})$ obey the standard commutation relations

$\left[\hat{\Psi}_{m}(\textbf{r}),\hat{\Psi}_{m^{\prime}}^{+}(\textbf{r})\right]=\delta_{m,m^{\prime}}\delta(\textbf{r}-\textbf{r}^{\prime}).$

Further,
$c_{0}=\frac{4\pi\hbar^{2}}{3M_{B}}\left(a^{(0)}_{B}+2a^{(2)}_{B}\right)$ and
$c_{2}=\frac{8\pi\hbar^{2}}{3M_{B}}\left(a^{(2)}_{B}-a^{(0)}_{B}\right)$ are
symmetric and spin-dependent interaction constants, respectively, where
$a^{(J)}_{B}$ is the scattering length in the state with spin $J$. The
operator
$\hat{n}_{B}(\textbf{r})=\sum_{m}\hat{\Psi}_{m}^{+}(\textbf{r})\hat{\Psi}_{m}(\textbf{r})$,
$m\in\\{-1,0,+1\\}$ denotes the number density for bosonic fields. As to the
ultracold gas of spin-1 bosonic atoms, we consider in this paper 23Na atoms in
the polar spinor Bose–Einstein condensate (BEC) state with
$c_{2}/c_{0}=6.25\cdot 10^{-2}$ and $a^{(0)}_{B}=50$, $a^{(2)}_{B}=55$ in
units of the Bohr radius $a_{B}$.

The second term $\mathcal{\hat{H}}_{F}(\textbf{r})$ in the Hamilton density
(1) describes the gas of Fermi atoms without external trapping potential

$\mathcal{\hat{H}}_{F}(\textbf{r})=\sum_{\nu}\hat{f}_{\nu}^{+}(\textbf{r})\left(-\frac{\hbar^{2}}{2M_{F}}\nabla^{2}-\mu\right)\hat{f}_{\nu}(\textbf{r})+\frac{1}{2}g_{F}:\hat{n}_{F}^{2}(\textbf{r}):$
(2)

with the field operators $\hat{f}$ obeying the standard anticommutation
relations

$\left\\{\hat{f}_{\nu}(\textbf{r}),\hat{f}_{\nu^{\prime}}^{+}(\textbf{r})\right\\}=\delta_{\nu,\nu^{\prime}}\delta(\textbf{r}-\textbf{r}^{\prime})$

and $\mu$ is the fermion chemical potential. The operator
$\hat{n}_{F}(\textbf{r})=\sum_{\nu}\hat{f}_{\nu}^{\dagger}(\textbf{r})\hat{f}_{\nu}(\textbf{r})$,
$\nu\in\\{\uparrow,\downarrow\\}$ denotes the number density of Fermi atoms,
$g_{F}=4\pi\hbar^{2}a_{F}/M_{F}$ is the fermion-fermion coupling constant,
where $a_{F}$ is the $s$-wave scattering length.

The last term $\mathcal{\hat{H}}_{BF}(\textbf{r})$ in the Hamilton density (1)
describes the interaction between the Bose and Fermi atoms

$\mathcal{\hat{H}}_{BF}(\textbf{r})=\frac{\alpha}{2}\,\hat{n}_{B}(\textbf{r})\hat{n}_{F}(\textbf{r})+\frac{\beta}{2}\,\textbf{S}\cdot\textbf{F},$
(3)

where
$\alpha=\frac{4\pi\hbar^{2}}{3M_{BF}}\left(2a^{(1/2)}_{BF}+a^{(3/2)}_{BF}\right)$
and
$\beta=\frac{8\pi\hbar^{2}}{3M_{BF}}\left(a^{(1/2)}_{BF}-a^{(3/2)}_{BF}\right)$
are coupling constants of the density-density and spin-spin interactions of
Bose and Fermi atoms, respectively. Further, $a^{(J)}_{BF}$ is the scattering
length in the state with spin $J$ and $M_{BF}=2M_{B}M_{F}/(M_{B}+M_{F})$ is
the reduced mass for a boson with mass $M_{B}$ and a fermion with mass
$M_{F}$.

The $j$-th component of the vector spin density operator for bosons F equals

$\textbf{F}_{j}=\sum_{m,m^{\prime}}\hat{\Psi}_{m}^{+}(\textbf{r})\left\\{\hat{F}_{j}\right\\}_{m,m^{\prime}}\hat{\Psi}_{m^{\prime}}(\textbf{r}),$

where $\hat{F}_{j}$ are the $3\times 3$ matrices of the vector representation
of the rotation group. The spin density operator for fermions has a similar
form

$\textbf{S}_{j}=\sum_{\nu,\nu^{\prime}}\hat{f}_{\nu}^{\dagger}(\textbf{r})\left\\{\hat{\sigma}_{j}\right\\}_{\nu,\nu^{\prime}}\hat{f}_{\nu^{\prime}}(\textbf{r}),$

where $\hat{\sigma}_{j}$ are the Pauli-matrices for spin-$\frac{1}{2}$
fermions.

To be specific, we fix in our study the fermionic and bosonic densities
$n_{F}=n_{B}=10^{14}$ cm-3 for 40K and 23Na atoms that gives the Fermi energy
$E_{F}=(3\pi^{2}n_{F})^{2/3}\hbar^{2}/(2M_{F})=\hbar^{2}k_{F}^{2}/(2M_{F})$.

##  III Gap equations

To derive the gap equations for both $s$\- and $p$-wave states, it is
convenient to obtain the effective action for Fermi atoms. Using Hamiltonian
(1), it is straightforward to show that the fermion-boson interactions lead to
the following effective fermion’s action:

$S_{\text{eff}}=S_{0}+S_{int},$ (4)

where

$S_{0}=\int
d\textbf{r}dt\Big{[}\sum_{\nu}\hat{f}^{+}_{\nu}\left(i\hbar\frac{\partial}{\partial
t}+\frac{\hbar^{2}}{2m}\nabla^{2}+\mu\right)\hat{f}^{+}_{\nu}-\frac{1}{2}g_{F}:\hat{n}_{F}^{2}(\textbf{r}):\Big{]}.$
(5)

is the free action of Fermi atoms and the interaction term is given by

$S_{int}=-\sum_{\textbf{k}\neq 0}\frac{\epsilon_{\textbf{k}}n_{B}}{2}\int
d\omega\Big{\\{}\frac{\beta^{2}}{E_{\textbf{k},1}^{2}-\hbar^{2}\omega^{2}}[\sigma_{+}(k)\sigma_{-}(-k)+\sigma_{-}(k)\sigma_{+}(-k)]+\frac{\alpha^{2}}{2(E_{\textbf{k},0}^{2}-\hbar^{2}\omega^{2})}n_{F}(k)n_{F}(-k)\Big{\\}}.\\\
$ (6)

Here $\epsilon_{\textbf{k}}=\hbar^{2}\textbf{k}^{2}/(2M_{B})$,
$k=(\omega,\mathbf{k})$, $r=(t,\mathbf{r})$,

$\sigma_{\pm}=\frac{1}{2\sqrt{\Omega}}\int
e^{i\textbf{kr}}\hat{f}^{+}(\textbf{r})(\hat{\sigma}_{x}\pm
i\hat{\sigma}_{y})\hat{f}(\textbf{r})d\textbf{r},\\\ $ (7)
$E_{\bm{k},0}=\sqrt{\epsilon_{\bm{k}}(\epsilon_{\bm{k}}+2c_{0}n_{B})},\quad
E_{\bm{k},1}=\sqrt{\epsilon_{\bm{k}}(\epsilon_{\bm{k}}+2c_{2}n_{B})},$ (8)

and $\Omega$ is the system’s volume.

The effective action (4) implies the following Dyson’s equations:

$(i\hbar\omega_{n}-\xi_{\textbf{k}})G(\textbf{k},\omega_{n})+\Delta(\textbf{k},\omega_{n})F^{+}(\textbf{k},\omega_{n})=\hbar,$
(9)
$(i\hbar\omega_{n}+\xi_{\textbf{k}})F^{+}(\textbf{k},\omega_{n})+\Delta^{+}(\textbf{k},\omega_{n})G(\textbf{k},\omega_{n})=0,$
(10)

for the fermion Green function $G(\textbf{k},\omega_{n})$ and the anomalous
Gor‘kov function $F(\textbf{k},\omega_{n})$ in the Matsubara formalism, where
$\xi_{\mathbf{k}}=\hbar^{2}\textbf{k}^{2}/(2M_{F})-\mu$. Here we introduced
the order parameter of general type

$\Delta=-\frac{k_{\textrm{B}}T}{\hbar}\sum_{\nu,\textbf{k}^{\prime}\neq
0}\Big{[}\frac{g_{F}}{\Omega}F(\textbf{k}-\textbf{k}^{\prime},\omega_{\nu})-\tilde{g}_{n}(\textbf{k}^{\prime},\omega_{\nu})F(\textbf{k}-\textbf{k}^{\prime},\omega_{n}-\omega_{\nu})-\tilde{g}_{s}(\textbf{k}^{\prime},\omega_{\nu})\sum_{s=+,-}\hat{\sigma}_{s}F(\textbf{k}-\textbf{k}^{\prime},\omega_{n}-\omega_{\nu})\hat{\sigma}_{s}\Big{]},$
(11)

where

$\displaystyle\tilde{g}_{n}(\textbf{k},\omega_{\nu})=\frac{c_{n}}{E_{\textbf{k},0}^{2}+\hbar^{2}\omega_{\nu}^{2}},$
$\displaystyle\qquad c_{n}=\frac{\alpha^{2}n\epsilon_{\textbf{k}}}{2\Omega},$
$\displaystyle\tilde{g}_{s}(\textbf{k},\omega_{\nu})=\frac{c_{s}}{E_{\textbf{k},1}^{2}+\hbar^{2}\omega_{\nu}^{2}},$
$\displaystyle\qquad c_{s}=\frac{\beta^{2}n\epsilon_{\textbf{k}}}{\Omega}.$

Since we wish to describe both $s$\- and $p$\- wave states, we consider the
following ansatz for the energy gap:

$\Delta_{\alpha\beta}(\textbf{k})=\Delta_{p}(i\bm{\sigma}\sigma_{y})_{\alpha\beta}\textbf{d}(\textbf{k})+\Delta_{s}(i\sigma_{y})_{\alpha\beta}.$
(12)

The vector d determines the form and symmetries of the p-wave superfluid
state. Further, $\alpha,\beta$ are the spin indices, $\bm{\sigma}$ are the
Pauli matrices, and $\Delta_{s}$ and $\Delta_{p}$ are scalar functions of
$\omega_{n}$ and $|\mathbf{k}|$. Using this ansatz, Eqs. (9) and (10) give

$G_{\alpha\beta}(\textbf{k},\omega_{n})=-(i\hbar\omega_{n}+\xi_{\textbf{k}})\cdot\frac{(\hbar^{2}\omega_{n}^{2}+\xi_{\textbf{k}}^{2}+\Delta^{2})\delta_{\alpha\beta}-\Delta^{2}\textbf{M}\bm{\sigma}_{\alpha\beta}}{(\hbar^{2}\omega_{n}^{2}+E_{+}^{2})(\hbar^{2}\omega_{n}^{2}+E_{-}^{2})},$
(13) $\displaystyle F_{\alpha\beta}(\textbf{k},\omega_{n})$ $\displaystyle=$
$\displaystyle\hbar\Delta_{p}(\bm{\sigma}i\sigma_{y})_{\alpha\beta}\cdot\frac{(\hbar^{2}\omega_{n}^{2}+\xi_{\textbf{k}}^{2}+\Delta^{2})\textbf{d}-i\Delta^{2}[\textbf{M}\cdot\textbf{d}]}{(\hbar^{2}\omega_{n}^{2}+E_{+}^{2})(\hbar^{2}\omega_{n}^{2}+E_{-}^{2})}$
(14)
$\displaystyle+\hbar\Delta_{s}\frac{(\hbar^{2}\omega_{n}^{2}+\xi_{\textbf{k}}^{2}+\Delta^{2}-2\Delta_{p}^{2}(Re[\textbf{d}]\cdot\textbf{d}))i\sigma_{y,\alpha\beta}}{(\hbar^{2}\omega_{n}^{2}+E_{+}^{2})(\hbar^{2}\omega_{n}^{2}+E_{-}^{2})}+\hbar\Delta_{s}\frac{\Delta^{2}(\bm{\sigma}i\sigma_{y})_{\alpha\beta}\textbf{M}}{(\hbar^{2}\omega_{n}^{2}+E_{+}^{2})(\hbar^{2}\omega_{n}^{2}+E_{-}^{2})},$

where

$\displaystyle\Delta^{2}=\Delta_{p}^{2}|\textbf{d}|^{2}+\Delta_{s}^{2},$
$\displaystyle\textbf{m}(\textbf{k})=i[\textbf{d}(\textbf{k})\times\textbf{d}^{*}(\textbf{k})],$
$\displaystyle\textbf{M}=\frac{\Delta_{p}^{2}\textbf{m}+2\Delta_{p}\Delta_{s}Re[\textbf{d}]}{\Delta^{2}},$
$\displaystyle E_{\pm}=\sqrt{\xi_{\textbf{k}}^{2}+\Delta^{2}(1\pm|M|)}.$

For the polar $p$-wave state with $\mathbf{d}=(0,0,k_{z})$, Eqs.(11) and (14)
lead to the following $p$-wave gap equation:

$k_{z}\Delta_{p}(\textbf{k},\omega_{k})=k_{\textrm{B}}T\,\sum_{n,\textbf{k}^{\prime}\neq
0}(\tilde{g}_{n}(\textbf{k}^{\prime},\omega_{n})+\tilde{g}_{s}(\textbf{k}^{\prime},\omega_{n}))\,\,\frac{(k_{z}-k^{\prime}_{z})\Delta_{p}(\textbf{k}-\textbf{k}^{\prime},\omega_{k}-\omega_{n})}{\hbar^{2}(\omega_{k}-\omega_{n})^{2}+\xi^{2}_{\textbf{k}-\textbf{k}^{\prime}}+(k_{z}-k^{\prime}_{z})^{2}\Delta^{2}_{p}}.$
(15)

The same equations for the s-wave gap $\Delta_{s}$ give the following gap
equation:

$\Delta_{s}(\textbf{k},\omega_{k})=k_{\textrm{B}}T\,\sum_{n,\textbf{k}^{\prime}\neq
0}\left(-\frac{g_{F}}{\Omega}+\tilde{g}_{n}(\textbf{k}^{\prime},\omega_{n})-\tilde{g}_{s}(\textbf{k}^{\prime},\omega_{n})\right)\,\,\frac{\Delta_{s}(\textbf{k}-\textbf{k}^{\prime},\omega_{k}-\omega_{n})}{\hbar^{2}(\omega_{k}-\omega_{n})^{2}+\xi^{2}_{\textbf{k}-\textbf{k}^{\prime}}+\Delta^{2}_{s}}.$
(16)

###  III.1 p-wave gap equation

Although we concentrate in this paper on the $s$-wave state, it is useful to
present here the details of the derivation of the $p$-wave gap equation
obtained in Ref. Matsyshyn et al. (2018). Setting $\omega_{k}=\omega_{0}=\pi
k_{\textrm{B}}T/\hbar$ in Eq.(15) and neglecting the dependence of the gap on
the Matsubara frequency, we obtain the following equation for the polar phase
gap $\Delta(\bm{k})=k_{z}\Delta_{p}(\mathbf{k},\omega_{0})$:

$\Delta(\bm{k})=k_{\textrm{B}}T\cdot\sum_{n,\textbf{k}^{\prime}\neq
0}(\tilde{g}_{n}(\textbf{k}^{\prime},\omega_{n}-\omega_{0})+\tilde{g}_{s}(\textbf{k}^{\prime},\omega_{n}-\omega_{0}))\cdot\frac{\Delta(\textbf{k}-\textbf{k}^{\prime})}{\hbar^{2}\omega_{n}^{2}+E^{2}(\textbf{k}-\textbf{k}^{\prime})}.$
(17)

By using the definition of $\tilde{g}_{n}$ and $\tilde{g}_{s}$ below Eq.(11),
we conclude that the summation over $n$ boils down to the calculation of the
following sum with certain $a$ and $b$:

$\displaystyle
S=k_{\textrm{B}}T\sum_{n}\frac{1}{\pi^{2}k_{\textrm{B}}^{2}T^{2}(2n+1)^{2}+a^{2}}\cdot\frac{1}{\pi^{2}k_{\textrm{B}}^{2}T^{2}(2n+1-k)^{2}+b^{2}}$
$\displaystyle=k_{\textrm{B}}T\sum_{n,s1,s2}\frac{s_{1}s_{2}}{4ab}\frac{1}{i\pi
k_{\textrm{B}}T(2n+1)-s_{1}a}\cdot\frac{1}{i\pi k_{\textrm{B}}T(2n+1)-i\pi
k_{\textrm{B}}Tk-s_{2}b},$ (18)

which can be explicitly calculated, where we introduced auxiliary variables
$s_{1},s_{2}$ with values $\pm 1$. Indeed, we have

$k_{\textrm{B}}T\sum\frac{1}{(i\hbar\omega_{n}-\xi_{1})(i\hbar\omega_{n}-\xi_{2})}=\frac{n_{F}(\xi_{1})-n_{F}(\xi_{2})}{\xi_{1}-\xi_{2}},$
(19)

where $\xi_{1}=s_{1}a$, and $\xi_{2}=s_{2}b+i\pi k_{\textrm{B}}Tk$. Thus, we
obtain

$S=\sum_{s_{1},s_{2}}\frac{s_{1}s_{2}}{4ab}\frac{n_{F}(s_{1}a)+n_{B}(s_{2}b)}{s_{1}a-s_{2}b-i\hbar\omega},$
(20)

where we used

$n_{F}(i\pi k_{\textrm{B}}Tk+x)=-n_{B}(x),\qquad k=2q+1.$ (21)

Replacing $s_{2}\to-s_{2}$, we have

$S=\sum_{s1,s2}-\frac{s_{1}s_{2}}{4ab}\frac{-n_{B}(-s_{2}b)-n_{F}(s_{1}a)}{i\hbar\omega-
s_{1}a-s_{2}b}=-\sum_{s1,s2}\frac{s_{1}s_{2}}{4ab}\frac{1+n_{B}(s_{2}b)-n_{F}(s_{1}a)}{i\hbar\omega-
s_{1}a-s_{2}b}.$ (22)

Using

$n_{F}(x)=\frac{1-\tanh\frac{x}{2k_{\textrm{B}}T}}{2},\qquad
n_{B}(x)=\frac{\coth\frac{x}{2k_{\textrm{B}}T}-1}{2},$ (23)

we can rewrite

$S=\sum_{s_{1},s_{2}}\frac{-1}{8ab}\frac{s_{1}\coth\frac{b}{2k_{\textrm{B}}T}+s_{2}\tanh\frac{a}{2k_{\textrm{B}}T}}{i\hbar\omega-
s_{1}a-s_{2}b}.$ (24)

Let us calculate first

$\displaystyle\sum_{s1,s2}\frac{-1}{8ab}\frac{(i\hbar\omega+s_{1}a-s_{2}b)s_{1}\coth\frac{b}{2k_{\textrm{B}}T}}{(i\hbar\omega-
s_{2}b)^{2}-a^{2}}=\sum_{s_{2}}\frac{-1}{4b}\frac{\coth\frac{b}{2k_{\textrm{B}}T}}{-a^{2}-\hbar^{2}\omega^{2}+b^{2}-2i\hbar\omega
s_{2}b}$
$\displaystyle=\frac{\coth\frac{b}{2k_{\textrm{B}}T}}{2b}\frac{\hbar^{2}\omega^{2}+a^{2}-b^{2}}{(\hbar^{2}\omega^{2}+a^{2}-b^{2})^{2}+4\hbar^{2}\omega^{2}b^{2}}.$
(25)

Similarly we have

$S=\frac{\coth\frac{b}{2k_{\textrm{B}}T}}{2b}\frac{\hbar^{2}\omega^{2}+a^{2}-b^{2}}{(\hbar^{2}\omega^{2}+a^{2}-b^{2})^{2}+4\hbar^{2}\omega^{2}b^{2}}+\frac{\tanh\frac{a}{2k_{\textrm{B}}T}}{2a}\frac{\hbar^{2}\omega^{2}+b^{2}-a^{2}}{(\hbar^{2}\omega^{2}+b^{2}-a^{2})^{2}+4\hbar^{2}\omega^{2}a^{2}}.$
(26)

According to Eq.(17),

$b=E_{k-k^{\prime}},\qquad a=E_{\mathbf{k},i},\qquad i=\\{0,1\\},$ (27)

therefore, we have

$S=\frac{Z_{i}}{2E_{\mathbf{k},i}}\tanh\frac{E_{\mathbf{k},i}}{2k_{\textrm{B}}T}+\frac{W_{i}}{2E_{k-k^{\prime}}}\coth\frac{E_{k-k^{\prime}}}{2k_{\textrm{B}}T},$
(28)

where

$\displaystyle
Z_{i}=\frac{\hbar^{2}\omega^{2}+E_{k-k^{\prime}}^{2}-E_{\mathbf{k},i}^{2}}{(E_{k-k^{\prime}}^{2}-E_{\mathbf{k},i}^{2})^{2}+2\hbar^{2}\omega^{2}\left(E_{\mathbf{k},i}^{2}+E_{k-k^{\prime}}^{2}+\frac{\hbar^{2}\omega^{2}}{2}\right)},$
(29) $\displaystyle
W_{i}=\frac{\hbar^{2}\omega^{2}+E_{\mathbf{k},i}^{2}-E_{k-k^{\prime}}^{2}}{(E_{k-k^{\prime}}^{2}-E_{\mathbf{k},i}^{2})^{2}+2\hbar^{2}\omega^{2}\left(E_{\mathbf{k},i}^{2}+E_{k-k^{\prime}}^{2}+\frac{\hbar^{2}\omega^{2}}{2}\right)}.$
(30)

Finally, we obtain

$\Delta(\bm{k})=\sum_{\bm{k^{\prime}}\neq
0}\Big{\\{}c_{n}\left(\frac{W_{0}}{2E_{k-k^{\prime}}}\coth{\frac{E_{k-k^{\prime}}}{2k_{\textrm{B}}T}}+\frac{Z_{0}}{2E_{\bm{k},0}}\tanh{\frac{E_{\bm{k},0}}{2k_{\textrm{B}}T}}\right)+c_{s}\left(\frac{W_{1}}{2E_{k-k^{\prime}}}\coth{\frac{E_{k-k^{\prime}}}{2k_{\textrm{B}}T}}+\frac{Z_{1}}{2E_{\bm{k},1}}\tanh{\frac{E_{\bm{k},1}}{2k_{\textrm{B}}T}}\right)\Big{\\}}\cdot\Delta(\bm{k}-\bm{k^{\prime}}).$
(31)

In what follows, we study the $s$-wave phase in detail because the $p$-wave
polar state was previously investigated in our recent work Matsyshyn et al.
(2018) including numerical solution of the corresponding gap equation.

###  III.2 s-wave gap equation

Similarly setting $\omega_{k}=\omega_{0}=\pi k_{\textrm{B}}T/\hbar$ in Eq.(16)
and neglecting the dependence of the gap on the Matsubara frequency, we obtain
the following equation for the $s$-wave gap:

$\Delta(\bm{k})=\sum_{\bm{k^{\prime}}\neq
0}\Big{\\{}-\frac{g_{F}}{\Omega}\frac{\tanh\frac{E_{k-k^{\prime}}}{2k_{\textrm{B}}T}}{2E_{k-k^{\prime}}}+c_{n}\left(\frac{W_{0}}{2E_{k-k^{\prime}}}\coth{\frac{E_{k-k^{\prime}}}{2k_{\textrm{B}}T}}+\frac{Z_{0}}{2E_{\bm{k},0}}\tanh{\frac{E_{\bm{k},0}}{2k_{\textrm{B}}T}}\right)\\\
-c_{s}\left(\frac{W_{1}}{2E_{k-k^{\prime}}}\coth{\frac{E_{k-k^{\prime}}}{2k_{\textrm{B}}T}}+\frac{Z_{1}}{2E_{\bm{k},1}}\tanh{\frac{E_{\bm{k},1}}{2k_{\textrm{B}}T}}\right)\Big{\\}}\cdot\Delta(\bm{k}-\bm{k^{\prime}}),$
(32)

where we used also the Poisson summation formula with the subsequent extension
of the contour to infinity)

$\sum^{+\infty}_{n=-\infty}\frac{1}{k_{\textrm{B}}^{2}T^{2}\pi^{2}(2n+1)^{2}+a^{2}}=\frac{\tanh(\frac{a}{2k_{\textrm{B}}T})}{2k_{\textrm{B}}Ta}$

and $W_{i}$, $Z_{i}$ are defined in Eqs.(29) and (30). Now let us find
solution to the gap equation of the $s$-wave state.

##  IV Results for the s-wave state

Before we find the solution to the gap equation for the $s$-wave state
Eq.(32), it is useful to check the correctness of the obtained $s$-wave gap
equation by considering the BCS limit.

###  IV.1 BCS limit

Equation (16) in the limiting case of vanishing bose-fermionic coupling
($\alpha\to 0$, $\beta\to 0$) yields the standard BCS gap equation for
attractive interactions in the fermionic sector ($g_{F}<0$)

$\Delta(\mathbf{k})=\sum\limits_{\mathbf{k^{\prime}}\neq
0}-\frac{g_{F}}{\Omega}\frac{\tanh\frac{E_{\mathbf{k}-\mathbf{k^{\prime}}}}{2k_{\textrm{B}}T}}{2E_{\mathbf{k}-\mathbf{k^{\prime}}}}\,\,\Delta(\mathbf{k}-\mathbf{k^{\prime}}).$
(33)

In the zero temperature limit $T\rightarrow 0$, we obtain the well-known BCS
gap equation at zero temperature

$\Delta=-\frac{\Delta
g_{F}}{2(2\pi)^{3}}\int\frac{d\textbf{k}}{\sqrt{\xi^{2}+\Delta^{2}}}=-\frac{\Delta
g_{F}M_{F}p_{F}}{2\pi^{2}\hbar^{3}}\ln\frac{\hbar\omega_{D}}{\Delta},$ (34)

where $\hbar\omega_{D}$ is the maximum exchange energy. For $g_{F}>0$, the
only possible solution is trivial solution $\Delta=0$. However, if $g_{F}<0$,
we obtain the well known BCS-type result

$\Delta=\hbar\omega_{D}e^{-\frac{1}{\zeta}},$

where $\zeta=|g_{F}|M_{F}p_{0}/(2\pi^{2}\hbar^{3}).$

For a repulsive fermion-fermion interaction ($g_{F}\geq 0$), the electron-
electron pairing is possible only due to the effective interaction induced by
Bogoluibov excitations in the bosonic sector. For $g_{F}>0$, the $s$-wave
pairing could be realized only if $\tilde{g}_{n}$ is sufficiently large. The
$s$-wave gap equation for the $s$-wave channel implies the following
inequalities for the interaction kernel:

$-\frac{g_{F}}{\Omega}+\tilde{g}_{n}(\textbf{k}^{\prime},\omega_{k}-\omega_{n})-\tilde{g}_{s}(\textbf{k}^{\prime},\omega_{k}-\omega_{n})\leq-\frac{g_{F}}{\Omega}+\tilde{g}_{n}(\textbf{k}^{\prime},\omega_{k}-\omega_{n})\leq-\frac{g_{F}}{\Omega}+\tilde{g}_{n}(0,\textbf{0})=-\frac{g_{F}}{\Omega}+\frac{\alpha^{2}}{4\Omega
nc_{0}}.$ (35)

It is obvious that the $s$-wave channel is suppressed in systems with
$g_{F}>0$.

It is remarkable that the density-density bose-fermi interaction (the term
proportional to $\alpha^{2}$ in Eq. (16)) produces effective attractive
interactions in the $s$-wave channel. At the same time the spin-spin bose-
fermi interactions (the term proportional to $\beta^{2}$ in Eq. (16)) lead to
an additional repulsive interaction for the $s$-wave channel. However, in
sharp contrast to the $s$-wave superfluidity, for the $p$-wave superfluidity,
both spin-spin and density-density bose-fermi interactions lead to an
effective attraction that results in the realization of the $p$-wave
superfluid.

Figure 1: The maximum value of the energy gap $\Delta_{\textrm{max}}$ of the
$s$-wave state (in units of the Fermi energy $E_{F}$) as a function of
temperature $T$ (in units of the Fermi temperature
$T_{F}=E_{F}/k_{\textrm{B}}$) found numerically in the BCS limit
$\alpha=\beta=0$. The inset represents the energy gap $\Delta$ (in units of
the Fermi energy $E_{F}$) as a function of momentum $k$ (in units of the Fermi
momentum $k_{F}=(3\pi^{2}n_{F})^{1/3}$) for temperature values indicated by
colored circles on the $\Delta_{\textrm{max}}$ curve. Figure 2: The same as in
Fig. 1 for $\alpha=0.3c_{0}$, $\beta=0.3c_{0}$.

###  IV.2 Numerical solution

We solve the gap equation for the $s$-wave state (32) by using the following
approximation for momentum-dependent functions:

$f(\mathbf{k}-\mathbf{k^{\prime}})=f(k)\mathbf{\theta}(k-k^{\prime})+f(k^{\prime})\mathbf{\theta}(k^{\prime}-k),$
(36)

where $\mathbf{\theta}$ is the Heaviside step function and

$E_{k}=\sqrt{\xi^{2}_{k}+\Delta(k)^{2}}.$ (37)

To solve the gap equation (32) numerically, we use the following stabilized
relaxation procedure:

$\Delta^{(n+1)}(\mathbf{k})=\int\mathcal{K}^{(n)}\times\Delta^{(n)}(\mathbf{k}-\mathbf{k^{\prime}})d\mathbf{k},$
(38)
$\mathcal{K}^{(n)}=-\frac{g_{F}}{\Omega}\frac{\tanh\frac{E_{\mathbf{k}-\mathbf{k^{\prime}}}}{2k_{\textrm{B}}T}}{2E_{\mathbf{k}-\mathbf{k^{\prime}}}}+c_{n}\left(\frac{W_{0}}{2E_{\mathbf{k}-\mathbf{k^{\prime}}}}\coth\frac{E_{\mathbf{k}-\mathbf{k^{\prime}}}}{2k_{\textrm{B}}T}+\frac{Z_{0}}{2E_{\mathbf{k},0}}\tanh\frac{E_{\mathbf{k},0}}{2k_{\textrm{B}}T}\right)\\\
-c_{s}\left(\frac{W_{1}}{2E_{\mathbf{k}-\mathbf{k^{\prime}}}}\coth\frac{E_{\mathbf{k}-\mathbf{k^{\prime}}}}{2k_{\textrm{B}}T}+\frac{Z_{1}}{2E_{\mathbf{k},1}}\tanh\frac{E_{\mathbf{k},1}}{2k_{\textrm{B}}T}\right).$
(39)

The subsequent iteration is calculated by using the above equation multiplying
$\Delta^{(n+1)}(\mathbf{k})$ by the ratio
$(\mathcal{N}_{n}/\mathcal{N}_{n+1})^{1/2}$, where
$\mathcal{N}_{n}=\int\Delta^{(n)}(\mathbf{k}-\mathbf{k^{\prime}})d\mathbf{k}$.
Note that there are no solutions of the $s$-wave gap equation for repulsive
fermi-fermi interactions ($g_{F}>0$) for considered here parameters $\alpha$
and $\beta$. Furthermore, for attractive fermi-fermi interactions ($g_{F}<0$)
there is the threshold value of the interaction constant $g_{F}$ which allows
the existence of $s$-wave superfluid state.

Our numerical results are presented in Figs.1 and 2. The inset in Fig. 1
illustrates typical gap profiles for the case $\alpha=\beta=0$ as a function
of momentum at different temperature. As one can expect, the maximum is
located near the Fermi momentum $k=k_{F}$. The gap amplitude decreases with
temperature and vanishes at the critical point. Figure 2 shows typical
solutions of the gap equation with nonzero density-density and spin-spin
interactions for $\alpha=\beta=0.3c_{0}$. Notice that the numerical method
used in the present work does not allow us to determine the gap values near
the $k=k_{F}$ as is illustrated in the inset of Fig. 2. Nevertheless, it is
still possible to find the amplitude value as a function of temperature and
determine the critical temperature when the pairing gap vanishes. Note that
spin-spin and density-density interactions decrease of the energy gap of the
$s$-wave state.

##  V Conclusions

The $s$\- and $p$-wave superfluidity is investigated in a mixture of fermionic
40K and $F=1$ spinor bosonic 23Na gases. As is known, in common superfluids,
the $s$-wave phase dominates because it has a higher critical temperature. In
the present work, we focused on the problem of the existence of the $s$-wave
superfluid phase of fermions in Bose-Fermi mixtures and determined the
physical conditions when the $p$-wave state is more favorable. The gap
equations for the $s$-wave and the polar phase of $p$-wave superfluid fermions
are derived. The gap equation for the $s$-wave state is solved numerically for
different values of temperature. As one can expect, the maximum value of the
gap is located near the Fermi momentum $k=k_{F}$ and the energy gap of the
$s$-wave state is suppressed by spin-spin interactions of bosonic and
fermionic atoms. The physical conditions necessary for the realization of the
$s$-wave and the polar phase in a well-controlled environment of atomic
physics are determined. Our analysis shows that a nontrivial solution to the
$s$-wave gap equation is absent for repulsive fermi-fermi interactions
($g_{F}>0$) in the case of physically relevant parameters. As to attractive
fermi-fermi interactions ($g_{F}<0$), the $s$-wave superfluid state exists
only for sufficiently strong fermion-fermion interaction strength $g_{F}$.
Thus, the $s$-wave state is suppressed compared to the $p$-wave state in Bose-
Fermi mixtures. We hope that the obtained results can be useful for the
practical realization of the $p$-wave state in ultracold atomic gases.

## ACKNOWLEDGMENTS

The authors are grateful to O.V. Bugaiko for the participation in the initial
stage and substantial contributions to this work. The work of E.V.G. was
partially supported by the Program of Fundamental Research of the Physics and
Astronomy Division of the National Academy of Sciences of Ukraine. Y.O.N. and
A.I.Y. acknowledge partial support from the National Research Foundation of
Ukraine through Grant No. 2020.02/0032.

## References

  * Nayak et al. (2008) C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. Das Sarma, Rev. Mod. Phys. 80, 1083 (2008).
  * Inotani et al. (2016) D. Inotani, R. Hanai, and Y. Ohashi, Phys. Rev. A 94, 043632 (2016).
  * Wu and Bruun (2016) Z. Wu and G. M. Bruun, Phys. Rev. Lett. 117, 245302 (2016).
  * Cooper and Shlyapnikov (2009) N. Cooper and G. Shlyapnikov, Phys. Rev. Lett. 103, 155302 (2009).
  * Zhang et al. (2008) C. Zhang, S. Tewari, R. M. Lutchyn, and S. D. Sarma, Phys. Rev. Lett. 101, 160401 (2008).
  * Juliá-Díaz et al. (2013) B. Juliá-Díaz, T. Graß, O. Dutta, D. Chang, and M. Lewenstein, Nature communications 4, 2046 (2013).
  * Han et al. (2009) Y.-J. Han, Y.-H. Chan, W. Yi, A. Daley, S. Diehl, P. Zoller, and L.-M. Duan, Phys. Rev. Lett. 103, 070404 (2009).
  * Bardyn et al. (2012) C.-E. Bardyn, M. Baranov, E. Rico, A. İmamoğlu, P. Zoller, and S. Diehl, Phys. Rev. Lett. 109, 130402 (2012).
  * Thompson et al. (2020) K. Thompson, J. Brand, and U. Zülicke, Phys. Rev. A 101, 013613 (2020).
  * Sato et al. (2009) M. Sato, Y. Takahashi, and S. Fujimoto, Phys. Rev. Lett. 103, 020401 (2009).
  * Leggett (1975) A. Leggett, Rev. Mod. Phys. 47, 331 (1975).
  * Sigrist and Ueda (1991) M. Sigrist and K. Ueda, Rev. Mod. Phys. 63, 239 (1991).
  * Stewart (1984) G. Stewart, Rev. Mod. Phys. 56, 755 (1984).
  * MacKenzie and Maeno (2003) A. P. MacKenzie and Y. Maeno, Rev. Mod. Phys. 75, 657 (2003).
  * Wu et al. (2020) Y. Wu, Z. Yan, Z. Lin, J. Lou, and Y. Chen, Scientific Reports 10, 10822 (2020), eprint 2001.00420.
  * Bazak and Petrov (2018) B. Bazak and D. S. Petrov, Phys. Rev. Lett. 121, 263001 (2018).
  * Wang et al. (2018) J.-B. Wang, W. Yi, and J.-S. Pan, Phys. Rev. A 98, 053630 (2018).
  * Efremov and Viverit (2002) D. V. Efremov and L. Viverit, Phys. Rev. B 65, 134519 (2002).
  * Regal et al. (2003) C. A. Regal, C. Ticknor, J. L. Bohn, and D. S. Jin, Phys. Rev. Lett. 90, 053201 (2003).
  * Levinsen et al. (2007) J. Levinsen, N. R. Cooper, and V. Gurarie, Phys. Rev. Lett. 99, 210402 (2007).
  * Levinsen et al. (2008) J. Levinsen, N. R. Cooper, and V. Gurarie, Phys. Rev. A 78, 063616 (2008).
  * Jona-Lasinio et al. (2008) M. Jona-Lasinio, L. Pricoupenko, and Y. Castin, Phys. Rev. A 77, 043611 (2008).
  * Matsyshyn et al. (2018) O. I. Matsyshyn, A. I. Yakimenko, E. V. Gorbar, S. I. Vilchinskii, and V. V. Cheianov, Phys. Rev. A 98, 043620 (2018).
  * Ho (1998) T.-L. Ho, Phys. Rev. Lett. 81, 742 (1998).
  * Kawaguchi and Ueda (2012) Y. Kawaguchi and M. Ueda, Phys. Rep. 520, 253 (2012).