# Vortical versus skyrmionic states in the topological phase of a twisted
bilayer with $d$-wave superconducting pairing

Leonardo R. Cadorim  Departamento de Física, Faculdade de Ciências,
Universidade Estadual Paulista (UNESP), Caixa Postal 473, 17033-360, Bauru-SP,
Brazil  Department of Physics, University of Antwerp, Groenenborgerlaan 171,
B-2020 Antwerp, Belgium  Edson Sardella  Departamento de Física, Faculdade de
Ciências, Universidade Estadual Paulista (UNESP), Caixa Postal 473, 17033-360,
Bauru-SP, Brazil  Milorad V. Milošević<EMAIL_ADDRESS>Department of Physics, University of Antwerp, Groenenborgerlaan 171, B-2020
Antwerp, Belgium  NANOlab Center of Excellence, University of Antwerp, Belgium

###### Abstract

It was recently shown that a chiral topological phase emerges from the
coupling of two twisted monolayers of superconducting Bi2Sr2CaCu2O8$+\delta$
for certain twist angles. In this work, we reveal the behavior of such twisted
superconducting bilayers with $d_{x^{2}-y^{2}}$ pairing symmetry in presence
of applied magnetic field. Specifically, we show that the emergent vortex
matter can serve as smoking gun for detection of topological superconductivity
in such bilayers. Moreover, we report two distinct skyrmionic states that
characterize the chiral topological phase, and provide full account of their
experimental signatures and their evolution with the twist angle.

††preprint: APS/123-QED

##  I Introduction

Chiral superconductivity [1] has been a topic of tremendous interest in the
recent literature due to its rich phenomenology [2, 3], including appearance
of nontrivial surface currents [4] and half-quantum vortices [5, 6, 7, 8, 9],
to name a few examples. Being mostly characterized by several Fermi surfaces,
chiral superconductors often present multiple superconducting gaps, and are
thereby prone to a plethora of interesting physics typical of multicomponent
superconductivity [10, 11, 12, 13, 14]. Arguably, chiral superconductors
gained a special relevance due to the increasing interest in topological
superconductivity [15] and its promise towards use in modern quantum
technologies [16]. With its highly non-trivial topology, the chiral state of
superconductors is known to present the uniquely associated phenomena, such as
the gapless edge states [17] and Majorana bound states localized in the vortex
cores [18], which obey the non-Abelian statistics [19] fundamental to future
applications in quantum computing.

Recently, Can et al. [20] showed that a twisted bilayer composed of two
monolayers of the high-temperature superconductor Bi2Sr2CaCu2O8$+\delta$ [21]
(Bi-2212) can display a chiral topological phase which breaks time-reversal
symmetry for twist angles near $45^{\circ}$. As they argued, at a twist angle
equal to $45^{\circ}$, the $d_{x^{2}-y^{2}}$ order parameter of each layer,
characteristic of Bi-2212, induces a significant $d_{xy}$ component in the
order parameter of the other layer. This results in a superconducting state
with $d+id^{\prime}$ pairing symmetry.

The above arguments were developed in Ref. 20 considering a homogeneous
superconducting state. In the present work we go beyond this premise and
investigate how such system responds to applied magnetic field, i.e. how the
vortex matter of such bilayers evolves with the twist angle between the
monolayers. As we will show, the emergent typical vortex configurations can be
used as a smoking gun for the detection of the chiral topological phase.
Namely, due to the broken time-reversal symmetry, skyrmionic vortex states
[22, 9, 23, 24] arise in the topological phase, and are clearly identifiable
by their magnetic signature. Moreover, we show that the vortex matter changes
even within the topological phase itself. Namely, as one varies the twist
angle in the range where topological phase is stable, two different skyrmionic
states are found. In one of them, states with unit topological charge are
favored (presenting as a lattice of vortex pairs), while in the other one
states with large topological charge become energetically favorable, causing
formation of extended vortex chains with a distinct appearance and magnetic
signature.

The outline of this work is as follows. In Sec. II we present our theoretical
formalism and show how we deal with the twisted bilayer system at hand within
the framework of the Ginzburg-Landau theory. In Sec. III we present and
discuss our main results. We start by showing the existence of a topological
phase for certain values of the twist angle in a homogeneous system, to
subsequently reveal and characterize the vortical and skyrmionic matter, as
well as transitions between them, inside the topological phase. Our concluding
remarks are given in Sec. IV.

##  II The theoretical model

The free energy density of our system can be described as a sum of three parts
$\mathcal{F}=F_{1}+F_{2}+F_{12}$, with

$\displaystyle F_{1}$ $\displaystyle=$
$\displaystyle-2\alpha_{s}|\Delta_{(s1)}|^{2}-|\Delta_{(d1)}|^{2}+\frac{4}{3}|\Delta_{(s1)}|^{4}+\frac{1}{2}|\Delta_{(d1)}|^{4}$
(1)
$\displaystyle+\frac{8}{3}|\Delta_{(s1)}|^{2}|\Delta_{(d1)}|^{2}+\frac{2}{3}(\Delta^{*2}_{(s1)}\Delta^{2}_{(d1)}+H.c.)$
$\displaystyle+2|\bm{\Pi}\Delta^{*}_{(s1)}|^{2}+|\bm{\Pi}\Delta^{*}_{(d1)}|^{2}+(\Pi_{x}\Delta^{*}_{(s1)}\Pi^{*}_{x}\Delta_{(d1)}$
$\displaystyle-\Pi_{y}\Delta^{*}_{(s1)}\Pi^{*}_{y}\Delta_{(d1)}+H.c.)$

being the free energy of the unrotated (non-twisted) layer. Here,
$\Delta_{(s1)}$ and $\Delta_{(d1)}$ are the order parameters corresponding to
the $s$ and $d$-wave pairings, respectively [25]. Hereafter, the subscript
$(1)$ indicates the order parameters of the unrotated layer, while the
subscript $(2)$ denotes the rotated (twisted) one. In this work, we add the
$s$ component of the order parameter in order to induce the correct fourfold
symmetry in the $d$ component. This is done through the mixed gradient terms
in the free-energy. The parameter $\alpha_{s}$ determines the relative
strength between the $s$ and $d$ order parameters. Once we are mainly
interested in the condensate with $d$-wave pairing, we use $\alpha_{s}=0.7$
that leads to a weak modulus for the $s$-wave order parameter (under 20% of
the $d$-wave order parameter). We also define the momentum operator
$\bm{\Pi}=i\bm{\nabla}-\bm{A}$ for compacter presentation of the formulae.
Here $\bm{A}$ stands for the magnetic vector potential due to the applied
magnetic field and the magnetic response of the superconducting layers.

The second contribution to the free energy stems from the rotated layer, and
reads

$\displaystyle F_{2}$ $\displaystyle=$
$\displaystyle-2\alpha_{s}|\Delta_{(s2)}|^{2}-|\Delta_{(d2)}|^{2}+\frac{4}{3}|\Delta_{(s2)}|^{4}+\frac{1}{2}|\Delta_{(d2)}|^{4}$
(2)
$\displaystyle+\frac{8}{3}|\Delta_{(s2)}|^{2}|\Delta_{(d2)}|^{2}+\frac{2}{3}(\Delta^{*2}_{(s2)}\Delta^{2}_{(d2)}+H.c.)$
$\displaystyle+2|\bm{\Pi}\Delta^{*}_{(s2)}|^{2}+|\bm{\Pi}\Delta^{*}_{(d2)}|^{2}+\cos(2\theta)\Pi_{x}\Delta^{*}_{(s2)}\Pi^{*}_{x}\Delta_{(d2)}$
$\displaystyle-\cos(2\theta)\Pi_{y}\Delta^{*}_{(s2)}\Pi^{*}_{y}\Delta_{(d2)}-\sin(2\theta)\Pi_{x}\Delta^{*}_{(s2)}\Pi^{*}_{y}\Delta_{(d2)}$
$\displaystyle-\sin(2\theta)\Pi_{y}\Delta^{*}_{(s2)}\Pi^{*}_{x}\Delta_{(d2)}+H.c..$

Here, $\theta$ is the twist angle and the expression for the mixed gradient
terms presented in Eq. (2) is obtained after transformation on the momentum
operator from the rotated coordinates to the unrotated ones.

The final contribution to the free energy captures the interaction between the
two layers [20]

$\displaystyle F_{12}=$ $\displaystyle A$
$\displaystyle|\Delta_{(d1)}|^{2}|\Delta_{(d2)}|^{2}-B\cos(2\theta)(\Delta_{(d1)}\Delta^{*}_{(d2)}+H.c.)$
$\displaystyle+$ $\displaystyle C$
$\displaystyle(\Delta^{2}_{(d1)}\Delta^{*2}_{(d2)}+H.c.),$ (3)

with $A$, $B$ and $C$ taken as phenomenological (free) parameters. The term
proportional to $B$ in Eq. (3) depends on $\cos(2\theta)$ due to symmetry
reasons, as discussed in Ref. [20], and represents the tunneling of Cooper
pairs between the two layers. Following same reasoning, one can interpret the
term proportional to $C$ as the coherent tunnelling of two Cooper pairs
between the layers.

In the above equations, all lengths are expressed in units of the coherence
length $\xi=(\nu_{F}/2)\sqrt{W/\ln(T_{d}/T)}$, with
$W=7\xi(3)/(8\pi^{2}T^{2})$, the order parameters are in units of
$\Delta_{0}=\sqrt{(4/3W)\ln(T_{d}/T)}$, the magnetic field is in units of
$H_{c2}=\Phi_{0}/(2\pi\xi^{2})$, where $\Phi_{0}=hc/2e$ is the magnetic flux
quantum, and the free energy density is in units of
$F_{0}=(4/3W)\ln(T_{d}/T)$. For details on the derivation of the free energy
for a single layer we refer to Refs. [25,26].

Minimizing the total energy $\mathcal{F}=F_{1}+F_{2}+F_{12}$ with respect to
the order parameters we arrive to the appropriate Ginzburg-Landau equations:

$\displaystyle-\alpha_{s}\Delta_{(s1)}+\frac{4}{3}|\Delta_{(s1)}|^{2}\Delta_{(s1)}+\frac{4}{3}|\Delta_{(d1)}|^{2}\Delta_{(s1)}+\frac{2}{3}\Delta^{2}_{(d1)}\Delta^{*}_{(s1)}$
$\displaystyle+\Pi^{*2}\Delta_{(s1)}+\frac{1}{2}(\Pi^{*2}_{x}-\Pi^{*2}_{y})\Delta_{(d1)}=0,$
(4)
$\displaystyle-\Delta_{(d1)}+|\Delta_{(d1)}|^{2}\Delta_{(d1)}+\frac{8}{3}|\Delta_{(s1)}|^{2}\Delta_{(d1)}+\frac{4}{3}\Delta^{2}_{(s1)}\Delta^{*}_{(d1)}$
$\displaystyle+A|\Delta_{(d2)}|^{2}\Delta_{(d1)}-B\cos(2\theta)\Delta_{(d2)}+2C\Delta^{2}_{(d2)}\Delta^{*}_{(d1)}$
$\displaystyle+\Pi^{*2}\Delta_{(d1)}+(\Pi^{*2}_{x}-\Pi^{*2}_{y})\Delta_{(s1)}=0,$
(5)
$\displaystyle-\alpha_{s}\Delta_{(s2)}+\frac{4}{3}|\Delta_{(s2)}|^{2}\Delta_{(s2)}+\frac{4}{3}|\Delta_{(d2)}|^{2}\Delta_{(s2)}$
$\displaystyle+\frac{2}{3}\Delta^{2}_{(d2)}\Delta^{*}_{(s2)}+\frac{1}{2}(\cos(2\theta)\Pi^{*2}_{x}-\cos(2\theta)\Pi^{*2}_{y})\Delta_{(d2)}$
$\displaystyle-\frac{1}{2}\sin(2\theta)\Pi^{*}_{x}\Pi^{*}_{y}\Delta_{(d2)}-\frac{1}{2}\sin(2\theta)\Pi^{*}_{y}\Pi^{*}_{x}\Delta_{(d2)}$
$\displaystyle+\Pi^{*2}\Delta_{(s2)}=0,$ (6)

and

$\displaystyle-\Delta_{(d2)}+|\Delta_{(d2)}|^{2}\Delta_{(d2)}+\frac{8}{3}|\Delta_{(s2)}|^{2}\Delta_{(d2)}+\frac{4}{3}\Delta^{2}_{(s2)}\Delta^{*}_{(d2)}$
$\displaystyle+A|\Delta_{(d1)}|^{2}\Delta_{(d2)}-B\cos(2\theta)\Delta_{(d1)}+2C\Delta^{2}_{(d1)}\Delta^{*}_{(d2)}$
$\displaystyle+\Pi^{*2}\Delta_{(d2)}+(\cos(2\theta)\Pi^{*2}_{x}-\cos(2\theta)\Pi^{*2}_{y})\Delta_{(s2)}$
$\displaystyle-\sin(2\theta)\Pi^{*}_{x}\Pi^{*}_{y}\Delta_{(s2)}-\sin(2\theta)\Pi^{*}_{y}\Pi^{*}_{x}\Delta_{(s2)}=0.$
(7)

Eqs. (4)-(7) are then solved for different twist angles $\theta$, assuming
periodic boundary conditions [27]. As the Ginzburg-Landau parameter $\kappa$
for Bi-2212 is typically much greater than 1, we disregard the contribution of
the supercurrents to the total magnetic field and use a vector potential
solely due to the applied magnetic field to solve the above set of equations.

For a given $\theta$, we initialize the calculations from dozens of different
initial conditions for the order parameters, while also varying the aspect
ratio of the unit cell of the simulation in order to identify the lowest-
energy solutions for the vortex states. In what follows, we display the
ground-state found for an external applied magnetic field that corresponds to
the flux of $24\Phi_{0}$ threading the shown unit cell, without loss of
generality. Namely, for other values of the applied magnetic field we obtained
qualitatively equivalent results.

##  III Results and Discussion

As shown in Ref. [20] in the case of homogeneous superconductivity, for a
certain range of $\theta$ the competition between the terms proportional to
$B$ and $C$ in $F_{12}$ yields a non-trivial phase difference between the
$d$-wave components of the order parameters of the two layers. Let us start by
discussing the homogeneous solutions of our free-energy model and show that it
analogously allows for the existence of a topological phase.

Figure 1: (Color online) Phase difference between the $d$-wave components of
the order parameters of two layers as a function of the twist angle $\theta$.
Blue and red curves represent the phase difference for $C=B/8$ and $C=B/5$,
respectively, while $A=B=0.1$. The nontrivial values of the phase difference
($\neq 0$ or $\pi$) indicate existence of a topological phase for a particular
twist angle.

To do this, we minimize the free energy density
$\mathcal{F}=F_{1}+F_{2}+F_{12}$ with respect to the modulus and phase of the
$s$ and $d$ components of the order parameter in both layers. Fig. 1 shows the
phase difference between $\Delta_{(d1)}$ and $\Delta_{(d2)}$ which minimizes
the free energy as a function of $\theta$. As can be seen from the figure, for
small twist angles up to a critical angle $\theta_{i}$, the free energy is
minimal when the order parameters have the same phase. For twist angles larger
than a critical value $\theta_{f}$, the phase difference that yields minimal
energy equals $\pi$.

On the other hand, for angles between $\theta_{i}$ and $\theta_{f}$, one
obtains a non-trivial phase difference between the condensates of the two
layers, which means a superconducting state that breaks time-reversal
symmetry. In particular, for $\theta=\pi/4$, the phase difference is equal to
$\pi/2$, i.e. a $d+id^{\prime}$ superconducting state is found [20]. The
values of $\theta_{i}$ and $\theta_{f}$ depend on the particular values chosen
for the parameters $A$, $B$ and $C$, as can be seen from the two examples
shown in Fig. 1. Nevertheless, the features of the superconducting state that
we discuss below are always present in the range
$\theta_{i}<\theta<\theta_{f}$, for any choice of the aforementioned
parameters. Therefore, without loss of generality of our results, in what
follows we will use the parameters correspondent to the blue curve in Fig. 1.
In that case, $\theta_{i}\approx 34^{\circ}$ and $\theta_{f}\approx
56^{\circ}$.

Figure 2: (Color online) Vortex configurations in the $d$-wave component of
the order parameter of the unrotated layer (first column) and the rotated
layer (second column), and the magnetic response of the system (third column),
for three selected twist angles. Figure 3: (Color online) Vortex
configurations for twist angle $\theta=36^{\circ}\gtrsim\theta_{i}$ at three
values of applied magnetic field. Each row corresponds to a different system
size (lateral sizes are shown), for fixed magnetic flux of $24\Phi_{0}$. From
left to right, the columns respectively show the magnitude of the $d$-wave
component of the order parameter for unrotated and rotated layers, the cosine
and sine of the phase difference between the order parameters in two layers,
and the magnetic field distribution across the system.

Let us now go beyond these results and show how the twist angle affects the
vortex matter of such bilayers. In Fig. 2 we show the spatial distribution of
the $d$-wave component of the order parameter for both layers, together with
the magnetic field distribution in the system, for $\theta=0$, $36^{\circ}$
and $45^{\circ}$. When plotting the magnetic field profile, we consider only
the contribution of the supercurrents, after subtracting the (strongly
dominating) homogeneous external field from the total one. For $\theta=0$, as
discussed above, the phase difference between the order parameters is locked
at zero. It is therefore energetically favorable for the vortices in two
layers to organize in a composite configuration, where the vortex cores are
vertically aligned between the unrotated and the rotated layer. In this case,
the magnetic field profile of the vortices exhibits well defined peaks at
vortex locations and one can clearly distinguish the fourfold symmetry
characteristic of $d_{x^{2}-y^{2}}$ superconductors. The vortex configurations
remain qualitatively unchanged for the non-zero twist angles outside the
topological phase, i.e. for $\theta<\theta_{i}$ or $\theta>\theta_{f}$.

Next we increase the twist angle to $36^{\circ}$, larger than $\theta_{i}$,
i.e. upon entry to the topological phase in the ground state of the system.
The composite vortex configuration is no longer the most energetically
favorable state, as the vortex cores in two layers no longer coincide. As a
consequence, the field of the vortex is now distributed over the two displaced
cores in two layers (cf. Fig. 1), reflecting a magnetic field profile of a
dimer rather than one clear peak. Notably, some vortices within the
configuration remain seemingly composite. The situation radically changes as
we increase the twist angle further, to $\theta=45^{\circ}$ (bottom row of
Fig. 1). Here one sees that not only vortex cores displace between the layers,
they also organize into extended closed vortex chains. As will be discussed in
the following paragraphs, such a vortex chain is formed along a domain wall
separating sample regions with different phase differences between the layers.
Moreover, such chains will exhibit skyrmionic topology, with an integer
topological charge equal to the total vorticity of the chain. Last but not
least, the overall shape and the magnetic signature of these chains are
uniquely distinct which facilitates their experimental observation.

These results demonstrate that the twist angle and the onset of a topological
phase strongly influence the vortex matter of the system, with detectable
consequences in the magnetic profile at the onset of the topological phase and
within the topological phase itself. This feature can therefore be used as a
smoking gun for the detection of topological superconductivity in such and
similar bilayers. In what follows, we further detail the vortex configurations
for $\theta=36^{\circ}$ and $\theta=45^{\circ}$, which are the representative
examples of two different types of behavior we encountered in the vortex
matter inside the topological phase.

Figure 4: (Color online) Vortex configurations for $\theta=45^{\circ}$, deep
inside the topological phase, for three values of applied magnetic field. Each
row corresponds to a different system size (lateral sizes are shown), for
fixed magnetic flux of $24\Phi_{0}$. From left to right, the columns
respectively show the magnitude of the $d$-wave component of the order
parameter for unrotated and rotated layers, the cosine and sine of the phase
difference between the order parameters in two layers, and the magnetic field
distribution across the system. Figure 5: (Color online) Energy of the
homogeneous system (without vortices) as a function of the phase difference
between the $d$-wave components of the order parameter in two twisted layers.
The blue, yellow and red curves show the energy for the twist angle
$\theta=36^{\circ}$, $\theta=42^{\circ}$ and $\theta=45^{\circ}$,
respectively. Figure 6: (Color online) Evolution of the vortex configuration
when adiabatically decreasing the twist angle from $\theta=45^{\circ}$ to
$36^{\circ}$. From left to right, the columns respectively show the magnitude
of the $d$-wave component of the order parameter for unrotated and rotated
layers, the cosine and sine of the phase difference between the order
parameters in two layers, and the magnetic field distribution across the
system. Figure 7: (Color online) Evolution of the vortex configuration when
adiabatically increasing the twist angle from $\theta=36^{\circ}$ to
$45^{\circ}$. From left to right, the columns respectively show the magnitude
of the $d$-wave component of the order parameter for unrotated and rotated
layers, the cosine and sine of the phase difference between the order
parameters in two layers, and the magnetic field distribution across the
system. Figure 8: (Color online) Energy as a function of the twist angle
$\theta$. The blue and red lines show the energy curves for the twist angle
being decreased from $\theta=45^{\circ}$ to $36^{\circ}$ and increased from
$36^{\circ}$ to $45^{\circ}$, respectively.

###  III.1 Vortex matter at the onset of the topological phase

We start the description of the vortex matter for an angle close to
$\theta_{i}$ (equivalent results are obtained for angles close to
$\theta_{f}$). In Fig. 3 we show vortex configurations found for
$\theta=36^{\circ}$, each row displaying minimum energy solutions for a
different size of the unit cell. As discussed previously, the vortex cores in
two layers are displaced from one another inside the topological phase, but
each vortex of a given layer remains connected to its counterpart in the other
layer. This is visible in the phase difference profile, suggesting existence
of a phase soliton between the two vortex cores. Hereafter, we refer to this
pair of connected vortices as the interlayer vortex pair. As can be seen in
Fig. 3, inside the interlayer vortex pair we find phase difference
$\theta_{d2}-\theta_{d1}=\pi$ between the condensates in two layers, and
$\theta_{d2}-\theta_{d1}=0$ outside of the pair.

As one object, the interlayer vortex pair displays skyrmionic properties,
which can be described by first defining the pseudo-spin [28]

$\bm{n}=\frac{\Delta_{d}^{\dagger}\bm{\sigma}\Delta_{d}}{\Delta_{d}^{\dagger}\Delta_{d}},$
(8)

with $\Delta_{d}=(\Delta_{(d1)},\Delta_{(d1)})$ and
$\bm{\sigma}=(\sigma_{1},\sigma_{2},\sigma_{3})$, where $\sigma_{i}$ is the
Pauli matrix. With such pseudo-spin profile, one then calculates the
topological charge of the system, defined as

$Q=\frac{1}{4\pi}\int\bm{n}\cdot\left(\frac{\partial\bm{n}}{\partial
x}\times\frac{\partial\bm{n}}{\partial y}\right)dxdy.$ (9)

If we calculate the total topological charge for the three different
configurations shown in Fig. 3, we obtain $Q=24$ for each of them. As $24$ is
also the number of flux quanta threading the shown unit cells, we conclude
that each interlayer vortex pair is actually a skyrmionic object with a
topological charge equal to $1$ (for a truly composite vortex, $Q=0$). We
confirmed this further by calculating the topological charge not over the
entire unit cell, but only around isolated interlayer pairs.

Within the dashed rectangles shown in Fig. 3 we highlight vortex
configurations characteristic of the topological phase for twist angles close
to $\theta_{i}$. As can be seen from the cosine and sine of the phase
difference between the condensate of each layer, two adjacent interlayer
vortex pairs typically organize themselves into a larger correlated object.
Inside the dashed rectangles in panels $(a1)-(a2)$ of Fig. 3, we can see that
in one of the interlayer vortex pair (the one near the top of the rectangle)
the vortex of the unrotated layer is on the right and the vortex of the
rotated layer on the left. In the interlayer pair near the bottom of the
rectangle, the vortex positions are interchanged. In other words, adjacent
pairs of vortices in two layers are twisted with respect to each other. Same
behavior can be easily verified in panels $(b1)-(b2)$ and $(c1)-(c2)$. After
such organization of interlayer vortex pairs, their interlayer phase solitons
become visibly connected, as seen in the dashed rectangle of the cosine of the
phase difference in Fig. 3. Due to proximity and strong overlap between
different phase domains, a supercurrent emerges surrounding the two interlayer
vortex pairs, yields rather weak magnetic field. As a consequence, while the
field profile of each pair is strong, and spatial correlation between them is
rather obvious, the consequence of the phase connection between them is
difficult to spot in the last column of Fig. 3.

###  III.2 Vortex matter deep inside the topological phase

Finally we reveal the evolution of the vortex matter when the twisted bilayer
is deeply inside the topological phase, i.e. for twist angle $\theta\approx
45^{\circ}$ in the present case. As previously shown, for this $\theta$ the
phase difference between the $d$-wave components of the order parameters in
two layers is $\pi/2$ and we have a $d+id^{\prime}$ superconducting state.
Despite this particularity, the encountered characteristics of the vortex
matter in this case can be related to the ones exhibited for other twist
angles within the topological phase in the vicinity of $\theta=45^{\circ}$.

Fig. 4 shows the order parameters of both layers, the cosine and the sine of
the phase difference between the layers and the magnetic field profile for
$\theta=45^{\circ}$. In this case, formation of interlayer vortex pairs with
topological charge $Q=1$ is still favorable, as highlighted by dashed
rectangles in panels $(c1)-(c5)$ of Fig. 4. However, as highlighted by dashed
rectangles in panels $(a1)-(a5)$, the organization of interlayer vortex pairs
into larger objects is preferable. As a consequence, a new vortex
configuration emerges - the skyrmionic chain. In this uniquely distinct state,
instead of the small interlayer phase domains within individual interlayer
vortex pairs, much larger domains are formed. Interlayer vortex pairs are
interconnected along the domain wall, with a vortex core from one layer being
in between two vortex cores of the other one, forming a closed chain of
interlayer vortex pairs. Dashed rectangles in panels $(a1)-(a5)$ of Fig. 4
exemplify one such structure, containing 11 interlayer vortex pairs in a
single closed chain. In panel $(a5)$, one can see that such vortex chain is a
very laterally extended object (nearly $40\xi\times 40\xi$ in this case), and
leaves a very clear and rather unusual signature in the magnetic field profile
of the system. Here, the peaks of the magnetic field along the vortex chain
are smaller in comparison with the isolated interlayer vortex pairs due to the
fact that the distance between the vortex cores in two layers is significantly
larger in the former case.

Once again, and as can be seen from the third and fourth columns in Fig. 4,
the closed vortex chain separates two regions with different phase differences
between the layers. Inside the vortex chain, the phase difference between the
condensates is $\theta_{d2}-\theta_{d1}=-\pi/2$, while outside the chain
$\theta_{d2}-\theta_{d1}=\pi/2$. The opposite is also possible: panels
$(b1)-(b5)$ show such an example, where $\theta_{d2}-\theta_{d1}=\pi/2$ inside
the chain, while $\theta_{d2}-\theta_{d1}=-\pi/2$ outside. We note the
difference from the case of the individual interlayer vortex pairs, harboring
phase difference $\pi$ within them, with zero phase difference away from the
pair.

If we now calculate the total topological charge around the vortex chains seen
in panels $(a1)-(a5)$ and $(b1)-(b5)$, we obtain $Q=11$ and $Q=3$,
respectively, reflecting the number of interlayer vortex pairs interconnected
in the chain. These large values for the topological charge of such a novel
object contrast the exclusively $Q=1$ of the individual interlayer vortex
pairs found for $\theta=36^{\circ}$. This very different behavior for
different twist angles emerges from the fact that, as discussed before, the
domain wall separates regions with a phase difference equal to $0$ and $\pi$
for $\theta=36^{\circ}$ and $-\pi/2$ and $\pi/2$ for $\theta=45^{\circ}$.

To understand how the value of the phase difference inside the domains affects
the topological charge, Fig. 5 shows the energy of the homogeneous system as a
function of $\theta_{d2}-\theta_{d1}$ for different values of the twist angle
$\theta$. At the onset of the topological phase (blue curve in Fig. 5), the
energy of the system is largest when the phase difference is $\pi$. Therefore,
larger splitting within each formed interlayer vortex pair costs energy, and
their interconnection into larger objects is not energetically favorable.
Notice that, as discussed in Fig. 1, the phase difference equal to $0$ does
not yield the free energy minimum in the homogeneous state of the system for
$\theta$ at which the topological state is stable. In the presence of magnetic
field, the formation of interlayer vortex pairs re-stabilizes the zero phase
difference in a large part of the superconductor for $\theta$ close to
$\theta_{i}$ and $\theta_{f}$.

On the other hand, deep in the topological phase (for $\theta=45^{\circ}$), a
degenerate lowest energy homogeneous state is found for phase difference equal
to either $-\pi/2$ or $\pi/2$, explaining the tendency to formation of
coexisting domains with such phase differences. The resulting long domain
walls would cost energy, but not in the presence of magnetic field when they
are decorated by the skyrmionic vortex chains.

For $\theta$ values in the vicinity of $45^{\circ}$, represented in Fig. 5 by
$\theta=42^{\circ}$, the free-energy minima no longer occur at $-\pi/2$ and
$\pi/2$ but shift to lower phase differences and become shallower (cf. Fig.
5). Nevertheless, the system still presents the vortex chains dividing the
superconductor in regions with phase differences $-\pi/2$ and $\pi/2$, since
vortices require a total phase difference $\pi$ across the domain wall on
which they reside. Once both values of the phase difference possess the same
free energy, the long domain walls described above are also present, with the
same size as the ones for $\theta=45^{\circ}$.

###  III.3 Transitions between the topological vortex matter with the
interlayer twist

Complementary, it seems relevant to discuss in which manner the above-
described characteristic skyrmionic states in the topological phase evolve as
one continuously varies the twist angle. To capture this behavior, we follow
two distinct procedures. In the first, we start deep in the topological phase,
i.e. at a twist angle $\theta=45^{\circ}$ and a skyrmionic vortex chain as the
initial state of the simulation. We then “adiabatically” decrease the twist
angle down to $36^{\circ}$, in decrements of $0.1^{\circ}$, recording the
evolution of the stable solution (which is no longer necessarily the lowest-
energy state). In Fig. 6 we show the selected vortex configurations obtained
during this procedure. Starting from the skyrmionic vortex chain (panels
$(a1)-(a5)$), we see that when the twist angle is decreased to
$\theta=42,1^{\circ}$ (panels $(b1)-(b5)$), the contour of the chain can still
be seen in the magnetic field profile of the system. At the same time, the
sine of the phase difference shows that the vortex chain still splits the
superconducting landscape in regions with interlayer phase difference equal to
either $\pi/2$ and $-\pi/2$. However, vortices in each layer start to group in
pairs, as reflected in double peaks appearing in the magnetic field profile
along the chain. This indicates the onset of the transition from the
skyrmionic vortex chain to the skyrmionic state with separate interlayer
vortex pairs. Such a transition is more apparent for $\theta=38,9^{\circ}$
(panels $(c1)-(c5)$). Although the magnetic field contour of the vortex chain
can still be visualized in this case, the separation of the superconductor
into regions with different interlayer phase differences becomes less clear.
Finally, at $\theta=36^{\circ}$ (panels $(d1)-(d5)$) the transition between
the two skyrmionic states is completed and the system displays an arrangement
of dissociated individual interlayer vortex pairs, surrounded by a landscape
of near-zero interlayer phase difference.

Along the opposite route, we start with the lowest-energy state with
interlayer vortex pairs for $\theta=36^{\circ}$ (shown in panels $(a1)-(a5)$
of Fig. 7) as the initial state of our system and then gradually increase the
twist angle up to $\theta=45^{\circ}$. As can be conveniently seen from the
third and fourth columns of Fig. 7, the state gradually changes from the
topological phase with regions of interlayer phase difference either $0$ or
$\pi$ to another one with regions of interlayer phase differences either
$\pi/2$ or $-\pi/2$. As the twist angle is increased, the vortices forming an
interlayer vortex pair slowly separate from each other. One sees this by
comparing the magnetic profile in panels $(b5)$ and $(c5)$, where the double
peak characteristic of a vortex pair becomes smeared. For twist angles in the
vicinity of $\theta=45^{\circ}$, this culminates in the formation of the
skyrmionic vortex chain state, as we display in panels $(d1)-(d5)$.

Both discussed transitions between the two different topological skyrmionic
vortex states occur through the second-order relocation of vortex cores in
each layer. This is further evidenced in Fig. 8, where we show the energy of
the system as a function of the twist angle for the cases where $\theta$ is
decreased from $45^{\circ}$ to $36^{\circ}$ (blue line, cf. Fig. 6), and
increased from $36^{\circ}$ to $45^{\circ}$ (red line, cf. Fig. 7). Notably,
the two detected characteristic skyrmionic vortex states in the topological
phase of the system cross in energy at a twist angle $\theta^{*}\approx
38^{\circ}$. To emphasize again, the interlayer vortex pairs are energetically
favorable for interlayer twist below this angle, whereas the states containing
skyrmionic vortex chains become favorable for $\theta>\theta^{*}$. Obviously
the exact value of $\theta^{*}$ will depend on the details of the simulation
(size of the unit cell, magnetic field), but we can safely generalize this
result to conclude that skyrmionic vortex chains should be observable in the
larger portion of the twist range where topological phase is expected.

##  IV Conclusion

To summarize, we have analyzed the vortex configurations emerging in a twisted
bilayer composed of superconducting monolayers with $d$-wave pairing -
motivated by prospects of such realizations using e.g. Bi2Sr2CaCu2O8+δ. In
such a system, the phase difference between the superconducting order
parameters in two layers depends on the twist angle $\theta$, with a
topological state with a non-trivial phase difference emerging for a range of
angles around $\theta=45^{\circ}$. In that topological phase, the
superconducting state exhibits broken time-reversal symmetry, giving rise to
skyrmionic vortex configurations with topological charge not equal to zero. We
revealed and characterized those nontrivial vortex states, and discussed their
detectable differences when compared to usual vortex lattice found for twist
angles outside the topological range. Based on those clearly discernible
differences, we argue that direct experimental observation of skyrmionic
vortex states can be used as a smoking gun to detect topological
superconductivity in such systems.

In addition, we showed that the skyrmionic vortex matter also evolves with the
twist angle inside the topological phase. Namely, we have identified two
distinct types of skyrmionic states. At the onset of the topological state,
the system prefers to preserve same phase of the order parameter in two
layers, so the broken reversal symmetry reflects solely in formation of the
interlayer vortex pairs. Each of this pairs carries a unit of topological
charge, and hosts phase difference $\pi$ between the coupled superconducting
layers. As the twist angle is varied towards $45^{\circ}$ and one is deeper in
the topological state, the phase difference of $\pm\pi/2$ becomes
energetically favorable. As a consequence, the interlayer vortex pairs
interconnect into extended closed chains, separating the regions of the sample
with phase difference either $-\pi/2$ or $\pi/2$. Such chains can easily
exhibit lateral extent on the micron scale, and carry topological charge equal
to the number of vortices interconnected in the chain. Once again, we
emphasize that each of the two types of skyrmionic flux objects leaves a clear
signature in the spatial profile of the magnetic field across the system, but
will also host uniquely related local density of states and bound states
detectable by e.g. Scanning Tunneling Microscopy. The calculation of such
states is left as a prospect for further work, being beyond the capability of
the present Ginzburg-Landau analysis (Bogolyubov-deGennes approach is a viable
alternative [29, 30]. Another interesting outlook is to adapt the here-
presented Ginzburg-Landau formalism to the cases of other pairing symmetries
that may arise in the twisted bilayers of present interest, so to classify the
emergent vortex matter according to the symmetries at hand - and thereby
enable conclusive identification of the pairing symmetry in experimental
systems through visualization of the vortex states - complementary to other
existing efforts (see e.g. [31]).

###### Acknowledgements.

This work has been supported by the Research Foundation-Flanders (FWO-
Vlaanderen), Special Research Funds of the University of Antwerp (BOF-UA), and
Brazilian Agency FAPESP (grant numbers 20/03947-2 and 20/10058-0).

## References

  * Kallin and Berlinsky [2016] C. Kallin and J. Berlinsky, Chiral superconductors, Reports on Progress in Physics 79, 054502 (2016).
  * Sigrist and Ueda [1991] M. Sigrist and K. Ueda, Phenomenological theory of unconventional superconductivity, Reviews of Modern physics 63, 239 (1991).
  * Vojta  _et al._ [2000] M. Vojta, Y. Zhang, and S. Sachdev, Quantum phase transitions in d-wave superconductors, Physical review letters 85, 4940 (2000).
  * Stone and Roy [2004] M. Stone and R. Roy, Edge modes, edge currents, and gauge invariance in p x+ i p y superfluids and superconductors, Physical Review B 69, 184511 (2004).
  * Read and Green [2000] N. Read and D. Green, Paired states of fermions in two dimensions with breaking of parity and time-reversal symmetries and the fractional quantum hall effect, Physical Review B 61, 10267 (2000).
  * Garaud and Babaev [2012] J. Garaud and E. Babaev, Skyrmionic state and stable half-quantum vortices in chiral p-wave superconductors, Physical Review B 86, 060514 (2012).
  * Jang  _et al._ [2011] J. Jang, D. Ferguson, V. Vakaryuk, R. Budakian, S. Chung, P. Goldbart, and Y. Maeno, Observation of half-height magnetization steps in sr2ruo4, Science 331, 186 (2011).
  * Zyuzin  _et al._ [2017] A. Zyuzin, J. Garaud, and E. Babaev, Nematic skyrmions in odd-parity superconductors, Physical review letters 119, 167001 (2017).
  * Becerra  _et al._ [2016] V. F. Becerra, E. Sardella, F. Peeters, and M. Milošević, Vortical versus skyrmionic states in mesoscopic p-wave superconductors, Physical Review B 93, 014518 (2016).
  * Milošević and Perali [2015] M. V. Milošević and A. Perali, Emergent phenomena in multicomponent superconductivity: an introduction to the focus issue, Superconductor Science and Technology 28, 060201 (2015).
  * Tanaka [2015] Y. Tanaka, Multicomponent superconductivity based on multiband superconductors, Superconductor Science and Technology 28, 034002 (2015).
  * Babaev and Speight [2005] E. Babaev and M. Speight, Semi-meissner state and neither type-i nor type-ii superconductivity in multicomponent superconductors, Physical Review B 72, 180502 (2005).
  * Lin and Hu [2011] S.-Z. Lin and X. Hu, Vortex states and the phase diagram of a multiple-component ginzburg-landau theory with competing repulsive and attractive vortex interactions, Physical Review B 84, 214505 (2011).
  * Tanaka [2001] Y. Tanaka, Soliton in two-band superconductor, Physical review letters 88, 017002 (2001).
  * Sato and Ando [2017] M. Sato and Y. Ando, Topological superconductors: a review, Reports on Progress in Physics 80, 076501 (2017).
  * Sarma  _et al._ [2015] S. D. Sarma, M. Freedman, and C. Nayak, Majorana zero modes and topological quantum computation, npj Quantum Information 1, 1 (2015).
  * Volovik [1997] G. Volovik, On edge states in superconductors with time inversion symmetry breaking, Journal of Experimental and Theoretical Physics Letters 66, 522 (1997).
  * Volovik [1999] G. Volovik, Fermion zero modes on vortices in chiral superconductors, Journal of Experimental and Theoretical Physics Letters 70, 609 (1999).
  * Ivanov [2001] D. A. Ivanov, Non-abelian statistics of half-quantum vortices in p-wave superconductors, Physical review letters 86, 268 (2001).
  * Can  _et al._ [2021] O. Can, T. Tummuru, R. P. Day, I. Elfimov, A. Damascelli, and M. Franz, High-temperature topological superconductivity in twisted double-layer copper oxides, Nature Physics 17, 519 (2021).
  * Yu  _et al._ [2019] Y. Yu, L. Ma, P. Cai, R. Zhong, C. Ye, J. Shen, G. D. Gu, X. H. Chen, and Y. Zhang, High-temperature superconductivity in monolayer bi2sr2cacu2o8+ $\delta$, Nature 575, 156 (2019).
  * Babaev  _et al._ [2002a] E. Babaev, L. D. Faddeev, and A. J. Niemi, Hidden symmetry and knot solitons in a charged two-condensate bose system, Physical Review B 65, 100512 (2002a).
  * Garaud  _et al._ [2013] J. Garaud, J. Carlström, E. Babaev, and M. Speight, Chiral $\mathbb{C}{P}^{2}$ skyrmions in three-band superconductors, Phys. Rev. B 87, 014507 (2013).
  * Benfenati  _et al._ [2022] A. Benfenati, M. Barkman, and E. Babaev, Demonstration of $\mathbb{C}{P}^{2}$ skyrmions in three-band superconductors by self-consistent solutions to a bogoliubov-de gennes model, arXiv preprint arXiv:2204.05242 (2022).
  * Zhang  _et al._ [2020] L.-F. Zhang, Y.-Y. Zhang, G.-Q. Zha, M. Milošević, and S.-P. Zhou, Skyrmionic chains and lattices in s+ i d superconductors, Physical Review B 101, 064501 (2020).
  * Ren  _et al._ [1995] Y. Ren, J.-H. Xu, and C. Ting, Ginzburg-landau equations and vortex structure of a d x 2- y 2 superconductor, Physical review letters 74, 3680 (1995).
  * Doria  _et al._ [1989] M. M. Doria, J. Gubernatis, and D. Rainer, Virial theorem for ginzburg-landau theories with potential applications to numerical studies of type-ii superconductors, Physical Review B 39, 9573 (1989).
  * Babaev  _et al._ [2002b] E. Babaev, L. D. Faddeev, and A. J. Niemi, Hidden symmetry and knot solitons in a charged two-condensate bose system, Physical Review B 65, 100512 (2002b).
  * Zhang  _et al._ [2016] L.-F. Zhang, V. F. Becerra, L. Covaci, and M. V. Milošević, Electronic properties of emergent topological defects in chiral $p$-wave superconductivity, Phys. Rev. B 94, 024520 (2016).
  * Su and Lin [2018] Y. Su and S.-Z. Lin, Pairing symmetry and spontaneous vortex-antivortex lattice in superconducting twisted-bilayer graphene: Bogoliubov-de gennes approach, Phys. Rev. B 98, 195101 (2018).
  * Wu and Das Sarma [2019] F. Wu and S. Das Sarma, Identification of superconducting pairing symmetry in twisted bilayer graphene using in-plane magnetic field and strain, Phys. Rev. B 99, 220507 (2019).


# Can supermassive black holes triggered by general QCD axion bubbles?

Hai-Jun Li<EMAIL_ADDRESS>Center for Advanced Quantum Studies,
Department of Physics, Beijing Normal University, Beijing 100875, China  Ying-
Quan Peng<EMAIL_ADDRESS>Center for Advanced Quantum Studies,
Department of Physics, Beijing Normal University, Beijing 100875, China  Wei
Chao<EMAIL_ADDRESS>Center for Advanced Quantum Studies, Department of
Physics, Beijing Normal University, Beijing 100875, China  Yu-Feng Zhou
<EMAIL_ADDRESS>Key Laboratory of Theoretical Physics, Institute of
Theoretical Physics, Chinese Academy of Sciences, Beijing 100190, China
School of Physical Sciences, University of Chinese Academy of Sciences,
Beijing 100049, China  School of Fundamental Physics and Mathematical
Sciences, Hangzhou Institute for Advanced Study, UCAS, Hangzhou 310024, China
International Centre for Theoretical Physics Asia-Pacific, Beijing/Hangzhou,
China

###### Abstract

The supermassive black holes (SMBHs) are ubiquitous in the center of galaxies,
although the origin of their massive seeds is still unknown. In this paper, we
investigate the SMBHs formation from the general QCD axion bubbles. In this
case, the primordial black holes (PBHs) are considered as the seeds of SMBHs,
which are generated from the QCD axion bubbles due to an additional Peccei-
Quinn (PQ) symmetry breaking after inflation. The QCD axion bubbles are formed
when the QCD axion starts to oscillate during the QCD phase transition (PT).
We consider the general case in which the QCD axion bubbles are formed with
the bubble effective angle $\theta_{\rm eff}\in(0,\,\pi]$, leading to the
minimum PBH mass $\sim\mathcal{O}(10^{4}-10^{7})M_{\odot}$ with $\theta_{\rm
eff}\sim\pi$ to $\pi/3$ for the axion scale
$f_{a}\sim\mathcal{O}(10^{16})\,\rm GeV$. The PBHs at this mass region may be
the seeds of SMBHs.

††preprint: BNU-23-043

##  I Introduction

The supermassive black holes (SMBHs) with the mass
$\sim\mathcal{O}(10^{6}-10^{9})M_{\odot}$ are ubiquitous in the center of
galaxies Kormendy and Richstone (1995); Richstone  _et al._ (1998); Kormendy
and Ho (2013), where $M_{\odot}$ is the solar mass. They are also considered
as the central engines of active galactic nuclei (AGN) and quasars Lynden-Bell
(1969); Ferrarese and Merritt (2000). However, the origin of SMBHs is not
clear. They may be formed from the stellar black holes (BHs) with the mass
$\sim\mathcal{O}(10)M_{\odot}$ through accretion and mergers Page and Thorne
(1974); Buonanno  _et al._ (2008). However, it is difficult to account for the
SMBHs at the high redshift $z\sim 7$ Mortlock  _et al._ (2011); De Rosa  _et
al._ (2014); Banados  _et al._ (2018) since the stellar BHs have no time to
grow into SMBHs. The another scenario is considering the primordial black
holes (PBHs) with the mass $\sim\mathcal{O}(10^{4}-10^{5})M_{\odot}$ as the
seeds of SMBHs Salpeter (1964); Duechting (2004). In this case, the PBHs could
subsequently grow up to $\sim\mathcal{O}(10^{9})M_{\odot}$ due to an efficient
accretion of matter on the massive seeds and mergings Choquette  _et al._
(2019); Dolgov and Porey (2019). (See $\rm e.g.$ Refs. Volonteri (2010);
Dolgov (2018) for reviews of the SMBHs formation.)

The PBHs have recently gathered significant interest as the candidates of dark
matter (DM) Bird  _et al._ (2016); Green and Kavanagh (2021); Carr and Kuhnel
(2020, 2022); Li (2022); Jho  _et al._ (2022); Agashe  _et al._ (2022); Xie
(2023); Gelmini  _et al._ (2023). They can be formed due to the collapse of
large density fluctuations produced during inflation Carr (1975); Ivanov  _et
al._ (1994); Garcia-Bellido  _et al._ (1996); Alabidi and Kohri (2009), the
cosmological phase transitions Kodama  _et al._ (1982); Hawking  _et al._
(1982); Baker  _et al._ (2021), and the collapse of topological defects
Hawking (1989); Polnarev and Zembowicz (1991); MacGibbon  _et al._ (1998);
Rubin  _et al._ (2001); Deng  _et al._ (2017); Ferrer  _et al._ (2019) and
false vacuum bubbles Deng and Vilenkin (2017); Deng (2020); Maeso  _et al._
(2022), etc. The another scenario is considering the PBHs formation from the
QCD axion bubbles in the early Universe, leading to the minimum PBH mass
$\sim\mathcal{O}(10^{4})M_{\odot}$, which was first proposed in Ref. Kitajima
and Takahashi (2020) and further studied by Ref. Li (2023). This scenario is
similar to the PBHs formation from the baryon bubbles in the Affleck-Dine (AD)
baryogenesis Dolgov  _et al._ (2009); Hasegawa and Kawasaki (2019); Kawasaki
and Murai (2019). The QCD axion Peccei and Quinn (1977a, b); Weinberg (1978);
Wilczek (1978) is a light pseudo-Nambu-Goldstone boson with a spontaneously
broken $U(1)$ Peccei-Quinn (PQ) symmetry. It is also the potential cold DM
candidate, which can be non-thermally generated through misalignment mechanism
in the early Universe Preskill  _et al._ (1983); Abbott and Sikivie (1983);
Dine and Fischler (1983).

The QCD axion bubbles can be generated due to an additional PQ symmetry
breaking after inflation with the multiple vacua, and are formed when the QCD
axion starts to oscillate during the QCD phase transition (PT) Kitajima and
Takahashi (2020); Li (2023). The axion acquires a light mass after inflation,
then starts to oscillate when this mass is comparable to the Hubble parameter
and settles down into different potential minima depending on its initial
position. Since the additional PQ symmetry is supposed to be temporarily
broken, the axion potential will disappear before the QCD PT. Therefore, the
misalignment mechanism is still valid to calculate the final QCD axion
abundance. However, the initial misalignment angle will be split into
different values by the multiple vacua. During the QCD PT, the QCD axions
start to oscillate with these different initial angles, in which the high
density axion bubbles can be produced with the large initial values. An
interesting phenomenon of the QCD axion bubbles is the PBHs formation, which
can be formed when the axions dominate the radiation in the bubbles.

In this paper, we focus our attention on the PBHs triggered by QCD axion
bubbles as the seeds of SMBHs. Here we assume that the axion acquires a light
mass due to an additional PQ symmetry breaking after inflation with a large
axion scale $f_{a}\sim\mathcal{O}(10^{16})\,\rm GeV$, and then settles down
into two potential minima, $\phi_{\rm min}^{0}$ and $\phi_{\rm min}^{1}$. We
consider the general case in which the QCD axion bubbles are formed with the
bubble effective angle $\theta_{\rm eff}\in(0,\,\pi]$. The $\phi_{\rm
min}^{0}$ accounts for the DM abundance with the initial misalignment angle
$\sim 0$, while the $\phi_{\rm min}^{1}$ forms the high energy density QCD
axion bubbles with the angle $\sim\theta_{\rm eff}$. The PBHs triggered by QCD
axion bubbles are formed after the QCD PT, which leads to the minimum PBH mass
$\sim\mathcal{O}(10^{4}-10^{7})M_{\odot}$ for
$f_{a}\sim\mathcal{O}(10^{16})\,\rm GeV$. Compared with the critical value of
the $\sim\mathcal{O}(10^{9})M_{\odot}$ SMBHs seeds, the PBHs at this mass
region may be the seeds of SMBHs.

The paper is structured as follows. In Sec. II, we introduce the general QCD
axion bubbles with the bubble effective angle $\theta_{\rm eff}\in(0,\,\pi]$.
In Sec. III, we investigate the PBHs as the seeds of SMBHs, which are
generated from the QCD axion bubbles. Finally, the conclusion is given in Sec.
IV.

##  II General QCD axion bubbles

In this section, we introduce the general QCD axion bubbles. We first briefly
review the misalignment mechanism. The QCD axion effective potential is given
by

$\displaystyle V_{\rm
QCD}(\phi)=m_{a}^{2}(T)f_{a}^{2}\left[1-\cos\left(\frac{\phi}{f_{a}}\right)\right]\,,$
(1)

with the axion field $\phi$, the axion decay constant $f_{a}$, the axion angle
$\theta=\phi/f_{a}$, and the temperature-dependent axion mass Borsanyi  _et
al._ (2016)

$\displaystyle\begin{aligned}
m_{a}(T)\simeq\begin{cases}m_{a,0}\left(\dfrac{T}{T_{\rm
QCD}}\right)^{-4.08}\,,&T\geq T_{\rm QCD}\\\ m_{a,0}\,,&T<T_{\rm
QCD}\end{cases}\end{aligned}$ (2)

where $m_{a,0}\simeq 5.70(7)\,{\mu\rm eV}(10^{12}\,{\rm GeV}/f_{a})$ is the
zero-temperature axion mass Grilli di Cortona  _et al._ (2016), $T$ is the
cosmic temperature, and $T_{\rm QCD}\simeq 150\,\rm MeV$ is the QCD scale.

During the QCD PT, the QCD axion starts to oscillate when its mass $m_{a}(T)$
is comparable to the Hubble parameter $H(T)$ with the oscillation temperature
$T_{a}$. The axion energy density at present can be described by

$\displaystyle\rho_{a}(T_{0})=\frac{m_{a,0}m_{a}(T_{a})s(T_{0})}{2s(T_{a})}f_{a}^{2}\left\langle\theta_{i}^{2}f(\theta_{i})\right\rangle\chi\,,$
(3)

with the entropy density $s(T)$, the present CMB temperature $T_{0}$, the
initial misalignment angle $\theta_{i}$, the numerical factor $\chi\simeq
1.44$ Turner (1986), and the anharmonic factor $f(\theta_{i})$ is given by
Lyth (1992); Visinelli and Gondolo (2009)

$\displaystyle
f(\theta_{i})\simeq\left[\ln\left(\frac{e}{1-\theta_{i}^{2}/\pi^{2}}\right)\right]^{1.16}\,.$
(4)

Then we derive the present QCD axion abundance

$\displaystyle\begin{aligned} \Omega_{a}h^{2}&\simeq
0.14\left(\frac{g_{*}(T_{a})}{61.75}\right)^{-0.42}\left(\frac{g_{*s}(T_{0})}{3.94}\right)\\\
&\times\left(\frac{f_{a}}{10^{12}\rm\,GeV}\right)^{1.16}\left\langle\theta_{i}^{2}f(\theta_{i})\right\rangle\,,\end{aligned}$
(5)

where $g_{*}(T)$ and $g_{*s}(T)$ are the numbers of effective degrees of
freedom of the energy density and the entropy density, respectively, and
$h\simeq 0.68$ is the reduced Hubble constant. In order to account for the
observed cold DM abundance ($\Omega_{\rm DM}h^{2}\simeq 0.12$ Aghanim  _et
al._ (2020)), the initial misalignment angle should be
$\theta_{i}\sim\mathcal{O}(1)$ for the scale
$f_{a}\sim\mathcal{O}(10^{12})\,\rm GeV$.

The QCD axion bubbles can be generated due to an additional PQ symmetry
breaking after inflation. There are many mechanisms for this additional PQ
symmetry breaking in the early Universe Nomura  _et al._ (2016); Kawasaki  _et
al._ (2016, 2018); Chiba  _et al._ (2004); Takahashi and Yamaguchi (2004);
Higaki  _et al._ (2014); Choi  _et al._ (1997); Banks and Dine (1997); Jeong
and Takahashi (2013); Takahashi and Yamada (2015). Then the axion field $\phi$
acquires the effective potential, which can be approximated as

$\displaystyle V_{\cancel{\rm
PQ}}(\phi)\simeq\frac{1}{2}m_{a}^{2}(\phi)\left(\phi-\phi_{\rm
min}^{n}\right)^{2}\,,$ (6)

where $m_{a}(\phi)$ is the axion effective mass, $\phi_{\rm min}^{n}$ is the
potential minimum, and $n$ is an integer. Note that the additional PQ symmetry
is assumed to be temporarily broken with the large scale
$f_{a}\sim\mathcal{O}(10^{16})\,\rm GeV$, which gives the axion a light
effective mass and forms the multiple vacua. Then the axion starts to
oscillate and settles down into the minimum $\phi_{\rm min}^{n}$ of
$V_{\cancel{\rm PQ}}(\phi)$.

The QCD axion bubbles are generated when the conventional axion potential
$V_{\rm QCD}(\phi)$ arises during the QCD PT. Therefore, one of the minimum
$\phi_{\rm min}^{0}$ of $V_{\cancel{\rm PQ}}(\phi)$ should be near the minimum
of $V_{\rm QCD}(\phi)$, ensuring the correct cold DM abundance with the
effective initial angle $\theta_{i,0}$. In addition, if the initial value is
greater than a critical value $\phi_{\rm crit}$, the axion will be stabilized
at the another minimum $\phi_{\rm min}^{1}$ with $\theta_{i,1}$. Here we
suppose that the QCD axion bubbles are produced at $\phi_{\rm min}^{1}$, $\rm
i.e.$, $\theta_{i,1}>\theta_{i,0}$, we define the axion bubble effective angle

$\displaystyle\boxed{\theta_{\rm eff}\in(0,\,\pi]\,.}$ (7)

This can be realized with the bubble effective angle, $\rm e.g.$, $\theta_{\rm
eff}=\pi$ Kitajima and Takahashi (2020) and $2\pi/3$ Li (2023). In Ref.
Kitajima and Takahashi (2020), they showed the Witten effect as an example of
the additional PQ symmetry breaking to produce the QCD axion bubbles with
$\theta_{\rm eff}=\pi$. In Ref. Li (2023), they discussed the axion bubbles
formation from the hidden $SU(N)$ gauge symmetry breaking with $\theta_{\rm
eff}=2\pi/3$. It should be emphasized that the lower limit of Eq. (7) is just
a hypothetical scope of this angle, not means all values in this region are
reasonable, which is related to the specific mechanism of the additional PQ
symmetry breaking. We do know the bubble effective angle cannot be small as
$\sim 0$, in this case, we have $\theta_{i,1}\simeq\theta_{i,0}$, $\rm i.e.$,
no QCD axion bubbles are produced. However, it is difficult to define a clear
lower limit of $\theta_{\rm eff}$ in this work. For discussion purposes, we
take a typical lower limit $\theta_{\rm eff}=\pi/3$ in what follows. Then the
effective initial misalignment angle are given by

$\displaystyle\theta_{i,n}=\begin{cases}0-\theta_{i}\,,&n=0\\\ \theta_{\rm
eff}-\theta_{i}\,,&n=1\end{cases}$ (8)

corresponding to the potential minima $\phi_{\rm min}^{0}$ and $\phi_{\rm
min}^{1}$, respectively. Note that $\theta_{i}$ in Eq. (8) is a small initial
misalignment angle

$\displaystyle\begin{aligned} \theta_{i}&\simeq 4.29\times
10^{-3}\left(\frac{g_{*}(T_{a})}{61.75}\right)^{0.21}\left(\frac{g_{*s}(T_{0})}{3.94}\right)^{1/2}\\\
&\times\left(\frac{f_{a}}{10^{16}\,\rm GeV}\right)^{-0.58}\,,\end{aligned}$
(9)

for $f_{a}\sim\mathcal{O}(10^{16})\,\rm GeV$, which is derived from Eq. (5).

Figure 1: The temperature $T_{B}$ as a function of the axion bubble effective
angle $\theta_{\rm eff}$. The red and blue lines represent the scales
$f_{a}=1\times 10^{16}\,\rm GeV$ and $1\times 10^{17}\,\rm GeV$, respectively.

When the axions dominate the radiation in the bubbles, the cosmic background
temperature is defined as $T_{B}$ Kitajima and Takahashi (2020). The axion
energy density at $T_{B}$ is equal to the radiation energy density,
$\rho_{a}(T_{B})=\rho_{R}(T_{B})$, where $\rho_{R}$ is the radiation energy
density. Since the axion energy density inside the bubbles at $T_{B}$ can be
described by

$\displaystyle\begin{aligned}
\rho_{a}(T_{B})&=\frac{m_{a,0}m_{a}(T_{a})s(T_{B})}{2s(T_{a})}\\\ &\times
f_{a}^{2}\left\langle\left(\theta_{\rm eff}-\theta_{i}\right)^{2}f(\theta_{\rm
eff}-\theta_{i})\right\rangle\chi\,,\end{aligned}$ (10)

then we can derive the temperature $T_{B}$

$\displaystyle\begin{aligned} T_{B}&\simeq 2.13\times 10^{-2}\,{\rm MeV}\\\
&\times\left(\frac{g_{*}(T_{a})}{61.75}\right)^{-0.42}\left(\frac{f_{a}}{10^{16}\,\rm
GeV}\right)^{1.16}\\\ &\times\left\langle\left(\theta_{\rm
eff}-\theta_{i}\right)^{2}f(\theta_{\rm
eff}-\theta_{i})\right\rangle\chi\,.\end{aligned}$ (11)

Using Eq. (9), we show $T_{B}$ as a function of the axion bubble effective
angle $\theta_{\rm eff}$ in Fig. 1. Two typical values of the scales
$f_{a}=1\times 10^{16}\,\rm GeV$ (red) and $1\times 10^{17}\,\rm GeV$ (blue)
are selected for comparisons. We find that the temperature
$T_{B}\sim\mathcal{O}(0.1-1)\,\rm MeV$ decreases slowly as the effective angle
$\theta_{\rm eff}$ decreases, and varies most rapidly when the angle
approaches a small value.

Figure 2: The minimum PBH mass $M_{\rm PBH}^{\rm min}$ (in the solar mass
$M_{\odot}$) as a function of the angle $\theta_{\rm eff}$. The red and blue
lines represent $f_{a}=1\times 10^{16}\,\rm GeV$ and $1\times 10^{17}\,\rm
GeV$, respectively. The gray shadow region represents the PBH mass region
$\sim\mathcal{O}(10^{4}-10^{7})M_{\odot}$.

The QCD axion bubbles abundance is related to the inflationary fluctuations.
The volume fraction of the bubbles can be described by Kitajima and Takahashi
(2020)

$\displaystyle\frac{{\rm d}\beta}{{\rm d}\ln k}\simeq\frac{\phi_{\rm
crit}-\phi_{i}}{2}P(\ln(k/k_{*}),\phi_{\rm crit})\,,$ (12)

where $\phi_{\rm crit}$ is the critical value, $k$ is the wave number,
$k_{*}=0.002\,\rm Mpc^{-1}$, and $P(\ln(k/k_{*}),\phi)$ is the probability
density function

$\displaystyle\begin{aligned}
P(\ln(k/k_{*}),\phi)&=\frac{1}{\sqrt{2\pi}\sigma(\ln(k/k_{*}))}\\\
&\times\exp{\left(-\frac{\left(\phi-\phi_{i}\right)^{2}}{2\sigma^{2}(\ln(k/k_{*}))}\right)}\,,\end{aligned}$
(13)

where $\sigma(\ln(k/k_{*}))=H_{\rm inf}\sqrt{\ln(k/k_{*})}/(2\pi)$ is the
variance, $P(\ln(k/k_{*}),\phi)=\delta(\phi-\phi_{i})$ is the initial
condition, and $H_{\rm inf}\sim\mathcal{O}(10^{13})\,\rm GeV$ is the Hubble
parameter during inflation. Here we take $\phi_{\rm crit}-\phi_{i}=4.5H_{\rm
inf}$ as a benchmark.

##  III SMBHs from the general QCD axion bubbles

The astronomical observations indicate that SMBHs are ubiquitous in the center
of galaxies with the mass $\sim 10^{9}M_{\odot}$ at $z\sim 7$ ($\sim 0.76\,\rm
Gyr$) Mortlock  _et al._ (2011); De Rosa  _et al._ (2014); Banados  _et al._
(2018). However, the origin of such BHs is still a mystery. The mass of an
accreting BH with the time $t$ is given by Salpeter (1964)

$\displaystyle
M(t)=M_{0}\exp{\left(\frac{1-\epsilon_{r}}{\epsilon_{r}}\frac{t}{t_{E}}\right)}\,,$
(14)

where $M_{0}$ is the seed BH mass, $\epsilon_{r}\simeq 0.1$ is the radiative
efficiency, and $t_{E}\simeq 0.45\,\rm Gyr$. In this case, a seed BH with the
mass $M_{0}\sim\mathcal{O}(10^{2}-10^{5})M_{\odot}$ would take at least $\sim
0.5\,\rm Gyr$ to grow up to a $\sim 10^{9}M_{\odot}$ SMBH Choquette  _et al._
(2019). Therefore, the SMBHs must either have the heavy seeds or have a
primordial origin.

Figure 3: The distributions of the PBH fractional abundance $\log(f_{\rm
PBH})$ on the “$\log(M_{\rm PBH}/M_{\odot})-\theta_{\rm eff}$” plane for
$f_{a}=1\times 10^{16}\,\rm GeV$ (left) and $1\times 10^{17}\,\rm GeV$
(right). The shadow regions represent the values of $\log(f_{\rm PBH})$. The
blank regions represent the limits set by $M_{\rm PBH}^{\rm min}$.

Here we investigate the PBHs as the seeds of SMBHs, which are produced from
the above general QCD axion bubbles after the QCD PT. The initial PBH mass at
the formation temperature $T_{f}$ is Carr  _et al._ (2010)

$\displaystyle M_{\rm
PBH}=\frac{4\pi}{3}\frac{\gamma\rho_{R}(T_{f})}{H^{3}(T_{f})}\,,$ (15)

where $\gamma\simeq 0.2$ is the gravitational collapse factor Carr (1975), and
$H(T_{f})$ is the Hubble parameter at $T_{f}$. Considering the PBHs triggered
by QCD axion bubbles, they will be produced when the axions dominate the
radiation in the bubbles, $\rm i.e.$, $T_{f}\lesssim T_{B}$. Using Eqs. (11)
and (15), we can derive the minimum PBH mass

$\displaystyle\begin{aligned} \frac{M_{\rm PBH}^{\rm min}}{M_{\odot}}&\simeq
6.58\times
10^{7}\left(\frac{\gamma}{0.2}\right)\left(\frac{g_{*}(T_{a})}{61.75}\right)^{0.84}\\\
&\times\left(\frac{g_{*}(T_{f})}{10.75}\right)^{-1/2}\left(\frac{f_{a}}{10^{16}\,\rm
GeV}\right)^{-2.33}\\\ &\times\left(\left\langle\left(\theta_{\rm
eff}-\theta_{i}\right)^{2}f(\theta_{\rm
eff}-\theta_{i})\right\rangle\chi\right)^{-2}\,.\end{aligned}$ (16)

We show $M_{\rm PBH}^{\rm min}$ as a function of the bubble effective angle
$\theta_{\rm eff}$ in Fig. 2. The red and blue lines represent $f_{a}=1\times
10^{16}\,\rm GeV$ and $1\times 10^{17}\,\rm GeV$, respectively. For
$f_{a}=1\times 10^{16}\,\rm GeV$, we find the minimum PBH mass
$\sim\mathcal{O}(10^{4}-10^{7})M_{\odot}$ with the angle $\theta_{\rm
eff}\sim\pi$ to $\pi/3$, which is shown with the gray shadow region. This mass
region would be smaller for $f_{a}=1\times 10^{17}\,\rm GeV$.

Here we discuss the PBH abundance. The energy density of PBH at the
temperature $T_{B}$ is given by

$\displaystyle\rho_{\rm PBH}(T_{B})=\frac{3}{4}T_{B}s(T_{B})\frac{{\rm
d}\beta}{{\rm d}\ln k}\,.$ (17)

Then the fractional abundance of PBH at present can be described by Li (2023)

$\displaystyle f_{\rm PBH}=\frac{3}{4}T_{B}s(T_{0})\frac{{\rm d}\beta}{{\rm
d}\ln k}\frac{1}{\rho_{c}}\frac{1}{\Omega_{\rm DM}}\,,$ (18)

where $\Omega_{\rm DM}\simeq 0.268$ is the total cold DM abundance Aghanim
_et al._ (2020), and $\rho_{c}$ is the critical energy density. Using Eq.
(12), we show the distributions of $\log(f_{\rm PBH})$ on the parameters
“$\log(M_{\rm PBH}/M_{\odot})-\theta_{\rm eff}$” plane in Fig. 3. The left and
right panels correspond to the scales $f_{a}=1\times 10^{16}\,\rm GeV$ and
$1\times 10^{17}\,\rm GeV$, respectively. For the same $\theta_{\rm eff}$, the
value of $f_{\rm PBH}$ decreases as the PBH mass $M_{\rm PBH}$ increases. Note
that the blank regions in the plots represent the limits set by $M_{\rm
PBH}^{\rm min}$.

Considering the PBHs as the seeds of SMBHs, the numerical simulations show
that PBHs with the mass $\sim\mathcal{O}(10^{4}-10^{5})M_{\odot}$ could
subsequently grow up to the $\sim\mathcal{O}(10^{9})M_{\odot}$ SMBHs Rosas-
Guevara  _et al._ (2016). In this case, the PBHs mass larger than a critical
value $M_{c}$ are considered as the seeds of SMBHs

$\displaystyle M_{c}=\left(10^{4}-10^{5}\right)M_{\odot}\,.$ (19)

Therefore, the seed PBH mass in the QCD axion bubbles scenario is

$\displaystyle M_{s}={\rm Max}\left[M_{\rm PBH}^{\rm min},\,M_{c}\right]\,,$
(20)

where $M_{\rm PBH}^{\rm min}$ is given by Eq. (16). We show the fractional
abundance of PBH as a function of the PBH mass $M_{\rm PBH}$ for $\theta_{\rm
eff}=\pi$ in Fig. 4. The red and blue lines represent $f_{a}=1\times
10^{16}\,\rm GeV$ and $1\times 10^{17}\,\rm GeV$, respectively. The gray cut-
offs represent the mass $M_{\rm PBH}^{\rm min}$. The red and blue cut-offs
represent the mass $M_{s}$ [Eq. (20)] with $M_{c}=10^{4}M_{\odot}$. We find
that the values of $M_{\rm PBH}^{\rm min}$ are smaller than that of $M_{c}$
for $\theta_{\rm eff}=\pi$. We also plot the case for $\theta_{\rm eff}=\pi/2$
in Fig. 4, which shows a different result. The PBHs at the mass region
$\sim\mathcal{O}(10^{4}-10^{7})M_{\odot}$ may be the seeds of SMBHs.

Figure 4: The PBH fractional abundance $f_{\rm PBH}$ as a function of the PBH
mass for $\theta_{\rm eff}=\pi$ (left) and $\pi/2$ (right). The red and blue
lines represent the scales $f_{a}=1\times 10^{16}\,\rm GeV$ and $1\times
10^{17}\,\rm GeV$, respectively. The gray cut-offs represent the minimum PBH
mass $M_{\rm PBH}^{\rm min}$. The red and blue cut-offs represent the values
of $M_{s}$ with $M_{c}=10^{4}M_{\odot}$.

Finally, we comment on the SMBHs formation from the baryon bubbles in the AD
baryogenesis Kawasaki and Murai (2019). The SMBHs formation in the QCD axion
bubbles scenario has some similarities with their mechanism. They considered
the high density baryon bubbles formation from the modified AD mechanism with
the Hubble induced mass and the finite temperature effect. The PBHs will be
formed when the baryon bubbles enter the horizontal plane with the
sufficiently large densities. They also showed the PBHs produced by that
mechanism may have a reasonable mass as the seeds of SMBHs.

##  IV Conclusion

In this paper, we have investigated the SMBHs formation from the general QCD
axion bubbles. The main idea is considering the PBHs triggered by QCD axion
bubbles as the seeds of SMBHs. We first introduce the general QCD axion
bubbles, which can be generated due to an additional PQ symmetry breaking
after inflation and are formed when the QCD axion starts to oscillate during
the QCD PT. In this work, we do not discuss the specific mechanism of the
additional PQ symmetry breaking. We consider the general case in which the QCD
axion bubbles are formed with the bubble effective angle $\theta_{\rm
eff}\in(0,\,\pi]$. The PBHs triggered by axion bubbles will be formed when the
axions dominate the radiation in the bubbles, which leads to a minimum PBH
mass $M_{\rm PBH}^{\rm min}$. We show the distributions of the temperature
$T_{B}$ with different $\theta_{\rm eff}$ for the scales $f_{a}=1\times
10^{16}\,\rm GeV$ and $1\times 10^{17}\,\rm GeV$. Then we find the minimum PBH
mass $M_{\rm PBH}^{\rm min}\sim\mathcal{O}(10^{4}-10^{7})M_{\odot}$ with the
angle $\theta_{\rm eff}\sim\pi$ to $\pi/3$ for $f_{a}=1\times 10^{16}\,\rm
GeV$. The PBH fractional abundance distributions on the “$\log(M_{\rm
PBH}/M_{\odot})-\theta_{\rm eff}$” plane are also shown, which are strongly
limited by $M_{\rm PBH}^{\rm min}$.

In order to explain the SMBHs in the galactic centers with the mass $\sim
10^{9}M_{\odot}$ at $z\sim 7$ with the PBH scenario, the PBH mass should be
larger than a critical value $\sim\mathcal{O}(10^{4}-10^{5})M_{\odot}$.
Combined with the limits set by $M_{\rm PBH}^{\rm min}$, we show the PBH
fractional abundance distributions for the effective angle $\theta_{\rm
eff}=\pi$ and $\pi/2$ for $f_{a}=1\times 10^{16}\,\rm GeV$ and $1\times
10^{17}\,\rm GeV$. We find that the PBHs triggered by QCD axion bubbles with
the mass $\sim\mathcal{O}(10^{4}-10^{7})M_{\odot}$ may be the seeds of SMBHs.
The SMBHs formation in the QCD axion bubbles scenario has some similarities
with that in the baryon bubbles scenario in the AD baryogenesis.

## Acknowledgments

WC is supported by the National Natural Science Foundation of China (NSFC)
(Grants No. 11775025 and No. 12175027). YFZ is supported by the National Key
R&D Program of China (Grant No. 2017YFA0402204), the CAS Project for Young
Scientists in Basic Research YSBR-006, and the NSFC (Grants No. 11821505, No.
11825506, and No. 12047503).

## References

  * Kormendy and Richstone (1995) John Kormendy and Douglas Richstone, “Inward bound: The Search for supermassive black holes in galactic nuclei,” Ann. Rev. Astron. Astrophys. 33, 581 (1995).
  * Richstone  _et al._ (1998) D. Richstone _et al._ , “Supermassive black holes and the evolution of galaxies,” Nature 395, A14–A19 (1998), arXiv:astro-ph/9810378 .
  * Kormendy and Ho (2013) John Kormendy and Luis C. Ho, “Coevolution (Or Not) of Supermassive Black Holes and Host Galaxies,” Ann. Rev. Astron. Astrophys. 51, 511–653 (2013), arXiv:1304.7762 [astro-ph.CO] .
  * Lynden-Bell (1969) D. Lynden-Bell, “Galactic nuclei as collapsed old quasars,” Nature 223, 690 (1969).
  * Ferrarese and Merritt (2000) Laura Ferrarese and David Merritt, “A Fundamental relation between supermassive black holes and their host galaxies,” Astrophys. J. Lett. 539, L9 (2000), arXiv:astro-ph/0006053 .
  * Page and Thorne (1974) Don N. Page and Kip S. Thorne, “Disk-Accretion onto a Black Hole. Time-Averaged Structure of Accretion Disk,” Astrophys. J. 191, 499–506 (1974).
  * Buonanno  _et al._ (2008) Alessandra Buonanno, Lawrence E. Kidder, and Luis Lehner, “Estimating the final spin of a binary black hole coalescence,” Phys. Rev. D 77, 026004 (2008), arXiv:0709.3839 [astro-ph] .
  * Mortlock  _et al._ (2011) Daniel J. Mortlock _et al._ , “A luminous quasar at a redshift of z = 7.085,” Nature 474, 616 (2011), arXiv:1106.6088 [astro-ph.CO] .
  * De Rosa  _et al._ (2014) G. De Rosa _et al._ , “Black hole mass estimates and emission-line properties of a sample of redshift z$>$6.5 quasars,” Astrophys. J. 790, 145 (2014), arXiv:1311.3260 [astro-ph.CO] .
  * Banados  _et al._ (2018) Eduardo Banados _et al._ , “An 800-million-solar-mass black hole in a significantly neutral Universe at redshift 7.5,” Nature 553, 473–476 (2018), arXiv:1712.01860 [astro-ph.GA] .
  * Salpeter (1964) E. E. Salpeter, “Accretion of Interstellar Matter by Massive Objects,” Astrophys. J. 140, 796–800 (1964).
  * Duechting (2004) Norbert Duechting, “Supermassive black holes from primordial black hole seeds,” Phys. Rev. D 70, 064015 (2004), arXiv:astro-ph/0406260 .
  * Choquette  _et al._ (2019) Jeremie Choquette, James M. Cline, and Jonathan M. Cornell, “Early formation of supermassive black holes via dark matter self-interactions,” JCAP 07, 036 (2019), arXiv:1812.05088 [astro-ph.CO] .
  * Dolgov and Porey (2019) A. D. Dolgov and S. Porey, “Massive primordial black holes in contemporary universe,” Bulg. Astron. J. 34, 2021 (2019), arXiv:1905.10972 [astro-ph.CO] .
  * Volonteri (2010) Marta Volonteri, “Formation of Supermassive Black Holes,” Astron. Astrophys. Rev. 18, 279–315 (2010), arXiv:1003.4404 [astro-ph.CO] .
  * Dolgov (2018) A. D. Dolgov, “Massive and supermassive black holes in the contemporary and early Universe and problems in cosmology and astrophysics,” Usp. Fiz. Nauk 188, 121–142 (2018), arXiv:1705.06859 [astro-ph.CO] .
  * Bird  _et al._ (2016) Simeon Bird, Ilias Cholis, Julian B. Muñoz, Yacine Ali-Haïmoud, Marc Kamionkowski, Ely D. Kovetz, Alvise Raccanelli, and Adam G. Riess, “Did LIGO detect dark matter?” Phys. Rev. Lett. 116, 201301 (2016), arXiv:1603.00464 [astro-ph.CO] .
  * Green and Kavanagh (2021) Anne M. Green and Bradley J. Kavanagh, “Primordial Black Holes as a dark matter candidate,” J. Phys. G 48, 043001 (2021), arXiv:2007.10722 [astro-ph.CO] .
  * Carr and Kuhnel (2020) Bernard Carr and Florian Kuhnel, “Primordial Black Holes as Dark Matter: Recent Developments,” Ann. Rev. Nucl. Part. Sci. 70, 355–394 (2020), arXiv:2006.02838 [astro-ph.CO] .
  * Carr and Kuhnel (2022) Bernard Carr and Florian Kuhnel, “Primordial black holes as dark matter candidates,” SciPost Phys. Lect. Notes 48, 1 (2022), arXiv:2110.02821 [astro-ph.CO] .
  * Li (2022) Hai-Jun Li, “Primordial black holes induced stochastic axion-photon oscillations in primordial magnetic field,” JCAP 11, 045 (2022), arXiv:2208.04605 [astro-ph.CO] .
  * Jho  _et al._ (2022) Yongsoo Jho, Tae-Geun Kim, Jong-Chul Park, Seong Chan Park, and Yeji Park, “Axions from Primordial Black Holes,” (2022), arXiv:2212.11977 [hep-ph] .
  * Agashe  _et al._ (2022) Kaustubh Agashe, Jae Hyeok Chang, Steven J. Clark, Bhaskar Dutta, Yuhsin Tsai, and Tao Xu, “Detecting Axion-Like Particles with Primordial Black Holes,” (2022), arXiv:2212.11980 [hep-ph] .
  * Xie (2023) Ke-Pan Xie, “Pinning down the primordial black hole formation mechanism with gamma-rays and gravitational waves,” (2023), arXiv:2301.02352 [astro-ph.CO] .
  * Gelmini  _et al._ (2023) Graciela B. Gelmini, Jonah Hyman, Anna Simpson, and Edoardo Vitagliano, “Primordial black hole dark matter from catastrogenesis with unstable pseudo-Goldstone bosons,” (2023), arXiv:2303.14107 [hep-ph] .
  * Carr (1975) Bernard J. Carr, “The Primordial black hole mass spectrum,” Astrophys. J. 201, 1–19 (1975).
  * Ivanov  _et al._ (1994) P. Ivanov, P. Naselsky, and I. Novikov, “Inflation and primordial black holes as dark matter,” Phys. Rev. D 50, 7173–7178 (1994).
  * Garcia-Bellido  _et al._ (1996) Juan Garcia-Bellido, Andrei D. Linde, and David Wands, “Density perturbations and black hole formation in hybrid inflation,” Phys. Rev. D 54, 6040–6058 (1996), arXiv:astro-ph/9605094 .
  * Alabidi and Kohri (2009) Laila Alabidi and Kazunori Kohri, “Generating Primordial Black Holes Via Hilltop-Type Inflation Models,” Phys. Rev. D 80, 063511 (2009), arXiv:0906.1398 [astro-ph.CO] .
  * Kodama  _et al._ (1982) Hideo Kodama, Misao Sasaki, and Katsuhiko Sato, “Abundance of Primordial Holes Produced by Cosmological First Order Phase Transition,” Prog. Theor. Phys. 68, 1979 (1982).
  * Hawking  _et al._ (1982) S. W. Hawking, I. G. Moss, and J. M. Stewart, “Bubble Collisions in the Very Early Universe,” Phys. Rev. D 26, 2681 (1982).
  * Baker  _et al._ (2021) Michael J. Baker, Moritz Breitbach, Joachim Kopp, and Lukas Mittnacht, “Primordial Black Holes from First-Order Cosmological Phase Transitions,” (2021), arXiv:2105.07481 [astro-ph.CO] .
  * Hawking (1989) S. W. Hawking, “Black Holes From Cosmic Strings,” Phys. Lett. B 231, 237–239 (1989).
  * Polnarev and Zembowicz (1991) Alexander Polnarev and Robert Zembowicz, “Formation of Primordial Black Holes by Cosmic Strings,” Phys. Rev. D 43, 1106–1109 (1991).
  * MacGibbon  _et al._ (1998) Jane H. MacGibbon, Robert H. Brandenberger, and U. F. Wichoski, “Limits on black hole formation from cosmic string loops,” Phys. Rev. D 57, 2158–2165 (1998), arXiv:astro-ph/9707146 .
  * Rubin  _et al._ (2001) Sergei G. Rubin, Alexander S. Sakharov, and Maxim Yu. Khlopov, “The Formation of primary galactic nuclei during phase transitions in the early universe,” J. Exp. Theor. Phys. 91, 921–929 (2001), arXiv:hep-ph/0106187 .
  * Deng  _et al._ (2017) Heling Deng, Jaume Garriga, and Alexander Vilenkin, “Primordial black hole and wormhole formation by domain walls,” JCAP 04, 050 (2017), arXiv:1612.03753 [gr-qc] .
  * Ferrer  _et al._ (2019) Francesc Ferrer, Eduard Masso, Giuliano Panico, Oriol Pujolas, and Fabrizio Rompineve, “Primordial Black Holes from the QCD axion,” Phys. Rev. Lett. 122, 101301 (2019), arXiv:1807.01707 [hep-ph] .
  * Deng and Vilenkin (2017) Heling Deng and Alexander Vilenkin, “Primordial black hole formation by vacuum bubbles,” JCAP 12, 044 (2017), arXiv:1710.02865 [gr-qc] .
  * Deng (2020) Heling Deng, “Primordial black hole formation by vacuum bubbles. Part II,” JCAP 09, 023 (2020), arXiv:2006.11907 [astro-ph.CO] .
  * Maeso  _et al._ (2022) David Navidad Maeso, Luca Marzola, Martti Raidal, Ville Vaskonen, and Hardi Veermäe, “Primordial black holes from spectator field bubbles,” JCAP 02, 017 (2022), arXiv:2112.01505 [astro-ph.CO] .
  * Kitajima and Takahashi (2020) Naoya Kitajima and Fuminobu Takahashi, “Primordial Black Holes from QCD Axion Bubbles,” JCAP 11, 060 (2020), arXiv:2006.13137 [hep-ph] .
  * Li (2023) Hai-Jun Li, “QCD axion bubbles from the hidden SU(N) gauge symmetry breaking,” (2023), arXiv:2303.04537 [hep-ph] .
  * Dolgov  _et al._ (2009) A. D. Dolgov, M. Kawasaki, and N. Kevlishvili, “Inhomogeneous baryogenesis, cosmic antimatter, and dark matter,” Nucl. Phys. B 807, 229–250 (2009), arXiv:0806.2986 [hep-ph] .
  * Hasegawa and Kawasaki (2019) Fuminori Hasegawa and Masahiro Kawasaki, “Primordial Black Holes from Affleck-Dine Mechanism,” JCAP 01, 027 (2019), arXiv:1807.00463 [astro-ph.CO] .
  * Kawasaki and Murai (2019) Masahiro Kawasaki and Kai Murai, “Formation of supermassive primordial black holes by Affleck-Dine mechanism,” Phys. Rev. D 100, 103521 (2019), arXiv:1907.02273 [astro-ph.CO] .
  * Peccei and Quinn (1977a) R.D. Peccei and Helen R. Quinn, “CP Conservation in the Presence of Instantons,” Phys. Rev. Lett. 38, 1440–1443 (1977a).
  * Peccei and Quinn (1977b) R.D. Peccei and Helen R. Quinn, “Constraints Imposed by CP Conservation in the Presence of Instantons,” Phys. Rev. D 16, 1791–1797 (1977b).
  * Weinberg (1978) Steven Weinberg, “A New Light Boson?” Phys. Rev. Lett. 40, 223–226 (1978).
  * Wilczek (1978) Frank Wilczek, “Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons,” Phys. Rev. Lett. 40, 279–282 (1978).
  * Preskill  _et al._ (1983) John Preskill, Mark B. Wise, and Frank Wilczek, “Cosmology of the Invisible Axion,” Phys. Lett. B 120, 127–132 (1983).
  * Abbott and Sikivie (1983) L.F. Abbott and P. Sikivie, “A Cosmological Bound on the Invisible Axion,” Phys. Lett. B 120, 133–136 (1983).
  * Dine and Fischler (1983) Michael Dine and Willy Fischler, “The Not So Harmless Axion,” Phys. Lett. B 120, 137–141 (1983).
  * Borsanyi  _et al._ (2016) Sz. Borsanyi _et al._ , “Calculation of the axion mass based on high-temperature lattice quantum chromodynamics,” Nature 539, 69–71 (2016), arXiv:1606.07494 [hep-lat] .
  * Grilli di Cortona  _et al._ (2016) Giovanni Grilli di Cortona, Edward Hardy, Javier Pardo Vega, and Giovanni Villadoro, “The QCD axion, precisely,” JHEP 01, 034 (2016), arXiv:1511.02867 [hep-ph] .
  * Turner (1986) Michael S. Turner, “Cosmic and Local Mass Density of Invisible Axions,” Phys. Rev. D 33, 889–896 (1986).
  * Lyth (1992) D. H. Lyth, “Axions and inflation: Sitting in the vacuum,” Phys. Rev. D 45, 3394–3404 (1992).
  * Visinelli and Gondolo (2009) Luca Visinelli and Paolo Gondolo, “Dark Matter Axions Revisited,” Phys. Rev. D 80, 035024 (2009), arXiv:0903.4377 [astro-ph.CO] .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Nomura  _et al._ (2016) Yasunori Nomura, Surjeet Rajendran, and Fabio Sanches, “Axion Isocurvature and Magnetic Monopoles,” Phys. Rev. Lett. 116, 141803 (2016), arXiv:1511.06347 [hep-ph] .
  * Kawasaki  _et al._ (2016) Masahiro Kawasaki, Fuminobu Takahashi, and Masaki Yamada, “Suppressing the QCD Axion Abundance by Hidden Monopoles,” Phys. Lett. B 753, 677–681 (2016), arXiv:1511.05030 [hep-ph] .
  * Kawasaki  _et al._ (2018) Masahiro Kawasaki, Fuminobu Takahashi, and Masaki Yamada, “Adiabatic suppression of the axion abundance and isocurvature due to coupling to hidden monopoles,” JHEP 01, 053 (2018), arXiv:1708.06047 [hep-ph] .
  * Chiba  _et al._ (2004) Takeshi Chiba, Fuminobu Takahashi, and Masahide Yamaguchi, “Baryogenesis in a flat direction with neither baryon nor lepton charge,” Phys. Rev. Lett. 92, 011301 (2004), [Erratum: Phys.Rev.Lett. 114, 209901 (2015)], arXiv:hep-ph/0304102 .
  * Takahashi and Yamaguchi (2004) Fuminobu Takahashi and Masahide Yamaguchi, “Spontaneous baryogenesis in flat directions,” Phys. Rev. D 69, 083506 (2004), arXiv:hep-ph/0308173 .
  * Higaki  _et al._ (2014) Tetsutaro Higaki, Kwang Sik Jeong, and Fuminobu Takahashi, “Solving the Tension between High-Scale Inflation and Axion Isocurvature Perturbations,” Phys. Lett. B 734, 21–26 (2014), arXiv:1403.4186 [hep-ph] .
  * Choi  _et al._ (1997) Kiwoon Choi, Hang Bae Kim, and Jihn E. Kim, “Axion cosmology with a stronger QCD in the early universe,” Nucl. Phys. B 490, 349–364 (1997), arXiv:hep-ph/9606372 .
  * Banks and Dine (1997) Tom Banks and Michael Dine, “The Cosmology of string theoretic axions,” Nucl. Phys. B 505, 445–460 (1997), arXiv:hep-th/9608197 .
  * Jeong and Takahashi (2013) Kwang Sik Jeong and Fuminobu Takahashi, “Suppressing Isocurvature Perturbations of QCD Axion Dark Matter,” Phys. Lett. B 727, 448–451 (2013), arXiv:1304.8131 [hep-ph] .
  * Takahashi and Yamada (2015) Fuminobu Takahashi and Masaki Yamada, “Strongly broken Peccei-Quinn symmetry in the early Universe,” JCAP 10, 010 (2015), arXiv:1507.06387 [hep-ph] .
  * Carr  _et al._ (2010) B. J. Carr, Kazunori Kohri, Yuuiti Sendouda, and Jun’ichi Yokoyama, “New cosmological constraints on primordial black holes,” Phys. Rev. D 81, 104019 (2010), arXiv:0912.5297 [astro-ph.CO] .
  * Rosas-Guevara  _et al._ (2016) Yetli Rosas-Guevara, Richard G. Bower, Joop Schaye, Stuart McAlpine, Claudio Dalla Vecchia, Carlos S. Frenk, Matthieu Schaller, and Tom Theuns, “Supermassive black holes in the EAGLE Universe. Revealing the observables of their growth,” Mon.Not.Roy.Astron.Soc. 462, 190–205 (2016), arXiv:1604.00020 [astro-ph.GA] .