# Long-range ballistic transport of Brown-Zak fermions in graphene
superlattices

Julien Barrier *  — 1,2 These authors contributed equally. Piranavan
Kumaravadivel *  — 1,211footnotemark: 1 Roshan Krishna Kumar *  — 1,2 L. A.
Ponomarenko *  — 1,3 Na Xin 1,2 Matthew Holwill *  — 2 Ciaran Mullan *  — 1
Minsoo Kim *  — 1 R. V. Gorbachev *  — 1,2 M. D. Thompson *  — 2,3 J. R.
Prance *  — 2,3 T. Taniguchi4 K. Watanabe *  — 4 I. V. Grigorieva *  — 1,2 K.
S. Novoselov *  — 1,2 A. Mishchenko *  — 1,2 V. I. Fal’ko *  — 1,2 A.K. Geim *
— 1,2 Correspondance for this work shold be addressed to A.I. Berdyugin
<EMAIL_ADDRESS>and A.K. Geim<EMAIL_ADDRESS>A.I.
Berdyugin *  — 1,2 22footnotemark: 2

(1 Department of Physics and Astronomy, University of Manchester, Manchester
M13 9PL, UK
2 National Graphene Institute, University of Manchester, Manchester M13 9PL,
UK
3 Department of Physics, University of Lancaster, Lancaster LA1 4YW, UK
4 National Institute for Materials Science, Ibaraki 305-0044, Japan
)

###### Abstract

In quantizing magnetic fields, graphene superlattices exhibit a complex
fractal spectrum often referred to as the Hofstadter butterfly. It can be
viewed as a collection of Landau levels that arise from quantization of Brown-
Zak minibands recurring at rational ($p/q$) fractions of the magnetic flux
quantum per superlattice unit cell. Here we show that, in graphene-on-boron-
nitride superlattices, Brown-Zak fermions can exhibit mobilities above 106
cm2V-1s-1 and the mean free path exceeding several micrometers. The
exceptional quality of our devices allows us to show that Brown-Zak minibands
are $4q$ times degenerate and all the degeneracies (spin, valley and mini-
valley) can be lifted by exchange interactions below 1K. We also found
negative bend resistance at $1/q$ fractions for electrical probes placed as
far as several micrometers apart. The latter observation highlights the fact
that Brown-Zak fermions are Bloch quasiparticles propagating in high fields
along straight trajectories, just like electrons in zero field.

## Introduction

Van der Waals assembly offers a possibility to create materials by stacking
atomically-thin layers of different crystals [1, 2, 3]. One of the simplest
and most studied van der Waals heterostructures is graphene encapsulated
between two hexagonal boron nitride (hBN) crystals. The encapsulation protects
graphene from extrinsic disorder [4, 5], allowing ultra-high electronic
quality and micrometer-scale ballistic transport often limited only by edge
scattering [6, 7]. A special case of encapsulated graphene heterostructures is
graphene superlattices where crystallographic axes of graphene and hBN are
intentionally aligned. A small (1.8%) mismatch between graphene and hBN
crystal lattices results in a periodic moiré potential acting on charge
carriers in graphene and leading to the formation of electronic minibands[1,
2, 3, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17].

A relatively large ($\approx$14 nm) periodicity of graphene-on-hBN
superlattices has also made it possible to study the regime of Hofstadter
butterflies, which requires the magnetic flux $\phi$ per superlattice unit
cell to be comparable to the flux quantum $\phi_{0}$ in experimentally-
accessible magnetic fields $B$. At $B=B_{p/q}$ corresponding to
$\phi=\phi_{0}p/q$, where $p$ and $q$ are integer, the translational symmetry
of the electronic system is restored (despite the presence of a quantizing
magnetic field) and the superlattice’s electronic spectrum can again be
described in terms of Bloch states [18, 19, 20, 21, 22, 23, 24, 25, 26, 27,
28, 29], just like in $B=0$. These high-field Bloch states are characterized
by their own miniband spectra [25, 26, 27, 28, 29] different from the zero-$B$
spectrum of graphene-on-hBN superlattices. The associated quasiparticles are
referred to as Brown-Zak (BZ) fermions. According to the group-theory
analysis, an electronic spectrum for each realization of BZ fermions should
have an additional $q$-fold degeneracy [18, 19, 20, 21, 22, 23, 24, 25, 26,
27, 30, 31] (that is, contains $q$ equivalent mini-valleys). This degeneracy
is additional to the 4-fold spin and valley degeneracy of graphene’s original
spectrum. Importantly, BZ fermions are Bloch quasiparticles, like electrons in
solids or Dirac fermions, and, at $B=B_{p/q}$, they move through the
superlattice as if the applied field is zero. Away from these exact values, BZ
fermions experience an effective magnetic field $B_{eff}$ = $B-B_{p/q}$ (refs
[18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 30, 31]).

## Results

### Experimental devices and measurement setup

We report the electronic properties of BZ fermions with different $p$ and $q$
using high-quality graphene superlattices. These devices were fabricated using
the standard dry transfer procedures, where the studied graphene crystal was
carefully aligned with one of the encapsulating hBN crystals using the
crystallographic edges[11]. The alignment was verified by Raman
spectroscopy[32] prior to encapsulation with the second hBN crystal. The
latter was intentionally misaligned to avoid competing moiré patterns[33, 34,
35]. The assembled stacks were placed on an oxidized Si wafer, which allowed
us to apply the back-gate voltage $V_{g}$ to control the carrier density $n$.
We studied six devices that were shaped into the multiterminal Hall bar
geometry and had the main channel widths $W$ ranging from 2 to 17 µm (see Fig.
1a and Supplementary Note 1). The devices were first characterized by
measuring their longitudinal resistivity $\rho_{xx}$ in zero $B$ and Hall
resistivity $\rho_{xy}$ in small non-quantising $B$ below 0.1 T. The latter
enabled us to find the $n(V_{g})$ dependences, except for gate voltages close
to the neutrality points (NPs) and van Hove singularities (vHS), where
$\rho_{xy}$ reversed its sign and could no longer be described by the standard
dependence $\rho_{xy}=B/ne$ ($e$ is the electron charge). All our devices
exhibited very high carrier mobilities $\mu=(\rho_{xx}ne)^{-1}$ of the order
of 106 cm2 V-1 s-1, which were still somewhat reduced by edge scattering
because of finite W (Fig. 1b). We corroborated the high quality of our devices
using transverse magnetic focusing measurements (Supplementary Note 2). In the
main text, we focus on two large-width devices (D1 and D2) exhibiting highest
$\mu$.

Figure 1: High-quality graphene superlattices and their transport properties.
a, Optical micrograph of one of our devices (D1; twist angle $\theta$ between
graphene and hBN of about 0.4∘). The Hall bar is seen in violet with golden
electrical contacts. b, Mobility and mean free path for D1 measured at zero
$B$ and 10 mK. Semitransparent vertical strips indicate the doping regions
around NPs and vHS where $n$ could not be extracted directly from Hall
measurements and charge inhomogeneity also plays a role. To calculate $\mu$
and $l$ within the shaded regions, we assumed a constant gate capacitance and
linearly extrapolated the $n(V_{g})$ dependences found sufficiently far from
NPs and vHS (Supplementary Note 3). The noisy behavior at large values of
$\mu$ and $l$ arises from $\rho_{xx}$ becoming small ($\approx$1 , about 4
orders of magnitude smaller than that at the NPs). The horizontal black line
indicates the device width $W$. c, $\sigma_{xx}(V_{g},B)$ measured by sweeping
$V_{g}$ and varying $B$ in small steps of 40 mT. $T$ = 10 and 250 mK below and
above 14 T, respectively. Indigo-to-yellow colors: Log scale truncated between
38 nS and 16 mS for $B<$ 14 T and between 4 nS and 0.4 mS above 14 T. White
rectangles: these regions are shown in finer detail in Figs. 3 and 4. d, Same
as in panel b but for $\phi/\phi_{0}=1/2$ ($B\approx$ 15 T); $T$ = 250 mK. In
addition to NPs and vHS, the grey strips also cover a wide region of the
quantum Hall regime ($|V_{g}|<$ 20 V), which is dominated by large cyclotron
gaps in the main graphene spectrum. The transport data used to calculate $\mu$
and $l$ in panel (b) and (d) are shown in Supplementary Note 3.

Fig. 1c shows a map of the longitudinal conductivity
$\sigma_{xx}=\rho_{xx}/(\rho_{xx}^{2}+\rho_{xy}^{2})$ as a function of $V_{g}$
and $B$ in fields up to 18 T. The main features on such maps have been well
understood in terms of the Hofstadter spectrum for Dirac fermions in moiré
superlattices [11, 12, 13, 14, 15, 16, 17, 22]. One can see numerous Landau
levels (LLs) fanning out from the main and secondary NPs (miniband edges)
which are located at $V_{g}$ near 0 and $\pm$45 V, respectively. There are
also pronounced funneling features near van Hove singularities at
$V_{g}\approx$ -60, -35, +40 and +55 V in both hole ($-$) and electron ($+$)
parts of the spectrum. Another important attribute of the map is horizontal
yellow streaks that occur at $\phi=\phi_{0}p/q$ (Fig. 1c). If temperature $T$
is increased above 100 K so that Landau quantization is strongly suppressed,
the horizontal streaks become the only dominant feature on such transport maps
[28, 29]. In the high-T regime, the streaks represent oscillations in both
$\rho_{xx}(B)$ and $\rho_{xy}(B)$ at a constant $V_{g}$. Their $1/B$ frequency
is independent of $n$, and they were named BZ oscillations28. The horizontal
streaks are seen in Fig. 1c correspond to maxima in $\sigma_{xx}$ and zeros in
$\rho_{xy}$ and, as explained in the introduction, reflect the recovery of
translational symmetry ($p$ flux quanta penetrate through $q$ superlattice
unit cells) and the emergence of Bloch states experiencing zero $B_{eff}$.
Along each horizontal streak, one can find numerous NPs and vHS, which reflect
different realizations of BZ fermions at each $B_{p/q}$. Landau mini-fans
radiate from these NPs in both directions along the $B$ axis (Fig. 1c),
representing Landau quantization of BZ electronic spectra by non-zero
$B_{eff}$ (see below) [18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29].
Another notable feature of the shown $\sigma_{xx}$ map is a repetitive
triangular-like pattern seen most clearly between 2 and 12 T, especially for
positive $V_{g}$. The yellow triangles are made of horizontal streaks at zero
$B_{eff}$, vertical streaks emerging from NPs for BZ fermions and slanted
streaks originating from vHS.

### Ballistic transport

The transport properties of BZ fermions can be analyzed in the same way as in
Fig. 1b for Dirac fermions. The results are plotted in Fig. 1d for the case of
$\phi/\phi_{0}=1/2$ and show that, away from NPs and vHS, BZ fermions in our
devices exhibit $\mu$ reaching a few 106 cm2 V-1 s-1. This is comparable to
$\mu$ of Dirac fermions in zero $B$. Another important characteristic of
charge carrier transport is the mean free path $l$. It can be evaluated for
both Dirac and BZ fermions, using the standard expression
$\sigma_{xx}=ge^{2}/h((k_{F}l)/2)$ where $h$ is the Planck constant, $k_{F}$
is the Fermi momentum and $g$ is the degeneracy. In zero $B$, $g=4$ because of
the spin and valley degeneracy of Dirac fermions. BZ fermions are expected to
have an additional, mini-valley degeneracy [18, 19, 20, 21, 22, 23, 24, 25,
26, 27, 30, 31], which is equal to $q$ (that is, $g=8$ for the case of Fig.
1d). Using these $g$, we calculated $l$ as shown in the bottom panels of Figs.
1b and 1d. The mean free path in zero $B$ reached $>$ 20 µm, implying that
ballistic transport in our devices is limited by $W$ rather than impurity
scattering. For some realizations of BZ fermions, their mean free path also
exceeds 10 µm, marginally smaller than $l$ for Dirac fermions (cf. Figs 1b and
1d). Supplementary Note 4 provides similar analysis for $\mu$ and $l$ at other
values of $q$.

We corroborate the existence of ballistic BZ fermions using so-called bend
resistance geometry [5, 36] sketched in Fig. 2a. The geometry allows one to
detect if charge carriers can move ballistically over the entire channel width
$W$, along straight trajectories connecting current and voltage contacts (see
caption of Fig. 2a). If this is the case, ballistic transport gives rise to
the negative sign of the bend resistance $R_{b}$(ref [36]) in contrast to its
conventional, positive sign for diffusive (ohmic) transport. As expected from
very long $l$ of Dirac fermions, $R_{b}$ was negative in zero $B$ everywhere
away from NPs and vHS (Fig. 2b), confirming further the high quality of our
superlattice devices. Finite $B$ bend Dirac fermion trajectories and, as
expected, $R_{b}$ rapidly reversed its sign [36] with increasing $B$ (inset of
Fig. 2b). Remarkably, our devices exhibited negative $R_{b}$ also in high
$B=B_{p/q}$, thus revealing straight trajectories over distances of several µm
(Figs. 2c,d and Supplementary Note 5). For comparison, the corresponding map
of $\rho_{xx}$ is provided in Supplementary Note 6, which shows that the
measured longitudinal resistance always remained positive. The profound
negative pockets in $R_{b}(V_{g},B)$ appeared only around $B$ corresponding to
$\phi/\phi_{0}=1/2$, $1/3$, $1/4$ and $1/5$. This is the regime where the
existence of BZ fermions, experiencing zero $B_{eff}$, was previously inferred
[28, 29] from maxima in $\sigma_{xx}$ and zeros in $\sigma_{xy}$. The
possibility of ballistic transport of BZ fermions was also suggested using
numerical simulations [15]. The observed negative $R_{b}$ prove the previous
conjectures unequivocally. The negative pockets in Fig. 2c are located between
NPs and vHS for BZ fermions, similar to the case of Dirac fermions. Away from
$B=B_{p/q}$, non-zero $B_{eff}$ bends the BZ fermion trajectories and the
negative signal disappeared (inset of Fig. 2d), again like in the case of
other Bloch quasiparticles (electrons and Dirac fermions). Our devices
exhibited long-range ballistic transport of BZ fermions only for unit
fractions, $\phi_{0}/q$. For example, Fig. 2c reveals a pronounced set of
negative $R_{b}$ pockets at $\phi/\phi_{0}=1/5$ but no negative signal was
observed for $2/5$ and $3/5$. This behavior probably stems from lower $\mu$ of
BZ fermions with $p>1$, which can be attributed to their larger effective
masses [29].

Figure 2: Ballistic transport of BZ fermions over micrometer distances. a,
Schematic of bend resistance measurements. Current $I$ is applied between
contacts 3 and 4, and voltage $V_{b}$ is measured between 2 and 1, yielding
the bend resistance, $R_{b}$ = $V_{b}/I$. The voltage is positive for
diffusive transport but becomes negative, if charge carriers move directly
from current injecting contact 3 into voltage probe 1 (as shown by the red
arrow). b, Bend resistance for Dirac fermions in zero $B$ (device D2 with $W$
= 4 µm and $\theta\approx$ 0∘). Inset: $R_{b}(B)$ taken at the minimum
indicated by the arrow in the main plot. c, Map $R_{b}(V_{g},B)$ for the same
device. $B$ was changed in steps of 50 mT. Pockets of negative $R_{b}$ appear
along $\phi/\phi_{0}=1/q$ and are seen in magenta. d, Cross-sections from
panel c for $q=$ 2 and 3. The inset shows sign reversals in $R_{b}$ plotted as
a function of $B_{eff}=B–B_{p/q}$ for the minima marked by the color-coded
arrows. $T$ = 2 K for all the plots.

### Lifting degeneracy in BZ minibands

The high mobility of BZ fermions resulted in their Landau quantization in
small $B_{eff}<$ 1 T, which allowed us to find experimentally the spectral
degeneracy $g$ for different $p/q$. To this end, Fig. 3a shows a high-
resolution map of $\sigma_{xx}$ between $\phi/\phi_{0}=1/2$ and $1/4$ (part of
Fig. 1c). One can clearly see many Landau mini-fans originating from NPs for
different realizations of BZ fermions. For example, there are three profound
mini-fans spreading from $B\approx$ 9.7 T ($\phi/\phi_{0}=1/3$) for both
negative and positive $B_{eff}$. For clarity, the conductivity map in Fig. 3a
is replotted in Fig. 3b by tracing all well-defined minima in $\sigma_{xx}$.
Each minimum can be described by its integer filling factor $\nu$, which we
calculated from the minimum’s slope. This representation of the Hofstadter
spectrum, where LLs are plotted as a function of $n$ or $V_{g}$ rather than
energy, is usually referred to as the Wannier diagram [11, 12, 13, 14, 15, 16,
17, 21, 22]. In the diagram of Fig. 3b, one can identify LLs for BZ-fermion
realizations at $q$ = 2, 3, 4, 5, 7, 8, 9 and 11, and for $p$ = 1, 2, 3, 4 and
5. The difference in $\nu$ between the nearest LLs yields directly their $g$.
For example, all LLs found for $\phi/\phi_{0}=1/2$ were separated by
$\Delta\nu$ = 2 whereas those at $\phi/\phi_{0}=1/3$ and $1/4$ by $\Delta\nu$
= 3 and 4, respectively (see Fig. 3b). Therefore, the observed degeneracies
were equal to $q$. Note that the measured Hall conductance also exhibited
quantized values in steps of $qe^{2}/h$ (Supplementary Note 7). Examining the
Wannier diagram further, we find that $\Delta\nu$ was equal to $q$,
independently of numerator $p$ and for all LLs stemming from the same
$B_{p/q}$. This observation agrees with the theoretical expectation that there
should be $q$ equivalent mini-valleys for each realization of BZ fermions [18,
19, 20, 21, 22, 23, 24, 25, 26, 27, 30, 31]. The agreement takes into account
that both spin and valley degeneracies for BZ fermions were lifted by exchange
interactions [14], as in the case of the main Dirac spectrum with its clearly
lifted degeneracies (see the LLs marked in black in Fig. 3b). Therefore, the
total degeneracy for BZ minibands in Fig. 3 was $g=4q$. This is further
supported by our measurements at a relatively high $T$ of 2 K that suppressed
LLs with lifted spin and valley degeneracies (Supplementary Note 6). The same
behavior (Landau mini-fans exhibiting the sequence $\Delta\nu=q$) was also
observed in other parts of the Wannier diagram for both electron and hole
doping (see, e.g., Fig. 4).

Figure 3: Landau quantization in BZ minibands at 10 mK. a, High-resolution map
$\sigma_{xx}(V_{g},B)$ for the electron-doped region indicated in Fig. 1c by
the white rectangle (device D1). To better resolve LLs around $q=3$, $B$ was
changed in steps of 10 to 20 mT whereas data in other regions were acquired
using 40 to 80 mT steps. This has resulted in the contrast discontinuities
seen in the map. Log color scale: indigo (230 nS) to yellow (7.8 mS) for the
entire map. b, Minima from panel (a) are shown schematically. The color-coded
numbers are the filling factors for the corresponding LLs. Thick black lines
correspond to the main sequence of LLs for graphene’s Dirac spectrum (spin and
valley degeneracy lifted). The green, red, navy, blue, orange, magenta, pink
and yellow lines correspond to $q=$ 2, 3, 4, 5, 7, 8, 9 and 11, respectively.
Dashed red lines: minima due to lifted mini-valley degeneracy.

There is however one notable exception, which was observed for
$\phi/\phi_{0}=1/3$ (left part of the Wannier diagram in Fig. 3b; dashed
lines). In this case, LLs of BZ fermions with $q=3$ are separated only by
$\Delta\nu=1$ so that all consecutive LLs from 3 to 12 are present on the fan
diagram. The indicated minima in $\sigma_{xx}$ were rather fragile, being
rapidly smeared by $T$ or excitation current (Supplementary Note 7). Landau
levels with $\Delta\nu$ = 1 at $\phi/\phi_{0}=1/q$ do not exist within the
single-particle Hofstadter-Wannier model [16, 17, 21, 22]. Such _extra_ LLs
have been observed previously and referred to as the fractional Bloch quantum
Hall effect (FBQHE) [16] or symmetry-broken Chern insulator (SBCI) [17].
Theoretically, the many-body states can arise from interaction effects that
lift the mini-valley degeneracy of BZ fermions. This may happen via mini-
valley mixing due to the formation of charge-density waves (commensurate with
the magnetic superlattice composed of $q$ unit cells of the underlying moiré
pattern) or a Wigner crystal [37] (states localized in a part of the magnetic
supercell). Alternatively, the fully lifted degeneracy may occur through
spontaneous mini-valley polarization of BZ LLs, analogous to spin/valley
ferromagnetism in graphene [14]. In the latter case, the BZ fermion states
localize in the momentum space around one of the mini-valleys but remain
delocalized across the magnetic supercell.

## Discussion

Figure 4: Anomalous behavior of Landau levels in BZ minibands. a, High-
resolution $\sigma_{xx}(V_{g},B)$ for the hole-doped region marked in Fig. 1c
(device D1). The measurements below $B=$ 12.5 T were carried out at 10 mK, and
at 250 mK in the fields above. $B$ was applied in steps of 16 and 40 mT above
and below 12.5 T, respectively, except for the low left corner where the data
were acquired at higher resolution of 10 mT (seen as discontinuities in the
contrast). Log scale: indigo (310 nS) to yellow (0.3 mS). b, Schematics for
the conductance minima found in (a). The same color coding as in Fig. 3b. The
solid lines indicate LLs evolving as expected, linearly in $B$ and $V_{g}$.
The thin curves indicate the anomalous bending whereas the dotted curves show
the staircase-like evolution observed for some LLs. All the anomalous features
were highly reproducible and, for example, did not depend on the step size in
$B$.

Some Landau mini-fans exhibited highly anomalous behavior at low $T$, which
cannot be understood either within the Hofstadter-Wannier model or by
considering LLs of non-interacting BZ fermions and, to the best of our
knowledge, was never reported before. From either perspective, individual LLs
must evolve linearly in $V_{g}$ and $B$, as indeed seen in most cases (e.g.,
Fig. 3). This is because the density of states on each BZ-fermion LL is
proportional to $B_{eff}=B-B_{p/q}$. Also, the long ($\approx$300 nm) distance
from graphene to the gate in our devices makes quantum capacitance corrections
[38, 39] negligible, leading to the $n\propto V_{g}$ dependence. Unexpectedly,
we found that some Landau fans of BZ fermions exhibited bending and staircase-
like features. Figure 4 shows examples of such behavior for the case of hole
doping. Both bending and staircases appear away from NPs, in the regions close
to BZ fermions’ vHS. The observed bending towards the gate-voltage axis
indicates that, in addition to the visible LLs, there can be some other
electronic states that are populated in parallel. Because electron-electron
interactions clearly play a significant role under the reported experimental
conditions (e.g., they lift all the spectral degeneracies), it is reasonable
to expect that interactions are also involved in the described anomalies.
However, the usual suspects, such as negative compressibility cannot possibly
explain our findings (in the latter case, LLs would bend toward the $B$ axis).
The anomalous behavior is possibly due to an interplay of BZ-fermion LLs with
quantized states originating from nearby vHS, which would lead to
redistribution of charge carriers between states with the light and heavy
effective masses. One possible scenario is localization of electrons within
some parts of the magnetic supercell, because its size becomes notably larger
than the magnetic length at $B>$ 10 T (Wigner crystallization mentioned
above). Unfortunately, we could not find any additional features that would
enable us to decipher origins of the described anomalies and, therefore, have
to leave them for further investigation.

#### Methods

The reported graphene-on-hBN superlattices were assembled using the standard
dry transfer procedures and polydimethylsiloxane (PDMS) stamps coated with a
polypropylcarbonate (PPC) layer7. The thicknesses of the used hBN crystals was
between 20 and 70 nm. After the trilayer stack was assembled, we used the
standard electron beam lithography and reactive ion etching to define trenches
for electrical contacts. Cr/Au films were evaporated into the trenches.
Finally, using ion etching again, the trilayer stack was shaped into a Hall
bar mesa. For electrical measurements, the standard low-frequency lock-in
technique was employed. At temperatures below 100 mK, we used excitation
currents of $\approx$10 nA, and between 0.1 to 1 µA at higher $T$. The carrier
concentration was controlled by applying a DC gate voltage between graphene
and the silicon substrate.

## References

  * [1] Geim, A. K. & Grigorieva, I. V.  Van der waals heterostructures.  _Nature_ 499, 419–425 (2013).
  * [2] Novoselov, K., Mishchenko, o. A., Carvalho, o. A. & Neto, A. C.  2d materials and van der waals heterostructures.  _Science_ 353 (2016).
  * [3] Yankowitz, M., Ma, Q., Jarillo-Herrero, P. & LeRoy, B. J.  van der waals heterostructures combining graphene and hexagonal boron nitride.  _Nature Reviews Physics_ 1, 112–125 (2019).
  * [4] Dean, C. R. _et al._ Boron nitride substrates for high-quality graphene electronics.  _Nature nanotechnology_ 5, 722–726 (2010).
  * [5] Mayorov, A. S. _et al._ Micrometer-scale ballistic transport in encapsulated graphene at room temperature.  _Nano letters_ 11, 2396–2399 (2011).
  * [6] Wang, L. _et al._ One-dimensional electrical contact to a two-dimensional material.  _Science_ 342, 614–617 (2013).
  * [7] Kumaravadivel, P. _et al._ Strong magnetophonon oscillations in extra-large graphene.  _Nature communications_ 10, 1–6 (2019).
  * [8] Chen, S.-C., Kraft, R., Danneau, R., Richter, K. & Liu, M.-H.  Electrostatic superlattices on scaled graphene lattices.  _Communications Physics_ 3, 1–7 (2020).
  * [9] Wallbank, J., Patel, A., Mucha-Kruczyński, M., Geim, A. & Fal’Ko, V.  Generic miniband structure of graphene on a hexagonal substrate.  _Physical Review B_ 87, 245408 (2013).
  * [10] Yankowitz, M. _et al._ Emergence of superlattice dirac points in graphene on hexagonal boron nitride.  _Nature Physics_ 8, 382–386 (2012).
  * [11] Ponomarenko, L. _et al._ Cloning of dirac fermions in graphene superlattices.  _Nature_ 497, 594–597 (2013).
  * [12] Dean, C. R. _et al._ Hofstadter’s butterfly and the fractal quantum hall effect in moiré superlattices.  _Nature_ 497, 598–602 (2013).
  * [13] Hunt, B. _et al._ Massive dirac fermions and hofstadter butterfly in a van der waals heterostructure.  _Science_ 340, 1427–1430 (2013).
  * [14] Yu, G. _et al._ Hierarchy of hofstadter states and replica quantum hall ferromagnetism in graphene superlattices.  _Nature physics_ 10, 525–529 (2014).
  * [15] Diez, M., Dahlhaus, J., Wimmer, M. & Beenakker, C.  Emergence of massless dirac fermions in graphene’s hofstadter butterfly at switches of the quantum hall phase connectivity.  _Physical Review Letters_ 112, 196602 (2014).
  * [16] Wang, L. _et al._ Evidence for a fractional fractal quantum hall effect in graphene superlattices.  _Science_ 350, 1231–1234 (2015).
  * [17] Spanton, E. M. _et al._ Observation of fractional chern insulators in a van der waals heterostructure.  _Science_ 360, 62–66 (2018).
  * [18] Brown, E.  Bloch electrons in a uniform magnetic field.  _Physical Review_ 133, A1038 (1964).
  * [19] Zak, J.  Magnetic translation group.  _Physical Review_ 134, A1602 (1964).
  * [20] Azbel, M. Y.  Energy spectrum of a conduction electron in a magnetic field.  _Sov. Phys. JETP_ 19, 634–645 (1964).
  * [21] Wannier, G.  A result not dependent on rationality for bloch electrons in a magnetic field.  _physica status solidi (b)_ 88, 757–765 (1978).
  * [22] Hofstadter, D. R.  Energy levels and wave functions of bloch electrons in rational and irrational magnetic fields.  _Physical review B_ 14, 2239 (1976).
  * [23] Streda, P.  Quantised hall effect in a two-dimensional periodic potential.  _Journal of Physics C: Solid State Physics_ 15, L1299 (1982).
  * [24] Rhim, J.-W. & Park, K.  Self-similar occurrence of massless dirac particles in graphene under a magnetic field.  _Physical Review B_ 86, 235411 (2012).
  * [25] Rammal, R.  Landau level spectrum of bloch electrons in a honeycomb lattice.  _Journal de Physique_ 46, 1345–1354 (1985).
  * [26] Delplace, P. & Montambaux, G.  Semi-dirac point in the hofstadter spectrum.  _Phys. Rev. B_ 82, 035438 (2010).
  * [27] Chen, X. _et al._ Dirac edges of fractal magnetic minibands in graphene with hexagonal moiré superlattices.  _Physical Review B_ 89, 075401 (2014).
  * [28] Krishna Kumar, R. _et al._ High-temperature quantum oscillations caused by recurring bloch states in graphene superlattices.  _Science_ 357, 181–184 (2017).
  * [29] Kumar, R. K. _et al._ High-order fractal states in graphene superlattices.  _Proceedings of the National Academy of Sciences_ 115, 5135–5139 (2018).
  * [30] Thouless, D. J., Kohmoto, M., Nightingale, M. P. & den Nijs, M.  Quantized hall conductance in a two-dimensional periodic potential.  _Physical review letters_ 49, 405 (1982).
  * [31] Landau, L. D., Lifshitz, E. M. & Pitaevskii, L.  _Statistical physics Part 2 - The theory of the condensed state_ , vol. 9 (Pergamon Press, 1980).
  * [32] Woods, C. _et al._ Commensurate–incommensurate transition in graphene on hexagonal boron nitride.  _Nature physics_ 10, 451–456 (2014).
  * [33] Wang, L. _et al._ New generation of moiré superlattices in doubly aligned hbn/graphene/hbn heterostructures.  _Nano letters_ 19, 2371–2376 (2019).
  * [34] Wang, Z. _et al._ Composite super-moiré lattices in double-aligned graphene heterostructures.  _Science Advances_ 5, eaay8897 (2019).
  * [35] Finney, N. R. _et al._ Tunable crystal symmetry in graphene–boron nitride heterostructures with coexisting moiré superlattices.  _Nature Nanotechnology_ 14, 1029–1034 (2019).
  * [36] Beenakker, C. & van Houten, H.  Quantum transport in semiconductor nanostructures.  In _Solid state physics_ , vol. 44, 1–228 (Elsevier, 1991).
  * [37] DaSilva, A. M., Jung, J. & MacDonald, A. H.  Fractional hofstadter states in graphene on hexagonal boron nitride.  _Physical Review Letters_ 117, 036802 (2016).
  * [38] Chen, Z. & Appenzeller, J.  Mobility extraction and quantum capacitance impact in high performance graphene field-effect transistor devices.  In _2008 IEEE international electron devices meeting_ , 1–4 (IEEE, 2008).
  * [39] Yu, G. L. _et al._ Interaction phenomena in graphene seen through quantum capacitance.  _Proceedings of the National Academy of Sciences_ 110, 3282–3286 (2013).

#### Acknowledgements

This work was supported by the Engineering and Physical Sciences Research
Council (EPSRC), Lloyd’s Register Foundation, Graphene Flagship, the Royal
Society and the European Research Council (grant Vander). A.I.B. acknowledges
support from the Graphene NOWNANO Doctoral Training Centre.

#### Author contributions

P.K., N.X., and M.H. fabricated devices. J.B., R.K.K., L.A.P., C.M., M.K.,
A.I.B. performed transport measurements with help from A.M, M.D.T and J.R.P.
J.B., P.K., R.K.K., L.A.P., V.I.F, A.K.G. and A.I.B. performed data analysis.
T.T. and K.W. provided the hBN crystals. J.B., P.K., V.I.F., A.I.B. and A.K.G.
wrote the manuscript. All authors contributed to the discussions.

#### Competing interests

The authors declare no competing interests


Gravitational waves from type II axion-like curvaton model and its implication
for NANOGrav result

Masahiro Kawasaki(a,b) and Hiromasa Nakatsuka(a)

(a)ICRR, The University of Tokyo, Kashiwa, Chiba 277-8582, Japan
(b)Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583,
Japan

##  1 Introduction

North American Nanohertz Observatory for Gravitational Waves (NANOGrav) is one
of the Pulsar-timing-array experiments, which has been trying to detect the
gravitational wave (GW) signal through the long-term observation of pulsars.
Recently the NANOGrav reported a hint of the stochastic GW signal in their
12.5–year observation [1]. This signal can be explained by the stochastic GW
with $\Omega_{\rm GW}\sim 10^{-9}$ at $f\sim 10^{-8}\mathrm{Hz}$, which can be
sourced by cosmic string [2, 3, 4, 5, 6, 7], phase transitions [8, 9, 10, 11,
12], large density fluctuations in primordial black hole (PBH) formation
models [13, 14, 15, 16, 17, 18, 19, 20, 21, 22] and other origin [23, 24, 25].

The PBH formation scenario is an attractive candidate for the NANOGrav signal
since it can simultaneously explain the observed $30{M_{\odot}}$ black holes
in the binary merger events in LIGO-Virgo collaboration [26, 27, 28, 29]. The
PBH formation requires the large curvature power spectrum
($\mathcal{P}_{\zeta}(k_{\rm pbh})\sim 0.02$) at a wavenumber $k_{\rm pbh}$
related to the mass of the formed PBHs. Since PBHs are formed when the high-
density regions enter the horizon and collapse, the PBH mass is roughly given
by the horizon mass at the collapse, which is related to the frequency of the
density fluctuation $f$ ($=k_{\text{pbh}}/2\pi$) as [30]

$\displaystyle M(f)$ $\displaystyle\simeq
30{M_{\odot}}\left(\frac{\gamma}{0.2}\right)\left(\frac{\text{g}_{*}}{10.75}\right)^{-1/6}\left(\frac{f}{5.3\times
10^{-10}\,\mathrm{Hz}}\right)^{-2},$ (1)

where $\gamma\sim 0.2$ is a ratio of PBH mass to the horizon mass [31] and
$\text{g}_{*}$ is the numbers of relativistic degrees of freedom at the PBH
formation. The same density fluctuations that collapse to form PBHs also
generate a stochastic background of GWs through the nonlinear coupling,
$\Braket{\zeta\zeta h}$ when they enter the horizon [32, 33, 34], with a
spectrum peaking at the frequency $f\sim 10^{-9}\mathrm{Hz}$, very close to
the one of the NANOGrav signal.

In this paper, we focus on the axion-like curvaton model of the PBH formation
scenario [35, 36, 37, 15, 38]. There are two types of the axion-like curvaton
models; in one type the complex field $\Phi$ (whose phase part is the axion)
rolls down toward origin during inflation [35, 36, 37, 15], and in the other
type $\Phi$ rolls down from the origin [38]. These two types lead to the
different dynamics and the power spectrum of induced GWs. The former type
(type I) has been already studied and successfully explains the NANOGrav
signal [15]. The type II model was first proposed in [38] where the model
parameters were chosen to achieve a narrow density power spectrum of the
density perturbations. However, a broader power spectrum is preferable to
account for the NANOGrav signal. The goal of this paper is to investigate if
such a broad power spectrum can be produced also in the type II model. It is
found that by choosing appropriate parameters of the potential term in Eq.
(15), we obtain the desired broader spectral shape. Moreover, since our model
produces large positive non-Gaussianity on the density power spectrum, we can
achieve the required abundance of PBHs by the smaller amplitude of the density
power spectrum than the Gaussian one. We improve the treatment of the non-
Gaussianity of this model to accurately estimate the required amplitude of
density power spectrum, and also the induced GW.

This paper is organized as follows. We calculate the curvature power spectrum
of the type II axion-like curvaton model in Sec. 2. Next, we calculate the PBH
abundance in Sec.3 and induced GWs in Sec. 4. We can achieve the broad power
spectrum of induced GWs and explain the NANOGrav signal. In Sec. 5, we
conclude this paper.

##  2 Type II axion-like curvaton model

We now briefly summarize the dynamics of background and fluctuations in the
type II axion-like curvaton model. (See [38] for the detailed calculation.) In
this model, large fluctuations are produced in the phase direction of a
complex scalar field $\Phi$, which we call “axion-like curvaton”. The
potential of $\Phi$ is given by [38]

$\displaystyle
V_{\Phi}=\frac{\lambda}{4}\left(|\Phi|^{2}-\frac{v^{2}}{2}\right)^{2}+gI^{2}|\Phi|^{2}-v^{3}\epsilon(\Phi+\Phi^{*}),$
(2)

where $I$ is the inflaton field, $\lambda$ and $g$ are coupling constants,
$v/\sqrt{2}$ is the vacuum expectation value after inflation. The bias term,
$v^{3}\epsilon(\Phi+\Phi^{*})$, is introduced to avoid the cosmic string
problem and the stochastic effect on dynamics of $\Phi$, which requires the
complicated treatment of fluctuations. Note that $\epsilon$ is naturally small
($\epsilon\ll 1$) in the sense of ’t Hooft’s naturalness [39] since $U(1)$
symmetry is restored for $\epsilon=0$.

The field value of $\Phi$ changes during inflation due to the first and second
term of Eq. (2). In the early stage of inflation
($I\gtrsim(\lambda/g)^{1/2}v$), the interaction term with the inflaton fixes
$\Phi$ near the origin. In the late stage, the inflaton field value becomes
smaller and $\Phi$ starts to roll down the Higgs-like potential towards
$|\Phi|=v$. In this dynamics the phase direction acquires large fluctuations
$\sim H/|\Phi|$ ($H$: Hubble parameter during inflation).

The PBH formation requires the large fluctuations with wavenumber $k=k_{\rm
pbh}\simeq 10^{5}\,\mathrm{Mpc}^{-1}$, which corresponds to the scale of
$30{M_{\odot}}$ PBHs. Thus, we determine the inflaton coupling $g$ so that
$\Phi$ starts to roll down potential at $t_{\rm pbh}$ satisfying
$H(t_{\text{pbh}})=k_{\text{pbh}}/a(t_{\text{pbh}})$, which leads to

$\displaystyle g\simeq\frac{\lambda}{4}\left(\frac{v}{I(t_{\rm
pbh})}\right)^{2}.$ (3)

We also assume that the effective mass of $\Phi$ is much larger than the
Hubble parameter ($\lambda v^{2}\gg H^{2}$) until $t_{\rm pbh}$ to suppress
the fluctuation at the larger wavelengths. This condition also ensures
independence on the initial condition since $\Phi$ settles near the origin of
the potential independently of the initial field value of $\Phi$.

Based on the above background field dynamics, we calculate the fluctuation of
$\Phi$, which is decomposed as

$\displaystyle\Phi=\frac{1}{\sqrt{2}}(\varphi_{0}+\varphi)e^{i\theta}$ (4)

with the homogeneous solution $\varphi_{0}$, and the perturbations $\varphi$
and $\theta$. The phase direction $\theta$ works as the curvaton, and the
canonical field of the phase direction is defined as
$\tilde{\sigma}\equiv\varphi_{0}\theta$. The inflation induces the quantum
fluctuations with amplitude $H/(2\pi)$ for $\tilde{\sigma}$ at the horizon
crossing. Thus, the power spectrum of the $\theta$ fluctuations is given by

$\displaystyle\mathcal{P}_{\theta}(k,t_{k})=\left(\frac{H(t_{k})}{2\pi\varphi_{0}(t_{k})}\right)^{2},$
(5)

where $t_{k}$ is the time when the fluctuation with $k$ crosses the horizon.
$\mathcal{P}_{\theta}(k,t_{k})$ is suppressed for $k>k_{\text{pbh}}$ since the
$\varphi_{0}(t_{k})$ quickly grows after $t_{\text{pbh}}$. While the $\Phi$ is
fixed near the origin before $t_{\text{pbh}}$, $\tilde{\sigma}$ has the large
effective mass given by

$\displaystyle\tilde{m}_{\sigma}^{2}$
$\displaystyle\equiv\frac{\partial^{2}\left(-v^{3}\epsilon(\Phi+\Phi^{*})\right)}{\partial\tilde{\sigma}^{2}}\bigg{|}_{\tilde{\sigma}=0}={\sqrt{2}\epsilon
v^{2}\frac{v}{\varphi_{0}}},$ (6)

where we choose the model parameter so as $m_{\tilde{\sigma}}^{2}>H^{2}$ for
$t<t_{\text{pbh}}$. Such large effective mass damps the fluctuations as
$\tilde{\sigma}\propto a^{-3/2}$, and we define the damping factor of the
phase direction with momentum $k$ as

$\displaystyle R_{k}$
$\displaystyle\equiv\left(\frac{\tilde{\sigma}_{k}(t_{\text{pbh}})}{\tilde{\sigma}_{k}(t_{k})}\right)\sim\left(\frac{k}{k_{\text{pbh}}}\right)^{3/2}\text{
for }k<k_{\text{pbh}}.$ (7)

($R_{k}=1$ for $k>k_{\text{pbh}}$.) Finally, the power spectrum of the
$\theta$ fluctuations at the end of inflation $t_{e}$ is given by

$P_{\theta}(k,t_{e})=R_{k}^{2}\mathcal{P}_{\theta}(k,t_{k}).$ (8)

Let us evaluate the density perturbations induced by the fluctuations of the
phase direction. We suppose that the curvaton obtains the axion-like potential
after inflation through some nonperturbative effect. The potential minimum of
the nonperturbative term does not coincide with that of the primordial one
determined by the bias term in our model. Suppose that nonperturbative
potential takes the minimum at $\theta=\theta_{i}$, then it is written as

$\displaystyle
V_{\sigma}=\Lambda^{4}\left[1-\cos\left(\Theta\right)\right]\simeq\frac{1}{2}m_{\sigma}^{2}\sigma^{2},$
(9)

where $\Theta\equiv\theta-\theta_{i}$, the curvaton $\sigma$ is defined as
$\sigma\equiv v\Theta$ and $m_{\sigma}=\Lambda^{2}/v$. We assume that
$m_{\sigma}$ is small enough to neglect the axion-like potential during
inflation, i.e. $H^{2}\gg m_{\sigma}^{2}$. The density fluctuation is given by
${\delta\rho_{\sigma}}/{\rho_{\sigma}}=2\delta\theta/\theta_{i}$.

Neglecting a small contribution from the inflaton ($\mathcal{P}_{\zeta}\sim
10^{-10}$) compared to the curvaton, the power spectrum of curvature
perturbations is given by

$\displaystyle\mathcal{P}_{\zeta}(k)$
$\displaystyle=\left(\frac{r}{4+3r}\right)^{2}\left(\frac{2}{\theta_{i}}\right)^{2}\mathcal{P}_{\theta}(k,t_{k})R_{k}^{2},$
(10)

where $r$ is the ratio of the energy density of the curvaton to that of the
inflaton (or radiation after reheating), $r=\rho_{\sigma}/\rho_{I}$. Assuming
the instant reheating at $t_{\rm reheat}$, the ratio is given by $r(t_{\rm
reheat})=(v^{2}\theta_{i}^{2})/(6M_{pl}^{2})$, which is chosen to be small to
ensure that $\Phi$ does not disturb the inflation. $r$ grows during radiation-
dominated era as $r\propto a$ due to the matter-like behavior of the curvaton,
and its growth ends at the curvaton decay, $t_{\rm decay}$. We require that
the curvaton decays into radiation before it overcloses the universe,
$r(t_{\rm decay})\sim 0.5$, at which the temperature of the universe is
$T_{\rm decay}\sim(r(t_{\rm reheat})/r(t_{\rm decay}))T_{\rm reheat}$. The
typical decay rate of $\Phi$ is related to the curvaton mass as
$\Gamma_{\sigma}=\kappa^{2}m_{\sigma}^{3}/(16\pi v^{2})$ where $\kappa$ is a
coupling constant. We confirm that $T_{\rm decay}\sim 10^{3}\mathrm{GeV}$ is
achieved for $m_{\sigma}\sim 10^{8}\mathrm{GeV}$ for our parameters in Eq.
(15). In the following, $r$ refers to the energy ratio after decay, that is
$r\equiv r(t_{\rm decay})$.

The PBH abundance highly depends on the eventual non-Gaussian distribution of
$\mathcal{P}_{\zeta}(k)$. It is known that the axion-like curvaton models
produce large non-Gaussianity with local type bispectrum characterized by the
following parameter: [37]

$\displaystyle f_{\rm
NL}=\frac{5}{12}\left(-3+\frac{4}{r}+\frac{8}{4+3r}\right).$ (11)

We discuss the enhancement of the PBH abundance by non-Gaussianity in Sec.3.
Non-Gaussianity also affects the power spectrum of the curvature perturbations
and induced gravitational wave through the higher-order correlations [40, 41],
whose contribution is characterized by $\mathcal{P}_{\zeta}f_{\rm NL}^{2}$. We
take non-Gaussianity into account only approximately since the Gaussian
contribution is the dominant for our choice of parameters given by Eq. (15),
$\mathcal{P}_{\zeta}f_{\rm NL}^{2}<0.04$. We estimate non-Gaussian
amplification on the curvature power spectrum as

$\displaystyle
Q^{\rm(NL)}_{\mathcal{P}}\equiv\frac{\mathcal{P}_{\zeta}^{\rm(NL)}(k_{*})}{\mathcal{P}_{\zeta}(k_{*})}$
$\displaystyle=1+\left(\frac{3}{5}f_{\rm
NL}\right)^{2}\frac{k_{*}^{3}}{2\pi\mathcal{P}_{\zeta}({\bm{k}_{*}})}\int\text{d}^{3}q\frac{\mathcal{P}_{\zeta}({\bm{q}})\mathcal{P}(|{\bm{k}_{*}-\bm{q}}|)}{q^{3}|\bm{k}_{*}-\bm{q}|^{3}},$
(12)

where $k_{*}$ is the wavenumber at the peak of $\mathcal{P}_{\zeta}$.
$Q^{\rm(NL)}_{\mathcal{P}}$ is about $1.07$ for $r=0.5$ and $1.01$ for $r=1.0$
in our parameter set.

Finally, the power spectrum of the density perturbations is given by

$\displaystyle\mathcal{P}_{\delta}(k,t)$
$\displaystyle=\left(\frac{2}{3}\frac{k}{a(t)H(t)}\right)^{4}T(k\eta(t))^{2}\mathcal{P}_{\zeta}(k),$
(13)

where $\eta(t)$ is the conformal time and $T(x)$ is the transfer function
during radiation dominated era, which includes the suppression of the density
perturbation in sub-horizon as

$\displaystyle T(x)\equiv$
$\displaystyle~{}3\frac{\sin(x/\sqrt{3})-(x/\sqrt{3})\cos(x/\sqrt{3})}{(x/\sqrt{3})^{3}}.$
(14)

In this paper, we use the following set of model parameters:

$\displaystyle\lambda$ $\displaystyle=7.5\times 10^{-6},\quad v=5.77\times
10^{-2}{M_{\text{pl}}},\quad$ $\displaystyle g$ $\displaystyle=9.56\times
10^{-11},\quad\epsilon=2.57\times 10^{-10},$
$\displaystyle\begin{cases}&r=0.5,\quad\theta_{i}=5.3\times 10^{-2}\\\
&r=1.0,\quad\theta_{i}=6.5\times 10^{-2}\end{cases},$ (15)

where we take $r$ and $\theta_{i}$ to achieve the PBH abundance required to
explains the LIGO-Virgo event rate. Here we remark some relations between the
parameters and the shape of the power spectrum. The larger $r$ or smaller
$\theta_{i}$ increases the abundance of ALP as Eq. (10), and smaller
$\epsilon$ changes the dynamics of $\varphi_{0}$ at $t\sim t_{\rm pbh}$, both
of which result in larger amplitude of the power spectrum. The peak wavenumber
of the power spectrum is determined by the ratio of potential terms,
$g/(\lambda v^{2})$, as discussed in Eq. (3). The width of the spectrum
depends on how slowly $\varphi_{0}$ changes since the fluctuation of $\theta$
with mode $k$ depends on the field value $\varphi_{0}$ at the horizon crossing
as Eq. (5). We flatten the potential by choosing the smaller $v$, $g$ and
$\lambda$ compared to the previous paper [38] to achieve a broad power
spectrum. We numerically evaluate the dynamics of $\Phi$ assuming the chaotic
inflation whose potential is given by $V_{I}=m_{I}^{2}I^{2}/2$ with
$m_{\phi}\simeq 10^{13}\mathrm{GeV}$. The similar dynamics also holds for
other inflation models. We show the curvature power spectrum
$\mathcal{P}_{\zeta}(k)$ [Eq. (10)] in Fig. 1. We also show the constraints on
$\mu$–distortion [42] by COBE/FIRAS [43] and BBN [44], and our model safely
avoids the current constraints.

Figure 1:  The curvature power spectrum
$Q^{\rm(NL)}_{\mathcal{P}}\mathcal{P}_{\zeta}(k)$ based on Eqs.(10) and (12).
The non-Gaussian contribution is included by $Q^{\rm(NL)}_{\mathcal{P}}$. The
orange regions are constraints on the curvature power spectrum from
$\mu$–distortion [42] by COBE/FIRAS [43] and BBN [44].

##  3 PBH formation

The PBHs are formed by the collapse of high-density regions when they re-enter
the horizon. The criterion on the PBH formation is estimated by the numerical
simulation [45], in which the threshold value of the averaged density
fluctuation is obtained. Since the threshold value is too large to neglect the
nonlinear contribution, we need to use the effective threshold value including
the nonlinear relation between density and curvature perturbations. The
detailed analysis show that the effective threshold value is $\delta_{\rm
th(eff)}\simeq\sqrt{2}\times 0.53$ for the averaged linear density
perturbation defined by [46]

$\displaystyle\bar{\delta}_{R}(\bm{x})$
$\displaystyle\equiv\int\text{d}^{3}yW(|\bm{x}-\bm{y}|,R)\delta(\bm{y})$
$\displaystyle=\int\frac{\text{d}^{3}p}{(2\pi)^{3}}\tilde{W}(pR)e^{i\bm{p}\cdot\bm{x}}\delta_{p},$
(16)

where $R=k^{-1}$ is the scale of the horizon corresponding to the PBH mass
[Eq.(1)], $W(x,R)$ and $\tilde{W}(z)$ are the window functions in the real and
Fourier spaces, and the $\delta_{p}$ is the density fluctuation. Although it
is known that the choice of window function causes a large uncertainly [47], a
natural choice, often used in the literature, is the real-space top-hat window
function used in the threshold value in the numerical simulation,

$\displaystyle\tilde{W}(z)$ $\displaystyle=3\frac{\sin(z)-z\cos(z)}{z^{3}}.$
(17)

Thus, a PBH is formed when a region has the averaged density
$\bar{\delta}_{R}$ larger than the threshold value $\delta_{\rm th(eff)}$.

We estimate the PBH formation rate based on the Press–Schechter formalism,
where the PBH formation rate is calculated by the probability distribution of
the averaged density perturbations. In our model, $\bar{\delta}_{R}$ follows
the non-Gaussian distribution due to $f_{\rm NL}$ given by Eq. (11), which
drastically changes the PBH abundance. The probability distribution of
$\bar{\delta}_{R}$ is characterized by the variance and skewness defined by

$\displaystyle\sigma_{R}^{2}$
$\displaystyle=\Braket{\bar{\delta}_{R}^{2}}-\Braket{\bar{\delta}_{R}}^{2},$
(18) $\displaystyle\mu_{R}$
$\displaystyle=\sigma_{R}^{-3}\left(\Braket{\bar{\delta}_{R}^{3}}-3\Braket{\bar{\delta}_{R}^{2}}\Braket{\bar{\delta}_{R}}+2\Braket{\bar{\delta}_{R}}^{3}\right),$
(19)

where $\Braket{...}$ describes the ensemble average of $\delta_{p}$. For
simplicity, we neglect the scale dependence of $\mu_{R}$ and evaluate it at
$R=R_{\rm pbh}=k_{\rm pbh}^{-1}$, which corresponds to $30{M_{\odot}}$ PBHs.
Using the formula in [46], the skewness is numerically given by

$\displaystyle\mu\equiv\mu_{R}|_{R=R_{\rm pbh}}\simeq 3.39f_{\rm
NL}\sigma_{R}|_{R=R_{\rm pbh}}.$ (20)

We construct the statistical variable which reproduces the probability
distribution of $\bar{\delta}_{R}$. Using the Gaussian variable $\chi$
characterized by $\Braket{\chi^{2}}=\sigma_{R}^{2}$, we define the statistical
variable as

$\displaystyle\bar{\delta}[\chi]\equiv\chi+\frac{\mu}{6\sigma_{R}}(\chi^{2}-\sigma_{R}^{2}),$
(21)

which has the same variance and skewness in Eqs. (18) and (19) up to
$\mathcal{O}(f_{\rm NL}\sigma_{R})$. The probability distribution of
$\bar{\delta}[\chi]$ is given by [48, 49, 50]

$\displaystyle P_{R}^{\rm(NG)}(\bar{\delta})$
$\displaystyle=\sum_{i=\pm}\left|\frac{\text{d}\chi_{i}(\bar{\delta})}{\text{d}\bar{\delta}}\right|P_{R}^{\rm(G)}(\chi_{i}(\bar{\delta})),$
(22)

where $P_{R}^{\rm(G)}(\chi)$ is the Gaussian distribution function,

$\displaystyle
P_{R}^{\rm(G)}(\chi)=\frac{1}{\sqrt{2\pi}\sigma_{R}}\exp\left(-\frac{1}{2}\frac{\chi^{2}}{\sigma_{R}^{2}}\right),$
(23)

and $\chi_{\pm}(\bar{\delta})$ are two solutions of
$\bar{\delta}=\bar{\delta}[\chi]$,

$\displaystyle\chi_{\pm}(\bar{\delta})=\frac{3\sigma_{R}}{\mu}\left(-1\pm\sqrt{1+\frac{2\mu}{3}\left(\frac{\mu}{6}+\frac{\bar{\delta}}{\sigma_{R}}\right)}\right).$
(24)

The PBH formation rate is given by the probability of
$\bar{\delta}>\delta_{\rm th(eff)}$, which leads to

$\displaystyle\beta(R)$ $\displaystyle=\int_{\bar{\delta}>\delta_{\rm
th(eff)}}P_{R}^{\rm(NG)}(\bar{\delta})\text{d}\bar{\delta}=\int_{\bar{\delta}[\chi]>\delta_{\rm
th(eff)}}P_{R}^{\rm(G)}(\chi)\text{d}\chi$
$\displaystyle\simeq\frac{\sigma_{R}}{\sqrt{2\pi}\chi_{+}(\delta_{\rm
th(eff)})}\exp\left(-\frac{\chi_{+}(\delta_{\rm
th(eff)})^{2}}{2\sigma_{R}^{2}}\right).$ (25)

Here we have used $\mu>0$ and $\chi_{+}(\delta_{\rm th(eff)})/\sigma_{R}\gg 1$
in the last line. The present PBH abundance is given by

$\displaystyle f(M)\equiv\frac{\text{d}\Omega_{\text{PBH}}}{\text{d}\ln
M}\frac{1}{\Omega_{\text{DM}}}$ (26)
$\displaystyle=\frac{\beta(R(M))}{1.8\times
10^{-8}}\left(\frac{\gamma}{0.2}\right)^{3/2}\left(\frac{10.75}{\text{g}_{*}}\right)^{1/4}\left(\frac{0.12}{\Omega_{\text{DM}}h^{2}}\right)\left(\frac{M}{{M_{\odot}}}\right)^{-1/2},$
(27)

where $\text{g}_{*}$ is the number of relativistic degrees of freedom at
$T\sim 30\mathrm{MeV}$.

We show the calculated mass spectrum of PBH in Fig. 2. We also plot the
relevant constrains by microlensing experiments “MACHO/EROS/OGLE” [51, 52, 53]
and energy injection into CMB through accretion around PBHs [54]. Since our
model predicts the broad mass distribution, the large mass part of the
distribution could conflict with the accretion constraints. However, it is
noticed that accretion constraint has a large uncertainty. In fact, two
different constraints are obtained depending on assumptions as shown in Fig.
2. Our mass distribution is consistent if we adopt the weaker accretion
constraint.

Figure 2:  The mass spectrum of PBHs (blue line) and the constraints (orange
regions) by microlensing experiments “MACHO/EROS/OGLE”(solid line) [51, 52,
53], CMB through spherical accretion (dashed line) and disk accretion (dotted
line) [54].

##  4 Induced gravitational waves

The large density fluctuations induce the gravitational waves through the
nonlinear coupling $\Braket{\zeta\zeta h}$ when they re-enter the horizon. The
current energy fraction of GWs is written as

$\displaystyle\Omega_{\rm GW}(t_{0},k)$
$\displaystyle=\left(\frac{a_{c}^{2}H_{c}}{a_{0}^{2}H_{0}}\right)^{2}\Omega_{\rm
GW}(\eta_{c},k)$ $\displaystyle\simeq
0.83\left(\frac{\text{g}_{c}}{10.75}\right)^{-1/3}\Omega_{r,0}\Omega_{\rm
GW}(t_{c},k),$ (28)

where $\Omega_{r,0}$ is the current energy fraction of radiation, the
subscript “c” denotes values when GW production effectively finishes.

The energy density of the induced GWs at $t_{c}$ is calculated by solving the
equation of motion of GWs with the source term of scalar fluctuations, and it
is given by [30]

$\displaystyle\Omega_{\rm
GW}(t_{c},k)=\frac{8}{243}\int^{\infty}_{0}\text{d}y\int^{1+y}_{\left|1-y\right|}\text{d}x\mathcal{P}_{\zeta}(kx)\mathcal{P}_{\zeta}(ky)\frac{y^{2}}{x^{2}}$
$\displaystyle\qquad\times\left(1-\frac{(1+y^{2}-x^{2})^{2}}{4y^{2}}\right)^{2}\overline{\mathcal{I}(x,y,\eta_{c})}^{2},$
(29)

where overline means the time average over $\eta_{c}$.
$\mathcal{I}(x,y,\eta_{c})$ is written as

$\displaystyle\mathcal{I}(x,y,\eta_{c})=\frac{k^{2}}{a(\eta_{c})}\int^{\eta_{c}}\text{d}\bar{\eta}a(\bar{\eta})g_{k}(\eta_{c};\bar{\eta})f(ky,kx,\bar{\eta}).$
(30)

Here $g_{k}$ is the Green function,

$\displaystyle g_{k}(\eta,\tilde{\eta})$
$\displaystyle=\frac{\sin(k(\eta-\bar{\eta}))}{k}\theta(\eta-\bar{\eta}),$
(31)

and $f(k_{1},k_{2},\eta)$ is given by

$\displaystyle f(k_{1},k_{2},\eta)$
$\displaystyle=\bigg{[}2T(k_{1},\eta)T(k_{2},\eta)$ $\displaystyle+$
$\displaystyle\left(\frac{\dot{T}(k_{1},\eta)}{H(\eta)}+T(k_{1},\eta)\right)\left(\frac{\dot{T}(k_{2},\eta)}{H(\eta)}+T(k_{2},\eta)\right)\bigg{]},$
(32)

where $T(x)$ is the transfer function of the scalar fluctuations given by
Eq.(14).

We comment on the non-Gaussian contribution on the induced gravitational waves
[55, 40, 41]. It is pointed out in [41] that non-Gaussianity of scalar
fluctuations amplifies the induced gravitational waves when
$\mathcal{P}_{\zeta}f_{\rm NL}^{2}$ is large. Since $\mathcal{P}_{\zeta}f_{\rm
NL}^{2}<0.04$ in our calculation, the effect of non-Gaussianity is expected to
be sub-dominant. Thus, we approximately include the effect of the non-
Gaussianity on $\Omega_{\rm GW}$ by multiplying the factor
$(Q^{\rm(NL)}_{\mathcal{P}})^{2}$ given by Eq. (12). This approximation
includes a part of non-Gaussian contributions, “Hybrid” and “Reducible”-type
terms discussed in [40]. Hybrid type is a product of the Gaussian and non-
Gaussian contribution of curvature perturbation, and Reducible type is that of
non-Gaussian and non-Gaussian contribution. Although there are other types of
sources of GWs, it is known that Hybrid-type is one of the largest
contributions among them. Thus, our calculation can estimate most of the
effects of non-Gaussianity.

The estimated GW spectra for $r=0.5$ (blue line) and $r=1.0$ (orange line) are
shown in Fig. 3. The GW spectrum observed by the NANOGrav experiment is fitted
by the power-law spectrum around $f\sim 10^{-8}$,

$\displaystyle\Omega_{\rm
GW}(f)h^{2}=\frac{2\pi^{2}}{3}\frac{h^{2}f^{2}}{H_{0}^{2}}h_{c}^{2}(f)=A_{\Omega}f^{5-\gamma},$
(33)

where $\gamma$ is the tilt of the spectrum. In Fig. 3, we show the observed
GWs for $\gamma=5$ and 6 with 2-$\sigma$ uncertainty on $A_{\Omega}$. We also
plot current constraints by other PTA experiments, EPTA (solid) [56] and PPTA
(dotted)[57], and future sensitivity by SKA (dashed) [58]. It is found that
the broad power spectrum of GWs in the present model can explain the reported
NANOGrav signal.

Figure 3:  The induced GW spectrum and the constraints by PTA experiments. We
include the contribution of non-Gaussianity by using the factor
$Q^{\rm(NL)}_{\mathcal{P}}$. The orange lines are current constraints by EPTA
(solid) [56] and PPTA (Dotted)[57], and future sensitivity by SKA (Dashed)
[58]. We show the reported NANOGrav signal with $\gamma=5$ (blue region) and
$\gamma=6$ (pink region) with 2-$\sigma$ uncertainty (see Eq. (33)).

##  5 Conclusion

The reported signal by the NANOGrav experiment indicates the various
cosmological phenomena like cosmic string, phase transition and PBH formation.
The PBH formation scenario is attractive among them since the reported
frequency $f\sim 10^{-8}\mathrm{Hz}$ is close to the scale of the density
fluctuations to produce $30{M_{\odot}}$ PBH, which can explain the binary
black holes observed by LIGO-Virgo collaboration.

The typical PBH formation models predict the induced GWs with smaller
frequency and larger amplitude compared to the NANOGrav signal. To avoid this
difficulty, one needs to modify the induced GW spectrum by the broad power
spectrum and the non-Gaussianity of density fluctuations, which can enhance
the PBH formation rate and give a good fit to the NANOGrav signal. The axion-
like curvaton model can achieve the required features.

There are two types of the axion-like curvaton models; in type I the complex
field rolls down toward origin during inflation [35, 36, 37, 15], and in type
II it rolls down from the origin[38]. In this paper, we focused on the type II
model and chose the appropriate parameters of the potential term, which leads
to a broader power spectrum of the density perturbations than that in [38]. As
a result it was found that the induced GW spectrum can explain the NANOGrav
signal as shown in Fig. 3. Moreover, the broad power spectrum of the density
fluctuation results in the broad mass spectrum of PBHs as shown in Fig. 2,
which can be tested by the accumulation of the binary merger observations. Our
model predicts large local-type non-Gauusianity which can be probed through
observation of GWs. The spectrum of induced gravitational waves is a useful
tool to distinguish our model from others. For example, cosmic string and
type-I axion-like curvaton models gererally predict much broader GW spectra
than our model.

## Acknowledgment

This work was supported by JSPS KAKENHI Grant Nos. 17H01131 (M.K.), 17K05434
(M.K.), 20H05851 (M.K.), 21K03567(M.K.), JP19J21974 (H.N.), Advanced Leading
Graduate Course for Photon Science (H.N.), and World Premier International
Research Center Initiative (WPI Initiative), MEXT, Japan (M.K.).

## References

  * [1] NANOGrav collaboration, Z. Arzoumanian et al., _The NANOGrav 12.5-year Data Set: Search For An Isotropic Stochastic Gravitational-Wave Background_ , 2009.04496.
  * [2] J. Ellis and M. Lewicki, _Cosmic String Interpretation of NANOGrav Pulsar Timing Data_ , 2009.06555.
  * [3] S. Blasi, V. Brdar and K. Schmitz, _Has NANOGrav found first evidence for cosmic strings?_ , 2009.06607.
  * [4] W. Buchmuller, V. Domcke and K. Schmitz, _From NANOGrav to LIGO with metastable cosmic strings_ , _Phys. Lett. B_ 811 (2020) 135914, [2009.10649].
  * [5] R. Samanta and S. Datta, _Gravitational wave complementarity and impact of NANOGrav data on gravitational leptogenesis: cosmic strings_ , 2009.13452.
  * [6] S. Datta, A. Ghosal, R. Samanta and R. Sinha, _Baryogenesis from ultralight primordial black holes and strong gravitational waves_ , 2012.14981.
  * [7] N. Ramberg and L. Visinelli, _The QCD Axion and Gravitational Waves in light of NANOGrav results_ , 12, 2020, 2012.06882.
  * [8] Y. Nakai, M. Suzuki, F. Takahashi and M. Yamada, _Gravitational Waves and Dark Radiation from Dark Phase Transition: Connecting NANOGrav Pulsar Timing Data and Hubble Tension_ , 2009.09754.
  * [9] A. Addazi, Y.-F. Cai, Q. Gan, A. Marciano and K. Zeng, _NANOGrav results and Dark First Order Phase Transitions_ , 2009.10327.
  * [10] A. Neronov, A. Roper Pol, C. Caprini and D. Semikoz, _NANOGrav signal from MHD turbulence at QCD phase transition in the early universe_ , 2009.14174.
  * [11] S.-L. Li, L. Shao, P. Wu and H. Yu, _NANOGrav Signal from First-Order Confinement/Deconfinement Phase Transition in Different QCD Matters_ , 2101.08012.
  * [12] B. Barman, A. Dutta Banik and A. Paul, _Implications of NANOGrav results and UV freeze-in in a fast-expanding Universe_ , 2012.11969.
  * [13] V. Vaskonen and H. Veermäe, _Did NANOGrav see a signal from primordial black hole formation?_ , 2009.07832.
  * [14] V. De Luca, G. Franciolini and A. Riotto, _NANOGrav Hints to Primordial Black Holes as Dark Matter_ , 2009.08268.
  * [15] K. Inomata, M. Kawasaki, K. Mukaida and T. T. Yanagida, _NANOGrav results and LIGO-Virgo primordial black holes in axion-like curvaton model_ , 2011.01270.
  * [16] K. Kohri and T. Terada, _Solar-Mass Primordial Black Holes Explain NANOGrav Hint of Gravitational Waves_ , _Phys. Lett. B_ 813 (2021) 136040, [2009.11853].
  * [17] G. Domènech and S. Pi, _NANOGrav Hints on Planet-Mass Primordial Black Holes_ , 2010.03976.
  * [18] S. Sugiyama, V. Takhistov, E. Vitagliano, A. Kusenko, M. Sasaki and M. Takada, _Testing Stochastic Gravitational Wave Signals from Primordial Black Holes with Optical Telescopes_ , 2010.02189.
  * [19] M. Braglia, D. K. Hazra, F. Finelli, G. F. Smoot, L. Sriramkumar and A. A. Starobinsky, _Generating PBHs and small-scale GWs in two-field models of inflation_ , _JCAP_ 08 (2020) 001, [2005.02895].
  * [20] N. Bhaumik and R. K. Jain, _Stochastic induced gravitational waves and lowest mass limit of primordial black holes with the effects of reheating_ , 2009.10424.
  * [21] M. Braglia, X. Chen and D. Kumar Hazra, _Probing Primordial Features with the Stochastic Gravitational Wave Background_ , _JCAP_ 03 (2021) 005, [2012.05821].
  * [22] V. Atal, A. Sanglas and N. Triantafyllou, _NANOGrav signal as mergers of Stupendously Large Primordial Black Holes_ , 2012.14721.
  * [23] S. Vagnozzi, _Implications of the NANOGrav results for inflation_ , _Mon. Not. Roy. Astron. Soc._ 502 (2021) L11, [2009.13432].
  * [24] S. Bhattacharya, S. Mohanty and P. Parashari, _Implications of the NANOGrav result on primordial gravitational waves in nonstandard cosmologies_ , 2010.05071.
  * [25] S. Kuroyanagi, T. Takahashi and S. Yokoyama, _Blue-tilted inflationary tensor spectrum and reheating in the light of NANOGrav results_ , _JCAP_ 01 (2021) 071, [2011.03323].
  * [26] LIGO Scientific, Virgo collaboration, B. P. Abbott et al., _Binary Black Hole Population Properties Inferred from the First and Second Observing Runs of Advanced LIGO and Advanced Virgo_ , 1811.12940.
  * [27] S. Bird, I. Cholis, J. B. Muñoz, Y. Ali-Haïmoud, M. Kamionkowski, E. D. Kovetz et al., _Did LIGO detect dark matter?_ , _Phys. Rev. Lett._ 116 (2016) 201301, [1603.00464].
  * [28] S. Clesse and J. García-Bellido, _The clustering of massive Primordial Black Holes as Dark Matter: measuring their mass distribution with Advanced LIGO_ , _Phys. Dark Univ._ 15 (2017) 142–147, [1603.05234].
  * [29] M. Sasaki, T. Suyama, T. Tanaka and S. Yokoyama, _Primordial Black Hole Scenario for the Gravitational-Wave Event GW150914_ , _Phys. Rev. Lett._ 117 (2016) 061101, [1603.08338].
  * [30] K. Inomata, M. Kawasaki, K. Mukaida and T. T. Yanagida, _Double inflation as a single origin of primordial black holes for all dark matter and LIGO observations_ , _Phys. Rev._ D97 (2018) 043514, [1711.06129].
  * [31] B. J. Carr, _The Primordial black hole mass spectrum_ , _Astrophys. J._ 201 (1975) 1–19.
  * [32] R. Saito and J. Yokoyama, _Gravitational wave background as a probe of the primordial black hole abundance_ , _Phys. Rev. Lett._ 102 (2009) 161101, [0812.4339].
  * [33] R. Saito and J. Yokoyama, _Gravitational-wave constraints on the abundance of primordial black holes_ , _Progress of theoretical physics_ 123 (2010) 867–886.
  * [34] C. Unal, E. D. Kovetz and S. P. Patil, _Multi-messenger Probes of Inflationary Fluctuations and Primordial Black Holes_ , 2008.11184.
  * [35] M. Kawasaki, N. Kitajima and T. T. Yanagida, _Primordial black hole formation from an axionlike curvaton model_ , _Phys. Rev._ D87 (2013) 063519, [1207.2550].
  * [36] M. Kawasaki, N. Kitajima and S. Yokoyama, _Gravitational waves from a curvaton model with blue spectrum_ , _JCAP_ 1308 (2013) 042, [1305.4464].
  * [37] K. Ando, K. Inomata, M. Kawasaki, K. Mukaida and T. T. Yanagida, _Primordial Black Holes for the LIGO Events in the Axion-like Curvaton Model_ , 1711.08956.
  * [38] K. Ando, M. Kawasaki and H. Nakatsuka, _Formation of primordial black holes in an axionlike curvaton model_ , _Phys. Rev. D_ 98 (2018) 083508, [1805.07757].
  * [39] G. ’t Hooft, _Naturalness, chiral symmetry, and spontaneous chiral symmetry breaking_ , _NATO Sci. Ser. B_ 59 (1980) 135–157.
  * [40] C. Unal, _Imprints of Primordial Non-Gaussianity on Gravitational Wave Spectrum_ , _Phys. Rev. D_ 99 (2019) 041301, [1811.09151].
  * [41] R.-g. Cai, S. Pi and M. Sasaki, _Gravitational Waves Induced by non-Gaussian Scalar Perturbations_ , _Phys. Rev. Lett._ 122 (2019) 201101, [1810.11000].
  * [42] J. Chluba, A. L. Erickcek and I. Ben-Dayan, _Probing the inflaton: Small-scale power spectrum constraints from measurements of the cosmic microwave background energy spectrum_ , _The Astrophysical Journal_ 758 (2012) 76.
  * [43] D. J. Fixsen, E. S. Cheng, J. M. Gales, J. C. Mather, R. A. Shafer and E. L. Wright, _The Cosmic Microwave Background spectrum from the full COBE FIRAS data set_ , _Astrophys. J._ 473 (1996) 576, [astro-ph/9605054].
  * [44] K. Inomata, M. Kawasaki and Y. Tada, _Revisiting constraints on small scale perturbations from big-bang nucleosynthesis_ , _Phys. Rev._ D94 (2016) 043527, [1605.04646].
  * [45] M. Shibata and M. Sasaki, _Black hole formation in the Friedmann universe: Formulation and computation in numerical relativity_ , _Phys. Rev. D_ 60 (1999) 084002, [gr-qc/9905064].
  * [46] M. Kawasaki and H. Nakatsuka, _Effect of nonlinearity between density and curvature perturbations on the primordial black hole formation_ , _Phys. Rev. D_ 99 (2019) 123501, [1903.02994].
  * [47] K. Ando, K. Inomata and M. Kawasaki, _Primordial black holes and uncertainties on choice of window function_ , 1802.06393.
  * [48] S. Young and C. T. Byrnes, _Primordial black holes in non-Gaussian regimes_ , _JCAP_ 1308 (2013) 052, [1307.4995].
  * [49] C. T. Byrnes, E. J. Copeland and A. M. Green, _Primordial black holes as a tool for constraining non-Gaussianity_ , _Phys. Rev._ D86 (2012) 043512, [1206.4188].
  * [50] S. Young, D. Regan and C. T. Byrnes, _Influence of large local and non-local bispectra on primordial black hole abundance_ , _JCAP_ 1602 (2016) 029, [1512.07224].
  * [51] Macho collaboration, R. A. Allsman et al., _MACHO project limits on black hole dark matter in the 1-30 solar mass range_ , _Astrophys. J._ 550 (2001) L169, [astro-ph/0011506].
  * [52] EROS-2 collaboration, P. Tisserand et al., _Limits on the Macho Content of the Galactic Halo from the EROS-2 Survey of the Magellanic Clouds_ , _Astron. Astrophys._ 469 (2007) 387–404, [astro-ph/0607207].
  * [53] L. Wyrzykowski et al., _The OGLE View of Microlensing towards the Magellanic Clouds. IV. OGLE-III SMC Data and Final Conclusions on MACHOs_ , _Mon. Not. Roy. Astron. Soc._ 416 (2011) 2949, [1106.2925].
  * [54] P. D. Serpico, V. Poulin, D. Inman and K. Kohri, _Cosmic microwave background bounds on primordial black holes including dark matter halo accretion_ , _Phys. Rev. Res._ 2 (2020) 023204, [2002.10771].
  * [55] J. Garcia-Bellido, M. Peloso and C. Unal, _Gravitational Wave signatures of inflationary models from Primordial Black Hole Dark Matter_ , _JCAP_ 09 (2017) 013, [1707.02441].
  * [56] L. Lentati et al., _European Pulsar Timing Array Limits On An Isotropic Stochastic Gravitational-Wave Background_ , _Mon. Not. Roy. Astron. Soc._ 453 (2015) 2576–2598, [1504.03692].
  * [57] R. M. Shannon et al., _Gravitational waves from binary supermassive black holes missing in pulsar observations_ , _Science_ 349 (2015) 1522–1525, [1509.07320].
  * [58] C. J. Moore, R. H. Cole and C. P. L. Berry, _Gravitational-wave sensitivity curves_ , _Class. Quant. Grav._ 32 (2015) 015014, [1408.0740].


# Linear MIM-width of the Square of Trees

Svein Høgemo

###### Abstract

Graph parameters measure the amount of structure (or lack thereof) in a graph
that makes it amenable to being decomposed in a way that facilitates dynamic
programming. Graph decompositions and their associated parameters are
important both in practice (as a tool for designing robust algorithms for NP-
hard problems) and in theory (relating large classes of problems to the graphs
on which they are solvable in polynomial time).

Linear MIM-width is a variant of the graph parameter MIM-width, introduced by
Vatshelle. MIM-width is a parameter that is constant for many classes of
graphs. Most graph classes which have been shown to have constant MIM-width
also have constant linear MIM-width. However, computing the (linear) MIM-width
of graphs, or showing that it is hard, has proven to be a huge challenge. To
date, the only graph class with unbounded linear MIM-width, whose linear MIM-
width can be computed in polynomial time, is the trees. In this follow-up, we
show that for any tree $T$ with linear MIM-width $k$, the linear MIM-width of
its square $T^{2}$ always lies between $k$ and $2k$, and that these bounds are
tight for all $k$.

##  1 Introduction

Divide-and-conquer is a tried and true strategy for obtaining fast algorithms.
In addition to its classical use for polynomial-time algorithms, it is also a
valuable tool for designing parameterized algorithms for NP-hard problems. For
problems on graphs, this strategy often manifests in the form of _graph
decompositions_.

A graph decomposition is a tree structure that relates to the graph in such a
way that one can solve any of a number of problems on the graph by executing
dynamical programming on the decomposition. However, not all graphs will be
equally suitable for a given decomposition (nor should we expect them to be,
as we assume NP-hard problems are hard to solve on some instances). The
_width_ of a decomposition is a number $k$ that measures how badly a graph
fits the decomposition – the lower the $k$, the better the structure of the
graph fits the particular decomposition technique. The goal is to solve the
problem in time polynomial in the size of the graph, but exponential in $k$.
The arguably most famous decomposition technique, _tree decomposition_ , and
its associated parameter _treewidth_ works well on many classes of sparse
graphs – trees have treewidth 1, and in general graphs with low treewidth have
a “tree-like” structure. Due to the large role treewidth has played in
structural graph theory, and its usefulness in designing algorithms for a
plethora of problems [10, 19, 20, 22], the problem of computing tree
decompositions of decent treewidth has received great attention, and several
algorithms, both exact and approximating, have been proposed, with varying
usability in practice [4, 8, 9, 17]. The different uses of treewidth are well
explained in Bodlaender’s series of surveys on the matter [5, 6, 7, 18].

Since tree decomposition only works well for sparse graphs, several other
decomposition techniques have been proposed. One quite versatile technique is
_branch decomposition_ , where the vertices of the graph are mapped to the
leaves of a tree of maximum degree 3. Each edge in the tree corresponds to a
cut of the graph (defined by which side of the edge the leaf corresponding to
each vertex lies on), and there have been defined several width measures on
these decompositions, based on functions on the cuts. One such width measure
that has garnered significant theoretical interest is _MIM-width_. “MIM”
stands for “maximum induced matching”, and MIM-width measures the the biggest
induced matching in any of the bipartite graphs induced by the cuts defined by
the decomposition. MIM-width was introduced by Vatshelle [24] where its
strongness and algorithmic properties were expounded.

The significance of MIM-width lies in how veritably strong it is. Several
important graph classes have constant MIM-width, among them interval graphs
and circular arc graphs, permutation graphs, complements of graphs with
constant degeneracy, graphs with constant treewidth [24], and powers of graphs
with constant MIM-width [14]. The trade-off is a worse dependence on the
parameter than in most other types of decomposition – typical running times of
algorithms parameterized by MIM-width are $n^{O(k)}$ or $n^{O(k^{2})}$ for
graphs with MIM-width $k$, given a decomposition (see e.g. the algorithms
given in [3, 15, 16]). In other words, for any constant $k$ a polynomial time
algorithm exists, but the larger $k$ is, the higher the degree of the
polynomial is (this is called an $XP$ algorithm). Still, obtaining an
algorithm parameterized by MIM-width means that the problem is solvable in
polynomial time for many graph classes.

For all width measures of branch decompositions, one can define a “linear”
variant, where the allowed decompositions are restricted to linear layouts of
the graph. Linear graph parameters are interesting, both as a sort of test-bed
for proving things about the decompositions (as they are easier to reason
about), but also because algorithms that work on linear layouts are faster in
practice than algorithms that work on trees. Furthermore, MIM-width and linear
MIM-width have the peculiar relationship that most graph classes that have
constant MIM-width also have constant linear MIM-width. In fact, of the graph
classes that have unbounded treewidth or clique-width, only the class of leaf
powers and some of its subclasses, such as the trees themselves, are known to
have bounded MIM-width and unbounded linear MIM-width [13].

Without doubt, the greatest challenge regarding (linear) MIM-width is the
problem of computing it on arbitrary graphs. Several hardness results exist:
In [23], it was proved that it is at least as hard to compute the (linear)
MIM-width of a graph, or an optimal branch decomposition, as it is to compute
the MIM of a graph. Maximal induced matching is itself a hard problem in
several ways; in addition to being NP-complete, it is hard to approximate to a
constant factor in polynomial time [11], and also it is hard for the
parameterized complexity class $W[1]$, meaning that we likely will never find
any algorithm that decides whether an arbitrary graph has a maximal induced
matching of size at least $k$ in time bounded by $f(k)\cdot n^{O(1)}$ [21].
All of these barriers carry over to the problem of computing (linear) MIM-
width. This is coupled with a lack of positive results: To date, no algorithm
has been found that in polynomial time decides whether an arbitrary graph has
(linear) MIM-width at most $k$, for any constant $k$. Even for (linear) MIM-
width 1, no algorithm has been found that in polynomial time outputs a
decomposition of _any_ constant width, or concludes that the graph has width
$>1$. It is not implausible that recognizing graphs with (linear) MIM-width
$k$ is NP-complete for some (small) $k$, especially in light of a recent,
similar result regarding twin-width, another strong parameter [2]. If this
turns out to be true, it would naturally be a huge disadvantage, although not
taking away from its proven usefulness for the many graph classes with bounded
MIM-width.

Regarding positive results, for all graph classes with proven bounded (linear)
MIM-width, there is an easy way of finding a good layout or decomposition; for
example, an optimal layout of an interval graph is found by ordering the
vertices in order of the left endpoints of their respective intervals [1, 24].
Beware that for graph classes with a higher bound than 1, the decompositions
may not be optimal, they are only guaranteed to be bounded. To date, the only
polynomial-time algorithm for computing (linear) MIM-width on a class of
unbounded width is the one given in [12] for the linear MIM-width of trees
(the MIM-width of trees is 1). Om the class of trees, the linear MIM-width is
at most logarithmic and within a constant factor of their pathwidth (the
linear variant of treewidth).

This discrepancy between the use of MIM-width in designing robust graph
algorithms and how little we know about its computability, makes it hard to
gauge the exact potential of this parameter. Studying the linear MIM-width of
simple graph classes, like the squares of trees, can help illuminating what
makes the problem difficult to solve in the general case, and find out when
the structural tools we have for analysing trees break down. Indeed, if
deciding (linear) MIM-width $\leq k$ for some constant $k$ actually proves NP-
complete, it will likely be proved by restricting the problem to some graph
class on which the problem is also NP-complete. This is not unusual, as
restricting the problem to a special case makes reductions easier to find (see
e.g. the classic result for the the NP-completeness of chordal completion
[26]).

Following the result in [14], any power of a graph $G^{k}$ has at most twice
the (linear) MIM-width of the original graph $G$. Therefore powers of trees
have at most logarithmic linear MIM-width. This is in stark contrast to e.g.
pathwidth, since the square of a star is a complete graph and therefore has
linear pathwidth. Furthermore, for any graph $G$, $G^{diam(G)}$ is a complete
graph and therefore has linear MIM-width 1. Squares of trees are a simple
graph class that nevertheless have more going on than trees themselves, and
therefore serves as a natural class to extend the research from [12] in.

The rest of the paper is organized as follows: Section 2 introduces the
notation necessary to follow the paper; Section 3 contains the main result and
the previous results that we use in order to prove the main result; and
finally, Section 4 concludes with a short discussion on its implications (if
any).

##  2 Preliminaries

All graphs considered are finite and simple. $V(G)$ and $E(G)$ denotes the
sets of vertices and edges in the graph $G$, respectively.

$N_{G}(v)$ denotes the _open neighborhood_ of the vertex $v$ in the graph $G$,
i.e. the set of vertices with which $v$ shares an edge. $N_{G}[v]$ denotes the
_closed neighborhood_ , i.e. $N(v)\cup\\{v\\}$. For a set of vertices
$S\subseteq V(G)$, $N[S]=\bigcup_{v\in S}N[v]$, and $N(S)=N[S]\setminus S$. If
$G$ is obvious from context, subscripts can be omitted.

Let $S,T$ be two disjoint subsets of $V(G)$. $G[S,T]$ is the bipartite graph
induced by $S$ and $T$, i.e. the subgraph consisting of all edges with one
endpoint in $S$ and the other in $T$.

If $T$ is a rooted tree, the notation $T[v]$ for some node $v$ refers to the
subtree rooted in $v$, i.e. the rooted tree consisting of $v$ (as root) and
all its descendants.

The _distance_ between two vertices $u,v\in V(G)$ is the length of a shortest
path between $u$ and $v$. Given two subgraphs $H,H^{\prime}\subseteq G$, the
distance between $H$ and $H^{\prime}$ is the minimum distance between any
vertex in $H$ and any vertex in $H^{\prime}$. For example, the distance
between $H$ and $H^{\prime}$ is at least 1 iff $H$ and $H^{\prime}$ are
disjoint. The _diameter_ of $G$, $diam(G)$, is the length of the longest
distance between any two vertices in $G$.

######  Definition 1 (graph power).

Given a graph $G$ and some integer $k\geq 1$, the $k$-th _power_ of $G$,
denoted $G^{k}$, is the unique graph that has $V(G^{k})=V(G)$, and the
additional property that for any two vertices $u,v$, $u$ and $v$ are adjacent
in $G^{k}$ if and only if they have distance at most $k$ in $G$. If we
consider every vertex in $G$ “adjacent to” itself, we have that $G^{1}=G$, and
furthermore, given the adjacency matrix of $G$ as a boolean matrix
$\mathcal{A}$, the adjacency matrix of $G^{k}$ is given by $\mathcal{A}^{k}$.

We define the _square_ of $G$ as $G^{2}$, the second power of $G$.

######  Definition 2 (linear layout).

A _linear layout_ of a graph $G$ is a bijection
$\sigma:V(G)\rightarrow\\{1,\ldots,|V(G)|\\}$, i.e. a total order on the
vertices of $G$. We will routinely use $v_{i}$ as a shorthand for
$\sigma^{-1}(i)$ when $G$ and $\sigma$ is given.

The next three definitions are found in [12]:

######  Definition 3 (maximum induced matching, mim-width of a layout).

For a graph $G$ on $n$ vertices, we denote by $\mathrm{mim}(G)$ the size of
its _maximum induced matching_ (MIM), the largest number of edges whose
endpoints induce a matching. Let $\sigma$ be a linear layout of $G$. For any
index $1\leq i<n$ we have a subset of $V(G)$,
$V^{\sigma}_{i}=\\{v_{1},\ldots,v_{i}\\}$. We call the partition
$(V^{\sigma}_{i},\overline{V^{\sigma}_{i}})$ a _cut_ of $G$. The _maximum
induced matching width_ , or _MIM-width_ of G under layout $\sigma$ is denoted
$\mathrm{mw}(\sigma,G)$, and is defined as the maximum, over all $1\leq i<n$,
of $mim(G[V^{\sigma}_{i},\overline{V^{\sigma}_{i}}])$.

######  Definition 4 (linear MIM-width).

The _linear maximum induced matching width_ – _linear MIM-width_ – of $G$ is
denoted $\mathrm{lmw}(G)$, and is the minimum value of $\mathrm{mw}(\sigma,G)$
over any possible layout $\sigma$ of the vertices of $G$.

Note: When talking about a cut $(V^{\sigma}_{i},\overline{V^{\sigma}_{i}})$ in
a graph $G$, a vertex $v$ is said to lie _to the left of_ the cut if and only
if $v\in V^{\sigma}_{i}$; otherwise, $v$ is said to lie _to the right of_ the
cut.

######  Definition 5 ($k$-neighbor).

Let $x$ be a node in the tree $T$ and $v$ a neighbor of $x$. If $v$ has a
neighbor $u\neq x$ such that the component of $T\setminus vu$ containing $u$
has linear MIM-width at least $k$, then we call $v$ a $k$-_neighbor_ of $x$.

##  3 Results

We repeat two structural results from [12] here:

######  Lemma 1 (path layout lemma ([12], Lemma 1)).

Let $T$ be a tree. If there exists a path $P=x_{1},\ldots,x_{p})$ in $T$ such
that every connected component of $T\setminus N[P]$ has linear MIM-width $\leq
k$, then $lmw(T)\leq k+1$.

######  Lemma 2 ([12], Theorem 1).

Let $T$ be a tree. $lmw(T)\geq k+1$ if and only if there is a node $x\in T$
that has at least three $k$-neighbors.

Lemma 2 will be our main tool for recursively generating trees with a certain
linear MIM-width. The next lemma is a generalization of its backwards
direction:

######  Lemma 3 (due to Vågset [25]).

Let $G$ be a graph, and let $C_{1}$, $C_{2}$ and $C_{3}$ be connected induced
subgraphs of $G$ with pairwise distance at least two. Let $k$ be the minimum
linear MIM-width of these three subgraphs. If, for each pair of subgraphs
$C_{i},C_{j}$ there exists a path $P_{i,j}$ that runs from $C_{i}$ to $C_{j}$
without intersecting the closed neighborhood of the third subgraph, then the
linear MIM-width of $G$ is strictly greater than $k$.

###### Proof.

To prove this lemma, we assume towards a contradiction that there exists a
linear layout $\sigma$ of $G$ with MIM-width $k$.

By definition of linear MIM-width, $\sigma$ must contain three cuts
$V^{\sigma}_{i_{1}}$, $V^{\sigma}_{i_{2}}$ and $V^{\sigma}_{i_{3}}$, such that
$G[V^{\sigma}_{i_{1}},\overline{V^{\sigma}_{i_{1}}}]$ contains an induced
matching $M_{1}$ with $k$ edges from $C_{1}$; likewise there exists an $M_{2}$
of size $k$ in $C_{2}$, and $M_{3}$ in $C_{3}$. Since the subgraphs have
distance at least 2, any edge from (say) $C_{2}$ can increase the size of the
matchings $M_{1}$ or $M_{3}$. Assuming that $(G,\sigma)$ has MIM-width $k$, it
must thus be the case that either, for every vertex $x\in C_{2}$,
$\sigma(x)\leq i_{1}$, or that, for every vertex $x\in C_{2}$,
$\sigma(x)>i_{1}$. These facts are obviously also true of every other pair of
subgraphs.

Now we assume w.l.o.g. that $i_{1}<i_{2}<i_{3}$. This implies that some vertex
in $C_{1}$ lies to the left of the cut $V^{\sigma}_{i_{2}}$ and that some
vertex in $C_{3}$ lies to the right of the cut. From the previous fact, we can
directly infer that _all_ vertices of $C_{1}$ (resp. $C_{3}$) lie to the left
(resp. right) of $V^{\sigma}_{i_{2}}$. This means that some edge $e$ on the
path $P_{1,3}$ must cross the cut. But the path has also distance at least 2
to $C_{2}$, thus $e$ can be taken into $M_{2}$; this implies that
$mw(G,\sigma)\geq k+1$. By contradiction, the above lemma is true.∎∎

The last lemma we will use is a special case of Theorem 5 of [14], stated in
terms of linear layouts; in the original paper this result is given in terms
of branch decompositions.

######  Lemma 4.

Given a graph $G$ and a layout $\sigma$ such that $G$ has MIM-width $k$ under
$\sigma$, then for any power $G^{m}$ of $G$, $G^{m}$ has MIM-width at most
$2k$ under $\sigma$. Therefore, $lmw(G^{m})\leq 2\cdot lmw(G)$.

###### Proof.

This follows directly from [14], Theorem 5, that states that the property
holds for any cut of the graph. Therefore it must hold also for linear
layouts.∎∎

######  Theorem 1.

For any tree $T$ with $lmw(T)=k$, then $k\leq lmw(T^{2})\leq 2k$. These bounds
are tight for any $k\geq 0$.

###### Proof.

To prove the first inequality, we use induction to prove that for every tree
$T$, if $lmw(T)\geq k$, then $lmw(T^{2})\geq k$.

For the base case, it is trivial to see that if $T$ has linear MIM-width at
least 1 (that is, if $T$ contains at least one edge), then $T^{2}$ must also
have linear MIM-width at least 1.

For the inductive step, we fix a $k\geq 1$, and assume that for every tree $T$
with $lmw(T)\geq k$, $lmw(T^{2})\geq k$. We show that for any tree
$T^{\prime}$ with $lmw(T^{\prime})\geq k+1$, $lmw(T^{\prime 2})\geq k+1$ as
follows:

Let $T$ be a tree with $lmw(T)\geq k+1$. From Lemma 2, we know that there must
exist a node $x\in T$ that has at least 3 $k$-neighbors in $T$, i.e. $x$ has
neighbors $v_{1},v_{2},v_{3}$ such that there are three subtrees
$S_{1},S_{2},S_{3}\subseteq T\setminus N[x]$ adjacent to $v_{1}$, $v_{2}$ and
$v_{3}$ respectively, with $lmw(S_{1}),lmw(S_{2}),lmw(S_{3})\geq k$.

By the inductive assumption, $lmw(S^{2}_{1}),lmw(S^{2}_{2}),lmw(S^{2}_{3})\geq
k$. $S^{2}_{1}$, $S^{2}_{2}$ and $S^{2}_{3}$ are three connected induced
subgraphs of $T^{2}$, all of distance at least two from each other.
Furthermore, between each two of the subtrees – say $S^{2}_{1}$ and
$S^{2}_{2}$ – there exists a path in $T^{2}$ that does not intersect the
closed neighborhood of the third subtree, in this case $N[S^{2}_{3}]$. By
Lemma 3, $T^{2}$ must have linear MIM-width at least $k+1$.

The second inequality follows directly from Lemma 4.

Next, we show the downward tightness of the bound; that is, that there exists
an infinite family of trees

$\mathcal{L}=(L(0),L(1),\ldots)$

where $lmw(L(k)^{2})=lmw(L(k))=k$ for every $k\geq 0$. For ease of notation,
we will define each tree in $\mathcal{L}$ as a rooted tree.

$L(0)$ is defined to be the singleton $K_{1}$. For every $k\geq 0$, $L(k+1)$
has a root $u$ with three children, $v_{1},v_{2},v_{3}$. Each $v_{i}$ in turn
has one child that is the root of a copy of $L(k)$, which we call $S_{i}$.
This recursive structure enables us to show that between any two trees $L(k)$
and $L(k+1)$, the linear MIM-width must increase with at least 1, due to Lemma
2. On the other hand, between any $L(k)^{2}$ and $L(k+1)^{2}$, the linear MIM-
width must increase with _at most_ 1\. This is shown by constructing a layout
of $L(k+1)^{2}$ of MIM-width $k+1$; this layout is identical to the one given
in the proof of the Path Layout Lemma (see [12] for detatils).

We prove the following claim by induction: For any tree $L(k)\in\mathcal{L}$,
$lmw(L(k)^{2})=lmw(L(k))=k$.

The base case is the trivial observation that $lmw(K^{2}_{1})=lmw(K_{1})=0$.

For the inductive step, we assume that $lmw(L(k)^{2})=lmw(L(k))=k$ for some
$k\geq 0$, and show that $lmw(L(k+1)^{2})=lmw(L(k+1))=k+1$. From the struc-
ture of $L(k+1)$ and the induction hypothesis, it is evident that the root $u$
has three k-neighbors, namely all its children.

By Lemma 2, we can conclude that $lmw(L(k+1))\geq k+1$. Regarding
$L(k+1)^{2}$, by the induction hypothesis there exists a layout of $L(k)^{2}$
that has MIM-width exactly $k$. We thus have optimal layouts
$\sigma_{S^{2}_{1}},\sigma_{S^{2}_{2}},\sigma_{S^{2}_{3}}$ available. We
construct a layout $\sigma_{L(k+1)}$ that has MIM-width exactly $k+1$ as
follows:

$\sigma_{L(k+1)^{2}}=(u)\oplus\sigma_{S^{2}_{1}}\oplus(v_{1})\oplus\sigma_{S^{2}_{2}}\oplus(v_{2})\oplus\sigma_{S^{2}_{3}}\oplus(v_{3})$

where $\oplus$ signifies concatenation.

For any cut $(V^{\sigma}_{i},\overline{V^{\sigma}_{i}})$, a maximum matching
contains at most $k$ edges from within some $S^{2}_{a}$. How many edges from
outside $S^{2}_{a}$, i.e. in the graph

$G^{\prime}_{i}:=L(k+1)^{2}[V^{\sigma}_{i},\overline{V^{\sigma}_{i}}]-E(S^{2}_{i})$

can be taken into a matching? We see that every vertex in
$V(G^{\prime}_{i})\cap V^{\sigma}_{i}$ is in
$\\{u\\}\cup\\{v_{1},\ldots,v_{a-1}\\}\cup(S_{a}\cap V^{\sigma}_{i})$, or has
no neighbors. Every vertex in $S_{a}\cap V^{\sigma}_{i}$ only has at most
$v_{a}$ as neighbor, $v_{1},\ldots,v_{a-1}$ have $v_{a},\ldots,v_{3}$ as
neighbors, and finally $u$ has all of these and possibly also one vertex in
each of $S_{a},\ldots,s_{3}$ as neighbors. This implies that $G^{\prime}_{i}$
is a bipartite chain and has MIM 1. Thus, no induced matching in
$L(k+1)^{2}[V^{\sigma}_{i},\overline{V^{\sigma}_{i}}]$ can have size more than
$k+1$.

Now we have that $lmw(L(k+1)^{2})\leq k+1\leq lmw(L(k+1))$. But, as we have
proven above, $lmw(L(k+1)^{2})\geq lmw(L(k+1))$. Thus,
$lmw(L(k+1)^{2})=lmw(L(k+1))=k+1$, and every tree in $\mathcal{L}$ has the
same linear MIM-width as its square.

Finally, we show the upward tightness of the bound; that is, that that there
exists an infinite family of trees $\mathcal{H}=(H(0),H(1),\ldots)$ where
$lmw(H(k)^{2})=2\cdot lmw(H(k))=2k$ for every $k\geq 0$. We will also define
each tree in $\mathcal{H}$ as a rooted tree.

$H(0)$ is again the singleton $K_{1}=u_{0}$. For every $k\geq 0$, $H(k+1)$ has
a root $u_{k}$ with three children, $v_{1},v_{2},v_{3}$. Each $v_{i}$ in turn
has three children, each of which is the root of a copy of $H(k)$. We call
these trees $S_{i,1}$, $S_{i,2}$ and $S_{i,3}$. This recursive structure
enables us to show that between any two trees $H(k)$ and $H(k+1)$, the linear
MIM-width must increase with at most 1, due to Lemma 1: Taking the path to be
$(u_{k+1})$, we see that all the subtrees that remain after removing the
neighborhood of $u_{k+1}$ are the $S_{i,a}$ for $1\leq i,a\leq 3$. Since, by
assumption, $S_{i,a}$ has linear MIM-width $k$, $lmw(H(k+1))\leq k+1$. And in
fact, it is exactly $k+1$ since $H(k)$ is a supertree of $L(k)$ for every $k$.

On the other hand, between any $H(k)^{2}$ and $H(k+1)^{2}$, the linear MIM-
width must increase with at least 2. This is shown by applying Lemma 3 twice.
To this end, we must show that the linear MIM-width of $H(k)^{2}$ does not
decrease when removing its root $u_{k}$; this trick is to ensure that the
subgraphs we consider are situated far enough apart that Lemma 3 is
applicable.

We will now prove the following claim by induction: Given that $uk$ is the
root of some tree $H(k)$,
$lmw(H(k)^{2})=lmw(H(k)^{2}\setminus\\{u_{k}\\})=2\cdot lmw(H(k))=2k$. Note
that the graph $H(k)^{2}\setminus\\{u_{k}\\}$ is still a connected graph.

For the base case, it is clear that
$lmw(k^{2}_{1})=lmw(u_{0})=lmw(\emptyset)=2\cdot lmw(K_{1})=0$.

For the induction step, we assume that
$lmw(H(k)^{2})=lmw(H(k)^{2}\setminus\\{u_{k}\\})=2\cdot lmw(H(k))=2k$ for some
tree $H(k)$ with root $u_{k}$, and show that
$lmw(H(k)^{2})=lmw(H(k+1)^{2}\setminus\\{u_{k+1}\\})=2\cdot
lmw(H(k+1)^{2})=2(k+1)$. We know from the induction hypothesis that for every
$(S_{i,a})^{2}$, the graph $S^{\prime}_{i,a}=(S_{i,a})2\setminus\\{u_{k}\\}$
is a connected subgraph of $H(k)^{2}$ with linear MIM-width $2k$. Thus, the
three graphs $S^{\prime}_{i,1}$, $S^{\prime}_{i,2}$ and $S^{\prime}_{i,3}$ are
three connected subgraphs with pairwise distance two in the graph
$T^{\prime}_{i}=(H(k+1)[v_{i}])^{2}\setminus\\{v_{i}\\}$. Furthermore, for
each pair of sub-graphs $S^{\prime}_{i,a}$ and $S^{\prime}_{i,b}$, there is a
path between them that does not intersect the neighborhood of $S_{i,c}$ (the
red path in the illustration below). By Lemma 3, the linear MIM-width of every
$T^{\prime}_{i}$ is at least $2k+1$. (In the case $k=0$, the notion of paths
between $S^{\prime}_{i,a}$ and $S^{\prime}_{i,b}$, which are empty sets, does
not really make sense. In this case, just note that $T^{\prime}_{i}$ contains
at least one edge and thus must have linear MIM-width at least 1.)

We use the same argument one more time: $T_{1}$, $T^{\prime}_{2}$ and
$T^{\prime}_{3}$ are three connected subgraphs with pairwise distance two in
the graph $H^{\prime}:=H(k+1)^{2}\setminus\\{u_{k+1}\\}$. For each pair
$T^{\prime}_{a}$ and $T^{\prime}_{b}$, there is a path between them that does
not intersect the neighborhood of the third subgraph $T^{\prime}_{c}$ (the
blue path in the illustration below). Thus, $lmw(H^{\prime})\geq 2k+2$. Since
$H^{\prime}$ is an induced subgraph of $H(k+1)^{2}$, we have the situation
that

$lmw(H(k+1)^{2})\geq lmw(H^{\prime})\geq 2(k+1)=2\cdot lmw(H(k+1))$

But, as we have noted above, $lmw(H(k+1)^{2})\leq 2\cdot lmw(H(k+1))$. Thus,

$lmw(H(k+1)^{2})=lmw(H(k+1)^{2}\setminus\\{u_{k+1}\\})=2\cdot
lmw(H(k+1))=2(k+1)$

and every tree in $\mathcal{H}$ has half the linear MIM-width of its square.∎∎

Figure 1: A subgraph of the graph $H(k+1)^{2}$. Dashed lines indicate power
edges.

##  4 Conclusion

We have shown that there is little connection between the linear MIM-width of
a tree and that of its square, except the fact that it cannot decrease. This
fact is interesting, as it implies that taking the square of a tree does not
makes its vertices more well-connected than in the original tree. As we know,
there must, for any tree $T$ with $lmw(T)\geq 2$, be an exponent $k$ such that
$lmw(T^{k})<lmw(T)$, since for any finite graph $G$, $G^{diam(G)}$ is a
complete graph and thus has linear MIM-width 1. How high this exponent must be
to decrease the linear MIM-width (or bring it down to 1), and whether there
exists a poly-time algorithm to evaluate the linear MIM-width of powers of
trees, must still be left in the open.

##  5 Acknowledgements

The author would like to thank O-joung Kwon for the initial discussion of this
topic.

## References

  * [1] Belmonte, R., Vatshelle, M.: Graph classes with structured neighborhoods and algorithmic applications. Theoretical Computer Science 511, 54–65 (2013)
  * [2] Bergé, P., Bonnet, É., Déprés, H.: Deciding twin-width at most 4 is np-complete. In: 49th EATCS International Colloquium on Automata, Languages and Programming (ICALP 2022) (2022)
  * [3] Bergougnoux, B., Kanté, M.M.: More applications of the d-neighbor equivalence: Acyclicity and connectivity constraints. SIAM Journal on Discrete Mathematics 35(3), 1881–1926 (2021)
  * [4] Bodlaender, H.L.: A linear time algorithm for finding tree-decompositions of small treewidth. In: Proceedings of the twenty-fifth annual ACM symposium on Theory of computing. pp. 226–234 (1993)
  * [5] Bodlaender, H.L.: Treewidth: Algorithmic techniques and results. In: International Symposium on Mathematical Foundations of Computer Science. pp. 19–36. Springer (1997)
  * [6] Bodlaender, H.L.: Treewidth: characterizations, applications, and computations. In: International Workshop on Graph-Theoretic Concepts in Computer Science. pp. 1–14. Springer (2006)
  * [7] Bodlaender, H.L.: Treewidth: Structure and algorithms. In: International Colloquium on Structural Information and Communication Complexity. pp. 11–25. Springer (2007)
  * [8] Bodlaender, H.L., Koster, A.M.: Treewidth computations i. upper bounds. Information and Computation 208(3), 259–275 (2010)
  * [9] Bodlaender, H.L., Koster, A.M.: Treewidth computations ii. lower bounds. Information and Computation 209(7), 1103–1119 (2011)
  * [10] De Givry, S., Schiex, T., Verfaillie, G.: Exploiting tree decomposition and soft local consistency in weighted csp. In: AAAI. vol. 6, pp. 1–6 (2006)
  * [11] Elbassioni, K., Raman, R., Ray, S., Sitters, R.: On the approximability of the maximum feasible subsystem problem with 0/1-coefficients. In: Proceedings of the 2009 Annual ACM-SIAM Symposium on Discrete Algorithms (SODA). pp. 1210–1219 (2009)
  * [12] Høgemo, S., Telle, J.A., Vågset, E.R.: Linear mim-width of trees. In: Graph-Theoretic Concepts in Computer Science: 45th International Workshop, WG 2019, Vall de Núria, Spain, June 19–21, 2019, Revised Papers 45. pp. 218–231. Springer (2019)
  * [13] Jaffke, L.: Bounded Width Graph Classes in Parameterized Algorithms. Ph.D. thesis, University of Bergen (2020)
  * [14] Jaffke, L., Kwon, O.j., Strømme, T.J., Telle, J.A.: Mim-width iii. graph powers and generalized distance domination problems. Theoretical Computer Science 796, 216–236 (2019)
  * [15] Jaffke, L., Kwon, O.j., Telle, J.A.: A unified polynomial-time algorithm for feedback vertex set on graphs of bounded mim-width. In: 35th Symposium on Theoretical Aspects of Computer Science. vol. 7148, p. 23 (2018)
  * [16] Jaffke, L., Kwon, O.j., Telle, J.A.: Mim-width i. induced path problems. Discrete Applied Mathematics 278, 153–168 (2020)
  * [17] Korhonen, T.: A single-exponential time 2-approximation algorithm for treewidth. In: 2021 IEEE 62nd Annual Symposium on Foundations of Computer Science (FOCS). pp. 184–192. IEEE (2022)
  * [18] Koster, A.M., Bodlaender, H.L., Van Hoesel, S.P.: Treewidth: computational experiments. Electronic Notes in Discrete Mathematics 8, 54–57 (2001)
  * [19] Madani, R., Sojoudi, S., Fazelnia, G., Lavaei, J.: Finding low-rank solutions of sparse linear matrix inequalities using convex optimization. SIAM Journal on Optimization 27(2), 725–758 (2017)
  * [20] Maniu, S., Senellart, P., Jog, S.: An experimental study of the treewidth of real-world graph data. In: 22nd International Conference on Database Theory (ICDT 2019). Schloss Dagstuhl-Leibniz-Zentrum fuer Informatik (2019)
  * [21] Moser, H., Sikdar, S.: The parameterized complexity of the induced matching problem. Discrete Applied Mathematics 157(4), 715–727 (2009)
  * [22] Ordyniak, S., Szeider, S.: Parameterized complexity results for exact bayesian network structure learning. Journal of Artificial Intelligence Research 46, 263–302 (2013)
  * [23] Sæther, S.H., Vatshelle, M.: Hardness of computing width parameters based on branch decompositions over the vertex set. Theoretical Computer Science 615, 120–125 (2016)
  * [24] Vatshelle, M.: New width parameters of graphs. Ph.D. thesis, The University of Bergen (2012)
  * [25] Vågset, E.R.: Unpublished result. (2018)
  * [26] Yannakakis, M.: Computing the minimum fill-in is np-complete. SIAM Journal on Algebraic Discrete Methods 2(1), 77–79 (1981)