# Wav-KAN: Wavelet Kolmogorov-Arnold Networks

Zavareh Bozorgasl, , Hao Chen Z. Bozorgasl<EMAIL_ADDRESS>and H. Chen<EMAIL_ADDRESS>are with the Department of Electrical and
Computer Engineering, Boise State University, Boise, ID, 83712.

###### Abstract

In this paper 111Another title for the paper can be ”Wavelet for Everybody”.,
we introduce Wav-KAN 222Wav-KAN refers to a family of continuous and discrete
wavelet transforms in KAN., an innovative neural network architecture that
leverages the Wavelet Kolmogorov-Arnold Networks (Wav-KAN) framework to
enhance interpretability and performance. Traditional multilayer perceptrons
(MLPs) and even recent advancements like Spl-KAN [1] face challenges related
to interpretability, training speed, robustness, computational efficiency, and
performance. Wav-KAN addresses these limitations by incorporating wavelet
functions into the Kolmogorov-Arnold network structure, enabling the network
to capture both high-frequency and low-frequency components of the input data
efficiently. Wavelet-based approximations employ orthogonal or semi-orthogonal
basis and also maintains a balance between accurately representing the
underlying data structure and avoiding overfitting to the noise. While
continuous wavelet transform (CWT) has a lot of potentials, we also employed
discrete wavelet transform (CWT) for multiresolution analysis which obviated
the need for recalculation of the previous steps in finding the details 333We
demonstrate its efficacy in capturing detailed signal characteristic. and
efficiently combines local detailed information where the data points are
dense with broader trends where the data points are sparse. Analogous to how
water conforms to the shape of its container, Wav-KAN adapts to the data
structure, resulting in enhanced accuracy, faster training speeds, and
increased robustness compared to Spl-KAN and MLPs. Our results highlight the
potential of Wav-KAN as a powerful tool for developing interpretable and high-
performance neural networks, with applications spanning various fields. This
work sets the stage for further exploration and implementation of Wav-KAN in
frameworks such as PyTorch, TensorFlow, and also it makes wavelet in KAN in
wide-spread usage like nowadays activation functions like ReLU and sigmoid in
universal approximation theory (UAT). The codes to replicate the simulations
are available at https://github.com/zavareh1/Wav-KAN.

###### Index Terms:

Kolmogorov-Arnold Networks (KAN), Wavelet, Wav-KAN, Neural Networks.

##  I Introduction

Advancements in artificial intelligence (AI) have led to the creation of
highly proficient AI systems that make decisions for reasons that are not
clear to us. This has raised concerns about the widespread deployment of
untrustworthy AI systems in the economy and our daily lives, introducing
several new risks, including the potential for future AIs to deceive humans to
achieve undesirable objectives [2, 3]. The tendency to unravel the black-box
behaviour of neural network has attracted a lot attention in recent years.
Interpretability of neural networks is crucial as it influences trust in these
systems and helps address ethical concerns such as algorithmic discrimination.
It is also essential for applying neural networks in scientific fields like
drug discovery and genomics, where understanding the model’s decisions is
necessary for validation and regulatory compliance [4, 5].
The multilayer feedforward perceptron (MLP) model is among the most widely
used and practical neural network models [6]. Despite their wide-spread usage,
MLPs have serious drwabacks like consuming almost all non-embedding parameters
in transformers and also have less interpretability relative to attention
layers [7, 1, 8].

Model renovation approaches aim to enhance interpretability by incorporating
more understandable components into a network. These components can include
neurons with specially designed activation functions, additional layers with
specific functionalities, modular architectures, and similar features [9].

Breaking down a neural network into individual neurons has a key challenge of
polysemantic neurons, i.e., each neuron activates for several unrelated types
of feature [10]. [11] explores the reasons behind the emergence of
polysemanticity and propose that it occurs because of superposition, i.e.,
models learn more distinct features than the available dimensions in a layer.
Due to the limitation that a vector space can only contain as many orthogonal
vectors as it has dimensions, the network ends up learning an overcomplete set
of non-orthogonal features. For superposition to be beneficial, these features
must activate sparsely; otherwise, the interference between non-orthogonal
features would negate any performance improvements [8].

Kolmogorov-Arnold Networks (KANs) which stand on the Kolmogorov-Arnold
representation theorem [12] can bring a lot of advantages including more
interpretability and accuracy [1]. KANs have univariate learnable activation
function on edges and nodes are summing those activation functions. The
integration of KAN, ensembles of probabilistic trees, and multivariate
B-spline representations is presented in [13]. However, most of the previous
works on KAN is for depth-2 representation; except, Liu et. al. [1] which
extends KANs to to arbitrary widths and depths. As [1] is based on B-Spline,
we call it Spl-KAN.
In this paper, we present an improved version of KAN, called Wav-KAN which
uses wavelet in KAN configuration. Figure 1 shows a Wav-KAN 444These fancy
plots are just to show the potential of wavelet for approximation; they are
not output of activation functions of a real network; hence, we see some of
them are not differentiable. with 2 input features, 3 hidden nodes, and 2
output nodes (Wav-KAN[2,3,2]). In general, it can have an arbitrary number of
layers. The structure is similar to MLPs with weights replaced by wavelet
functions, and nodes are doing summation of those Wavelet functions.

Figure 1: Wav-KAN with arbitrary number of layers (here is Wav-KAN[2,3,2])

Wavelet has been extensively employed in multiresolution analysis [14]. There
are some studies of wavelet in neural networks in Universal Approximation
Therory (UAT). For example, the complex Gabor wavelet activation function is
renowned for its optimal concentration in the space-frequency domain and its
exceptional ability to represent images accurately. Hence, [15] used Gabor
wavelet in implicit neural representation (INR) in application including image
denoising, image inpainting, super-resolution, computed tomography
reconstruction, image overfitting, and novel view synthesis with neural
radiance fields. To the best of our knowledge, our proposed framework is the
first work which uses wavelet in Kolmogorov-Arnold representation theorem for
an arbitrary widths and depth neural netwroks.
In comparison with the recent proposed Spl-KAN [1] and also MLPs, the proposed
configuration is faster, more accurate, and more robust, effectively
addressing the existing issues and significantly enhancing its performance.
Also, in statistical learning, we look for models which are more flexible and
more interpretable [16]. Wav-KAN will be a powerful member of a family of KAN
glass-box models which facilitates the interpretability of neural networks.
This work introduces Wav-KAN and we believe with the potential of Wav-KAN
which combines the potential of wavelet and KAN, it will be widely used in all
the fields and it will be implemented in Pytorch, Tensor-flow, R, to name but
few.

Incorporating Wav-KAN in neural networks makes increasingly explainable
components that can achieve similar state-of-the-art performance across
diverse tasks 555Some of the proofs and experiments will be presented in the
final version of the paper..

Section II discusses KAN and its generalization to multi-layers KAN with
wavelet as function approximation (i.e., activation function). Section III
presents Continuous Wavelet Transform, especially the definition and criteria
for being a mother wavelet which can be used as basis for function
approximation. The comparison of Wav-KAN, Spl-KAN and MLPs is given in V. Some
experiment results will be given in Section VI. Finally, Section VII concludes
the paper.

##  II Kolmogorov-Arnold Networks

Kolmogorov-Arnold Networks (KANs) represent a novel twist on neural network
design that challenges traditional concepts like the Multi-Layer Perceptron
(MLP). At the heart of KANs is a beautiful and somewhat abstract mathematical
theorem from Kolmogorov and Arnold.

###  II-A The Kolmogorov-Arnold Representation Theorem

Let’s start with the theorem that inspires KANs:

Theorem: Kolmogorov-Arnold Representation [12] For any continuous function $f$
of $n$ variables defined on a cube $[0,1]^{n}$, there exist $2n+1$ functions
$\phi_{q}$ and $2n+1\times n$ functions $\psi_{q,p}$, all univariate and
continuous, such that

$f(x_{1},\ldots,x_{n})=\sum_{q=1}^{2n+1}\Phi_{q}\left(\sum_{p=1}^{n}\phi_{q,p}(x_{p})\right).$
(1)

This theorem tells us that any multivariate function can essentially be
decomposed into the sum of functions of sums. The magic here is that these
inner functions are univariate, meaning they each take a single input. As we
use mother wavelet as our basis, in Wav-KAn notation, we use
$\psi_{i,j}(x_{j})$ instead of $\phi_{i,j}(x_{j})$ and $\Psi_{i}$ instead of
$\Phi_{i}$.

###  II-B From Theory to Networks

How do we translate this theorem into a neural network architecture? Imagine
this: Instead of weights and biases adjusting linear combinations of inputs at
each node, KANs modify this process to work with functions. One generalized
version of KAN theorem that corresponds to deeper KANs is the recently
published work of Spl-KAN[1].

  * •

In KANs, every ”weight” is actually a small function on its own. Each node in
a KAN does not apply a fixed non-linear activation function. Each learnable
activation function in edges, gets the input and gives an output.

Suppose we have an MLPs neural network with $n$ input and $m$ output in a
fully connected layer (between layer $l$ and $l+1$). The equation in matrix
form is given by:

$\mathbf{x}^{(l+1)}=\mathbf{W_{\mathit{l+1,l}}}\mathbf{x}^{(l)}+\mathbf{b}^{(l+1)}$
(2)

where:

  * •

$\mathbf{W_{\mathit{l+1,l}}}$ is the weight matrix connecting layer $l$ and
layer $l+1$,

  * •

$\mathbf{x}^{(l)}$ is the input vector,

  * •

$\mathbf{x}^{(l+1)}$ is the output vector,

  * •

$\mathbf{b}^{(l+1)}$ is the bias vector for layer $l+1$.

The weight matrix $\mathbf{W_{\mathit{l+1,l}}}$ is expanded as follows:

$\mathbf{W_{\mathit{l+1,l}}}=\begin{pmatrix}w_{1,1}&w_{1,2}&\cdots&w_{1,n}\\\
w_{2,1}&w_{2,2}&\cdots&w_{2,n}\\\ \vdots&\vdots&\ddots&\vdots\\\
w_{m,1}&w_{m,2}&\cdots&w_{m,n}\end{pmatrix}$ (3)

where $w_{i,j}$, $i=1,2,...,m$ and $j=1,2,...,n$ is the weight between $i$-th
node in $l+1$-th layer, and $j$-th node in $l$-th .The bias vector
$\mathbf{b}^{(l+1)}$ is:

$\mathbf{b}^{(l+1)}=\begin{pmatrix}b_{1}^{(l+1)}\\\ b_{2}^{(l+1)}\\\ \vdots\\\
b_{m}^{(l+1)}\end{pmatrix}$ (4)

Thus, the complete equation becomes:

$\begin{pmatrix}x_{1}^{(l+1)}\\\ x_{2}^{(l+1)}\\\ \vdots\\\
x_{m}^{(l+1)}\end{pmatrix}=\begin{pmatrix}w_{1,1}&w_{1,2}&\cdots&w_{1,n}\\\
w_{2,1}&w_{2,2}&\cdots&w_{2,n}\\\ \vdots&\vdots&\ddots&\vdots\\\
w_{m,1}&w_{m,2}&\cdots&w_{m,n}\end{pmatrix}\begin{pmatrix}x_{1}^{(l)}\\\
x_{2}^{(l)}\\\ \vdots\\\
x_{n}^{(l)}\end{pmatrix}+\begin{pmatrix}b_{1}^{(l+1)}\\\ b_{2}^{(l+1)}\\\
\vdots\\\ b_{m}^{(l+1)}\end{pmatrix}$ (5)

Now, suppose we have $L$ layers, each of them having the structure described
above. Let $\sigma(\cdot)$ be the activation function. The compact formula for
the whole network, $f(\mathbf{x})$, where $\mathbf{x}$ is the input vector and
$f(\cdot)$ is the neural network, is given by:

$f(\mathbf{x})=\mathbf{x}^{(L)}$ (6)

where

$\mathbf{x}^{(l+1)}=\sigma(\mathbf{W_{\mathit{l+1,l}}}\mathbf{x}^{(l)}+\mathbf{b}^{(l+1)})$
(7)

for $l=0,1,2,\ldots,L-1$, and $\mathbf{x}^{(0)}=\mathbf{x}$ is the input
vector.

$f(\mathbf{x})=\sigma\left(\mathbf{W_{L}}\sigma\left(\mathbf{W_{L-1}}\cdots\sigma\left(\mathbf{W_{2}}\sigma\left(\mathbf{W_{1}}\mathbf{x}+\mathbf{b_{1}}\right)+\mathbf{b_{2}}\right)\cdots+\mathbf{b_{L-1}}\right)+\mathbf{b_{L}}\right)$
(8)

In KAN, the relationship between layers turns into: Let $\mathbf{x}^{(l)}$ be
a vector of size $n$. We transpose $\mathbf{x}^{(l)}$ and place it in a matrix
$\mathbf{X}$ with $m$ rows and $n$ columns:

$\mathbf{x}^{(l)}\in\mathbb{R}^{n}\quad\Rightarrow\quad(\mathbf{x}^{(l)})^{T}\in\mathbb{R}^{1\times
n}$

We construct the matrix $\mathbf{X}_{l}$ as follows:

$\mathbf{X}_{l}=\begin{pmatrix}(\mathbf{x}^{(l)})^{T}\\\
(\mathbf{x}^{(l)})^{T}\\\ \vdots\\\
(\mathbf{x}^{(l)})^{T}\end{pmatrix}\in\mathbb{R}^{m\times n}$

where each row of $\mathbf{X}_{l}$ is the transposed vector
$(\mathbf{x}^{(l)})^{T}$.

We define the operator $T_{o}$ which acts on the matrix
$\Psi_{l+1,l}(\mathbf{X}_{l})$. This operator sums the elements of each row of
the matrix and outputs the resulting vector $\mathbf{v}$. The definition is as
follows:

$T_{o}\left(\Psi_{l+1,l}(\mathbf{X}_{l})\right)=\mathbf{v}$

where $\mathbf{v}$ is a vector with elements given by:

$v_{i}=\sum_{j}\left[\Psi_{l+1,l}(\mathbf{X}_{l})\right]_{ij}=\sum_{j=1}^{n}\psi_{i,j}(x_{j}^{(l)}),\quad\text{for
}i=1,2,\ldots,m$

In this expression, $\left[\Psi_{l+1,l}(\mathbf{X}_{l})\right]_{ij}$
represents the element in the $i$-th row and $j$-th column of the matrix
$\Psi_{l+1,l}(\mathbf{x}^{(l)})$.

Thus, the operator $T_{o}$ can be written as:

$T_{o}\left(\Psi_{l+1,l}(\mathbf{X}_{l})\right)=\left(\sum_{j}\left[\Psi_{l+1,l}(\mathbf{X}_{l})\right]_{ij}\right)_{i}$

In this definition, $T_{o}$ takes the matrix $\Psi_{l+1,l}(\mathbf{X}_{l})$,
sums the elements of each row, and outputs the resulting vector $\mathbf{v}$.

Indeed $\Psi_{l+1,l}$ acts on the input vector $\mathbf{x}^{(l)}$ and gives an
output where each element of $\Psi_{l+1,l}$ takes one corresponding element of
$\mathbf{x}^{(l)}$ as input, sums up the results, and produces one element of
the output:

$\mathbf{X}_{l+1}=\Psi_{l+1,l}(\mathbf{X}_{l})$ (9)

where:

$\Psi_{l+1,l}(\mathbf{X}_{l})=\begin{pmatrix}\psi_{1,1}(x_{1}^{(l)})&\psi_{1,2}(x_{2}^{(l)})&\cdots&\psi_{1,n}(x_{n}^{(l)})\\\
\psi_{2,1}(x_{1}^{(l)})&\psi_{2,2}(x_{2}^{(l)})&\cdots&\psi_{2,n}(x_{n}^{(l)})\\\
\vdots&\vdots&\ddots&\vdots\\\
\psi_{m,1}(x_{1}^{(l)})&\psi_{m,2}(x_{2}^{(l)})&\cdots&\psi_{m,n}(x_{n}^{(l)})\end{pmatrix}$
(10)

Here, $\Psi_{l+1,l}$ represents the activation functions connecting layer $l$
and layer $l+1$. Each element $\psi_{i,j}(\cdot)$ denotes the activation
function that connects the $j$-th neuron in layer $l$ to the $i$-th neuron in
layer $l+1$. Instead of multiplication, equation (9) computes a function of
input with distinct learnable parameters.

Hence, if $\mathbf{X}_{0}$ be considered as input which just contains input
vector as its rows, for the entire network, the output after $L$ layers is:

$f_{KAN}(\mathbf{X}_{0})=\mathbf{x}^{(L)}=T_{o}\left(\Psi_{L,L-1}\left(\mathbf{X}_{L-1}\right)\right)=T_{o}\left(\Psi_{L,L-1}\left(\begin{pmatrix}(T_{o}(\Psi_{L-1,L-2}(\mathbf{X}_{L-2})))^{T}\\\
(T_{o}(\Psi_{L-1,L-2}(\mathbf{X}_{L-2})))^{T}\\\ \vdots\\\
(T_{o}(\Psi_{L-1,L-2}(\mathbf{X}_{L-2})))^{T}\end{pmatrix}\right)\right)$

$=\cdots=T_{o}\left(\Psi_{L,L-1}\left(\begin{pmatrix}(T_{o}(\Psi_{L-1,L-2}\cdots(T_{o}(\Psi_{1,0}(\mathbf{X}_{0}))))^{T}\\\
(T_{o}(\Phi_{L-1,L-2}\cdots(T_{o}(\Psi_{1,0}(\mathbf{X}_{0}))))^{T}\\\
\vdots\\\
(T_{o}(\Psi_{L-1,L-2}\cdots(T_{o}(\Psi_{1,0}(\mathbf{X}_{0}))))^{T}\end{pmatrix}\right)\right)$
(11)

In summary, traditional MLPs use fixed nonlinear activation functions at each
node and linear weights (and biases) to transform inputs through layers. The
output at each layer is computed by a linear transformation followed by a
fixed activation function. During backpropagation, gradients of the loss
function with respect to weights and biases are calculated to update the model
parameters. In contrast, KANs replace linear weights with learnable univariate
functions placed on edges rather than nodes. In the nodes, we just have a
summation of some univariate functions from previous layers. Each function is
adaptable, allowing the network to learn both the activation and
transformation of the inputs. This change leads to improved accuracy and
interpretability, as KANs can better approximate functions with fewer
parameters. During backpropagation in KANs, the gradients are computed with
respect to these univariate functions, updating them to minimize the loss
function. This results in more efficient learning for complex and high-
dimensional functions.

###  II-C Why Bother with KANs?

The flexibility of KANs allows for a more nuanced understanding and adaptation
to data. By learning the functions directly involved in data relationships,
KANs aim to provide a more accurate and interpretable model:

  * •

Accuracy: They can fit complex patterns in data more precisely with
potentially fewer parameters.

  * •

Interpretability: Since each function has a specific, understandable role,
it’s easier to see what the model is ”thinking.”

In summary, KANs leverage deep mathematical insights to offer a fresh
perspective on how neural networks can understand and interact with the world.
By focusing on functions rather than mere weights, they promise a richer and
more intuitive form of machine learning.

##  III Continuous Wavelet Transform

The Continuous Wavelet Transform (CWT) is a method mostly used in signal
processing to analyze the frequency content of a signal as it varies over time
[14]. It acts like a microscope, zooming in on different parts of a signal to
determine its constituent frequencies and their variations.

CWT utilizes a base function known as a “mother wavelet,” which serves as a
template that can be scaled and shifted to match various parts of the signal.
The shape of the mother wavelet is critical as it dictates which features of
the signal are highlighted.

Let $\psi\in L^{2}(\mathbb{R})$ be the mother wavelet and $g(t)\in
L^{2}(\mathbb{R})$666
$L^{2}(\mathbb{R})=\left\\{f(x)\mid\int|f(x)|^{2}\,dx<\infty\right\\}$ . be
the function that we want to express in the wavelet basis. Then, a mother
wavelet must satisfy certain criteria [17, 18].

  1. 1.

Zero Mean: The integral of the wavelet over its entire range must equal zero:

$\int_{-\infty}^{\infty}\psi(t)\,dt=0$ (12)

  2. 2.

Admissibility Condition: The wavelet must have finite energy, which means the
integral of the square of the wavelet must be finite.

$C_{\psi}=\int_{0}^{+\infty}\frac{|\hat{\psi}(\omega)|^{2}}{\omega}\,d\omega<+\infty$
(13)

where $\hat{\psi}$ is the Fourier transform of the wavelet $\psi(t)$.

The CWT of a signal/function is represented by wavelet coefficients,
calculated as follows:

$C(s,\tau)=\int_{-\infty}^{+\infty}g(t)\frac{1}{\sqrt{s}}\psi\left(\frac{t-\tau}{s}\right)dt$
(14)

where:

  * •

$g(t)$ is the signal/function that we want to approximate by Wavelet basis.

  * •

$\psi(t)$ is the mother wavelet.

  * •

$s\in\mathbb{R}^{+}$ is the scale factor which is greater than zero.

  * •

$\tau\in\mathbb{R}$ is the shift factor.

  * •

$C(s,\tau)$ measures the match between the wavelet and the signal at scale $s$
and shift $\tau$.

A signal can be reconstructed from its wavelet coefficients using the inverse
CWT:

$g(t)=\frac{1}{C_{\psi}}\int_{-\infty}^{+\infty}\int_{0}^{+\infty}C(s,\tau)\frac{1}{\sqrt{s}}\psi\left(\frac{t-\tau}{s}\right)\frac{ds\,d\tau}{s^{2}}$
(15)

where $C_{\psi}$ is a constant that depends on the wavelet, ensuring the
reconstruction’s accuracy.

##  IV Discrete Wavelet Transform

The Discrete Wavelet Transform (DWT) is a widely used method in signal
processing that decomposes a signal into different frequency components with
high efficiency [14]. Unlike the Continuous Wavelet Transform (CWT), which
analyzes the signal at every possible scale and translation, the DWT uses
discrete sampling intervals which is suitable for digital applications. When
the data points are close together (high sampling rate), the wavelet basis can
zoom in to capture the high-frequency details (fine details) in that local
area. When the data points are spread out (low sampling rate), the wavelet
basis can zoom out to capture the low-frequency trends (overall shape) using
global information. We see some examples uniform and irregular samples of
data. Indeed, by employing wavelet we combine local detailed information where
the data points are dense with broader trends where the data points are
sparse.

The DWT employs a set of base functions derived from a single mother wavelet
through scaling and translation. These base functions are orthonormal,
ensuring that the DWT provides a compact and non-redundant representation of
the signal. By applying high-pass and low-pass filters iteratively to the
signal, the DWT produces approximation and detail coefficients at various
levels of resolution, allowing for a multi-resolution analysis.

Let $\psi\in L^{2}(\mathbb{R})$ be the mother wavelet and $\phi\in
L^{2}(\mathbb{R})$ be the scaling function. The signal $g(t)\in
L^{2}(\mathbb{R})$777
$L^{2}(\mathbb{R})=\left\\{f(x)\mid\int|f(x)|^{2}\,dx<\infty\right\\}$ . is
decomposed into approximation coefficients $a_{j}(k)$ and detail coefficients
$d_{j}(k)$ at level $j$ as follows:

$a_{j}(k)=\sum_{n}g(n)\phi_{j,k}(n)$ (16)

$d_{j}(k)=\sum_{n}g(n)\psi_{j,k}(n)$ (17)

where:

  * •

$g(n)$ is the discrete signal.

  * •

$\phi_{j,k}(n)$ is the scaling function at scale $j$ and position $k$.

  * •

$\psi_{j,k}(n)$ is the wavelet function at scale $j$ and position $k$.

  * •

$a_{j}(k)$ represents the approximation coefficients.

  * •

$d_{j}(k)$ represents the detail coefficients.

The original signal can be reconstructed from its wavelet coefficients using
the inverse DWT, ensuring that no information is lost during the
transformation process. This reconstruction is given by:

$g(n)=\sum_{k}a_{j}(k)\phi_{j,k}(n)+\sum_{j=1}^{J}\sum_{k}d_{j}(k)\psi_{j,k}(n)$
(18)

where $J$ is the number of decomposition levels.

The DWT’s ability to provide both time and frequency localization, coupled
with its computational efficiency, makes it an indispensable tool in various
fields such as image compression, noise reduction, and feature extraction.

Multi-resolution analysis (MRA) using wavelets, specifically discrete wavelet
transforms (DWT), is a powerful technique in signal processing and data
analysis. MRA allows the decomposition of a signal into different levels of
detail by using wavelets, which are localized functions. The DWT achieves this
by recursively applying high-pass and low-pass filters to the data/signal,
creating approximations and details at various resolutions. This process
provides a hierarchical framework, where each level captures different
frequency components of the signal, making it easier to analyze localized
features and transient phenomena. The ability to zoom in on fine details while
retaining the broader structure of the signal is a key advantage of MRA,
making it useful in applications such as image compression, noise reduction,
and feature extraction.

In the context of discrete wavelet transforms, the data/signal is represented
as a sum of wavelet functions at different scales and positions. In wavelet
analysis, once we compute the initial set of coefficients, we do not need to
recompute them when calculating additional coefficients at finer resolutions.
This efficiency arises because wavelet transforms leverage previously computed
coefficients to generate further details. Consequently, the hierarchical
nature of the multi-resolution analysis ensures that the earlier computations
remain valid and reusable, streamlining the process and significantly reducing
the computational overhead. This property makes wavelet-based methods
particularly advantageous for iterative and real-time signal processing
applications. The multi-resolution property of wavelets is particularly
beneficial for analyzing non-stationary signals, which have frequency
components that vary over time. By breaking down the signal into different
resolution levels, MRA enables a more comprehensive understanding of its
structure and characteristics. Moreover, the DWT’s ability to provide both
time and frequency localization offers a significant advantage over
traditional Fourier transforms, which only provide frequency information 888We
will add some simulations to this section..

In summary, wav-kan by employing wavelet adapts to the varying density of data
points, using more detailed information where there are more data points and
less detailed information where there are fewer data points. This approach is
particularly useful for real-world data that often suffers from irregular
sampling and data dropouts.

##  V Wav-KAN or Spl-KAN or MLPs?

Wavelets and B-splines are two prominent methods used for function
approximation, each with distinct advantages and limitations, particularly
when applied in neural networks. B-splines provide smooth and flexible
function approximations through piecewise polynomial functions defined over
control points. They offer local control, meaning adjustments to a control
point affect only a specific region, which is advantageous for precise
function tuning. This smoothness and local adaptability make B-splines
suitable for applications requiring continuous and refined approximations,
such as in CAD and computer graphics. However, the computational complexity
increases significantly with higher dimensions, making them less practical for
high-dimensional data. Managing knot placement can also be intricate and
affect the overall shape and accuracy of the approximation. While B-splines
can be used in neural networks to approximate activation functions or smooth
decision boundaries, their application is generally less suited to feature
extraction tasks compared to wavelets due to their limited ability to handle
multi-resolution analysis and sparse representations.

Wavelets excel in multi-resolution analysis enabling different levels of
detail to be represented simultaneously which makes it a precious tool to
decompose data into various frequency components. This capability is highly
beneficial for feature extraction in neural networks, as it enables capturing
of both high-frequency details and low-frequency trends. Additionally,
wavelets offer sparse representations, which can lead to more efficient neural
network architectures and faster training times. They are also well-suited for
handling non-stationary signals and localized features, making them ideal for
applications such as image recognition and signal classification. However, the
choice of the wavelet function is crucial and can significantly impact
performance, and edge effects can introduce artifacts that need special
handling.

In function approximation, wavelets excel by maintaining a delicate balance
between accurately representing the underlying data structure and avoiding
overfitting to the noise. Unlike traditional methods that may overly smooth
the data or fit to noise, wavelets achieve this balance through their inherent
ability to capture both local and global features. By decomposing the data
into different frequency components, wavelets can isolate and retain
significant patterns while discarding irrelevant noise. This multiresolution
analysis ensures that the approximation is robust and reliable, providing a
more accurate and nuanced representation of the original data without the
pitfalls of noise overfitting. On the other side, while Spl-KAN is powerful in
capturing the changes in data, it also captures the noise in training data.
Indeed, the strength of Spl-KAN is its weakness, too.
The advantages of Spl-KAN which are mentioned in [1] including
interpretability and/or accuracy with respect to MLPs exist for Wav-KAN. More
importantly Wav-KAN has solved the major disadvantage of Spl-KANs which was
slow training speed. In terms of number of parameter, we compared Wav-KAN with
Spl-KAN and MLPs for a hypothetical neural network which has $N$ input nodes
and $N$ output node, with $L$ layers. As we see in Table I and by considering
the value of G, Wav-KAN has less number of parameters than Spl-KAN (k should
be at least 2 to have good results in Spl-KAN, especially in complex tasks).
The coefficient 3 is becuase Wav-KAN has a learnable weight, a translation and
a scaling. Learnable parameters in each neural network is in column of
parameters. While the order of MLPs are less for the hypothetical neural
network, in practice, Wav-KAN needs less number of parameters to learn the
same task. Indeed, this originates because of capacity of Wavelet for
capturing both low frequency and high frequency functions.

TABLE I: Comparison of MLPs, Spl-KAN and Wav-KAN Neural network With $L$
layers, each layer has $N$ nodes
---
Neural Network Structure | Order | Parameters
MLPs |  $O(N^{2}L)$ or $O(N^{2}L+NL)$ | weights and biases
Spl-KAN | $O(N^{2}L(G+k+1))\sim O(N^{2}LG)$ | weights
Wav-KAN | $O(3N^{2}L)$ | weight, translation, scaling

Regarding the implementation, Spl-KAN requires a smooth function like $b(x)$
in the equation (2.10) [1] attempting to catch on some global features.
Because of the inherent scaling property of Wavelet, Wav-KAN does not need an
additional term in its activation functions. This helps Wav-KAN to be faster.
For example, wavelets equal to the second derivative of a Gaussian which are
called Mexican hats 999The following equation is normalized Mexican hat
wavelet. Indeed, pywavelets has a minus sign behind it which because of the
weight we put behind the mother wavelet in our activation functinos, the sign
doesn’t matter(https://pywavelets.readthedocs.io/en/latest/ref/cwt.html). and
first used in computer vision to detect multiscale edges [19, 14] has the
following form

$\psi(t)=\frac{2}{\pi^{1/4}\sqrt{3\sigma}}\left(\frac{t^{2}}{\sigma^{2}}-1\right)\exp\left(\frac{-t^{2}}{2\sigma^{2}}\right).$
(19)

Where $\sigma$ shows the adjustable standard deviation of Gaussian. In our
experiments, $\psi_{exp}(t)$

$\psi_{exp}(t)=w\psi(t)$ (20)

Indeed, $w$ plays the role of CWT coefficients which is multiplied by the
mother wavelet formula; as $w$ is a learnable parameter, it helps adapting the
shape of the mother wavelet to the function that it tries to approximate .

Moreover, Spl-KAN heavily depends on the grid spaces, and for better
performance, it requires increasing the number of grids 101010Which is
equivalent to decrease the spaces between the grids. This helps Spl-KAN
accuracy to be improved., though, it brings two disadvantages. First it needs
curvefitting which is a cumbersome and computational expensive operation, and
also while we increase the number of grids, loss has some jumps 111111Loss
will be increased for some training, then, it works better. See Figure 2.3 [1]
for better illustration.. Fortunately, wavelet is safe from such computations
and deficiencies and if one wants to capture more details, DWT efficiently
does that without recalculation of the previous steps.

Last but no least, we found that batch normalization [20] significantly
improves accuracy and speeds up training of both Wav-KAN and Spl-KAN; hence,
we included batch normalization in both of these methods 121212In [1], the
authors did not mention and didn’t apply batch normalization..

##  VI Simulation Results

In this section, we present the results of our experiments conducted using the
KAN (Kernel-based Artificial Neural network) model with various continuous
wavelet transformations (CWTs) 131313We add more simulations including some
discrete wavelet ones for MRA, etc. on the MNIST dataset, utilizing a training
set of 60,000 images and a test set of 10,000 images. It is important to note
that our objective was not to optimize the parameters to their best possible
values, but rather to demonstrate that Wav-KAN performs well in terms of
overall performance. We have incorporated batch normalization into both Spl-
KAN and Wav-KAN, resulting in improved performance. The wavelet types
considered in our study include Mexican hat, Morlet, Derivative of Gaussian
(DOG), and Shannon (see Table II). For each wavelet type and also Spl-KAN, we
performed five trials, training the model for 50 epochs per trial.

TABLE II: Mother Wavelet Formulas and Parameters Wavelet Type | Formula of Mother Wavelet | Parameters
---|---|---
Mexican hat | $\psi(t)=\frac{2}{\sqrt{3}\pi^{1/4}}\left(t^{2}-1\right)e^{-\frac{t^{2}}{2}}$ |  $\tau$, $s$
Morlet | $\psi(t)=\cos(\omega_{0}t)e^{-\frac{t^{2}}{2}}$ |  $\tau$, $s$ and scaling, $\omega_{0}$= 5
Derivative of Gaussian (DOG) | $\psi(t)=-\frac{d}{dt}\left(e^{-\frac{t^{2}}{2}}\right)$ |  $\tau$, $s$
Shannon | $\psi(t)=\text{sinc}(t/\pi)\cdot w(t)$ |  $\tau$, $s$, and $w(t)$: window function

The results were averaged across these trials to ensure the robustness and
reliability of our findings.
Both Wav-KAN and Spl-KAN have the structure (number of nodes) of [28*28,32,10]
141414[first layer nodes, middle layer nodes, output nodes]. Although we
enhanced Spl-KAN by using spline order 3 and a grid size of 5, this approach
is computationally much more expensive compared to Wav-KAN. We employed AdamW
optimizer [21, 22], with learning rate of 0.001 with weight decay of
$10^{-4}$. Loss is cross entropy.

Figure 2:  Training accuracy of Wav-KAN [28*28,32,10] versus Spl-KAN
[28*28,32,10] Figure 3:  Test accuracy of Wav-KAN [28*28,32,10] versus Spl-KAN
[28*28,32,10]

Figures 2 and 3 show the result of training accuracy and test accuracy of Wav-
KAN in comparison to Spl-KAN. To not clutter up, we just show the result of
Derivative of Gaussian (DOG) and Mexican hat wavelet. These have been shown as
a sample. By fine tuning, and using all the freedom of wavelet (like frequency
of sinusoid and variance of Gaussian), Wavelet shows significant superiority
151515We have done a lot of such experiments and we will publish them soon..
While Spl-KAN has better performance in training, which is because of
overfitting to data, a lot of wavelet types have shown superior performance
with respect to Spl-KAN. Indeed, for this experiment we set the variance of
the Gaussian in wavelets to be 1; though, we can find better variances by grid
search or making it a learnable parameter. Also, one may like to consider
$\omega_{0}$ in Morlet and parameters of frequency and/or windowing in Shannon
wavelet learnable. Indeed, we can use different degrees of freedom that we
have to make the shape of wavelets more flexible. Wavelet makes a balance by
not fitting to the noise in the data.

We evaluated the performance of each wavelet type in terms of training loss,
training accuracy, validation loss, and validation accuracy. Figures 1 and 2
summarize the results, depicting the training and validation metrics averaged
over the five trials for each wavelet type.
The simulation results indicate that the choice of wavelet significantly
impacts the performance of the KAN model. This suggests that these wavelets
are particularly effective at capturing the essential features of the MNIST
dataset while maintaining robustness against noise. On the other hand,
wavelets like Shannon and Bump did not perform as well, highlighting the
importance of wavelet selection in designing neural networks with wavelet
transformations.

##  VII Conclusion

In this paper, we have introduced Wav-KAN, a novel neural network architecture
that integrates wavelet functions within the Kolmogorov-Arnold Networks (KAN)
framework to enhance both interpretability and performance. By leveraging the
multiresolution analysis capabilities of wavelets, Wav-KAN effectively
captures complex data patterns and provides a robust solution to the
limitations faced by traditional multilayer perceptrons (MLPs) and recently
proposed Spl-KANs.

Our experimental results demonstrate that Wav-KAN not only achieves superior
accuracy but also benefits from faster training speeds compared to Spl-KAN.
The unique structure of Wav-KAN, which combines the strengths of wavelet
transforms and the Kolmogorov-Arnold representation theorem, allows for more
efficient parameter usage and improved model interpretability.

Wav-KAN represents a significant advancement in the design of interpretable
neural networks. Its ability to handle high-dimensional data and provide clear
insights into model behavior makes it a promising tool for a wide range of
applications, from scientific research to industrial deployment. Future work
will focus on further optimizing the Wav-KAN architecture, exploring its
applicability to other datasets and tasks, and implementing the framework in
popular machine learning libraries such as PyTorch and TensorFlow.

Overall, Wav-KAN stands out as a powerful and versatile model, paving the way
for the development of more transparent and efficient neural network
architectures. Its potential to combine high performance with interpretability
marks a crucial step forward in the field of artificial intelligence.

## References

  * [1] Z. Liu, Y. Wang, S. Vaidya, F. Ruehle, J. Halverson, M. Soljačić, T. Y. Hou, and M. Tegmark, “Kan: Kolmogorov-arnold networks,” _arXiv preprint arXiv:2404.19756_ , 2024.
  * [2] D. Hendrycks, M. Mazeika, and T. Woodside, “An overview of catastrophic ai risks,” _arXiv preprint arXiv:2306.12001_ , 2023.
  * [3] R. Ngo, L. Chan, and S. Mindermann, “The alignment problem from a deep learning perspective,” _arXiv preprint arXiv:2209.00626_ , 2022.
  * [4] Y. Zhang, P. Tiňo, A. Leonardis, and K. Tang, “A survey on neural network interpretability,” _IEEE Transactions on Emerging Topics in Computational Intelligence_ , vol. 5, no. 5, pp. 726–742, 2021.
  * [5] F. Doshi-Velez and B. Kim, “Towards a rigorous science of interpretable machine learning,” _arXiv preprint arXiv:1702.08608_ , 2017.
  * [6] A. Pinkus, “Approximation theory of the mlp model in neural networks,” _Acta numerica_ , vol. 8, pp. 143–195, 1999.
  * [7] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, “Attention is all you need,” _Advances in neural information processing systems_ , vol. 30, 2017.
  * [8] H. Cunningham, A. Ewart, L. Riggs, R. Huben, and L. Sharkey, “Sparse autoencoders find highly interpretable features in language models,” _arXiv preprint arXiv:2309.08600_ , 2023.
  * [9] F.-L. Fan, J. Xiong, M. Li, and G. Wang, “On interpretability of artificial neural networks: A survey,” _IEEE Transactions on Radiation and Plasma Medical Sciences_ , vol. 5, no. 6, pp. 741–760, 2021.
  * [10] C. Olah, N. Cammarata, L. Schubert, G. Goh, M. Petrov, and S. Carter, “Zoom in: An introduction to circuits,” _Distill_ , vol. 5, no. 3, pp. e00 024–001, 2020.
  * [11] N. Elhage, R. Lasenby, and C. Olah, “Privileged bases in the transformer residual stream, 2023,” _URL https://transformer-circuits. pub/2023/privilegedbasis/index. html. Accessed_ , pp. 08–07, 2023.
  * [12] A. N. Kolmogorov, _On the representation of continuous functions of several variables by superpositions of continuous functions of a smaller number of variables_. American Mathematical Society, 1961.
  * [13] D. Fakhoury, E. Fakhoury, and H. Speleers, “Exsplinet: An interpretable and expressive spline-based neural network,” _Neural Networks_ , vol. 152, pp. 332–346, 2022.
  * [14] S. Mallat, _A wavelet tour of signal processing_. Elsevier, 1999.
  * [15] V. Saragadam, D. LeJeune, J. Tan, G. Balakrishnan, A. Veeraraghavan, and R. G. Baraniuk, “Wire: Wavelet implicit neural representations,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2023, pp. 18 507–18 516.
  * [16] G. James, D. Witten, T. Hastie, R. Tibshirani, and J. Taylor, _An introduction to statistical learning: With applications in python_. Springer Nature, 2023.
  * [17] A. Calderón, “Intermediate spaces and interpolation, the complex method,” _Studia Mathematica_ , vol. 24, no. 2, pp. 113–190, 1964.
  * [18] A. Grossmann and J. Morlet, “Decomposition of hardy functions into square integrable wavelets of constant shape,” _SIAM journal on mathematical analysis_ , vol. 15, no. 4, pp. 723–736, 1984.
  * [19] A. P. Witkin, “Scale-space filtering,” in _Readings in computer vision_. Elsevier, 1987, pp. 329–332.
  * [20] S. Ioffe and C. Szegedy, “Batch normalization: Accelerating deep network training by reducing internal covariate shift,” in _International conference on machine learning_. pmlr, 2015, pp. 448–456.
  * [21] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.
  * [22] I. Loshchilov and F. Hutter, “Decoupled weight decay regularization,” _arXiv preprint arXiv:1711.05101_ , 2017.


# Moiré-like superlattice generated van Hove singularities in strained CuO2
double layer

A.O. Sboychakov  Institute for Theoretical and Applied Electrodynamics,
Russian Academy of Sciences, 125412 Moscow, Russia  K.I. Kugel  Institute for
Theoretical and Applied Electrodynamics, Russian Academy of Sciences, 125412
Moscow, Russia  Institute of Metal Physics, Ural Branch, Russian Academy of
Sciences, 620990 Ekaterinburg, Russia  A. Bianconi  Rome International Center
for Materials Science Superstripes RICMASS, via dei Sabelli 119A, 00185 Rome,
Italy  Institute of Crystallography, CNR, Via Salaria Km 29.300, 00015
Monterotondo Rome, Italy  National Research Nuclear University, MEPhI (Moscow
Engineering Physics Institute), Kashirskoe sh. 31, 115409 Moscow, Russia

###### Abstract

While it is known that the double-layer Bi2Sr2CaCu2O8+y (BSCCO) cuprate
superconductor exhibits a one-dimensional (1D) incommensurate superlattice
(IS), the effect of IS on the electronic structure remains elusive. Following
the recent shift of the interest from underdoped to optimum and overdoped
phase in BSCCO by increasing the hole doping $x$, controlled by the oxygen
interstitials concentration $y$, here we focus on the multiple splitting of
the density of states (DOS) peaks and emergence of higher order van Hove
singularities (VHS) due to the 1D incommensurate superlattice. It is known
that 1D incommensurate wave vector $\textbf{q}=\epsilon\textbf{b}$ (where b is
the reciprocal lattice vector of the orthorhombic lattice) is controlled by
the misfit strain between different atomic layers in the range $0.209$ –
$0.215$ in BSCCO and in the range $0.209$ – $0.25$ in Bi2Sr2Ca1-xYxCu2O8+y
(BSCYCO). This work reports the theoretical calculation of a complex pattern
of VHS due to the 1D incommensurate superlattice with large 1D quasi-
commensurate supercells with the wave vector $\epsilon=9/\eta$ in the range
$36>\eta>43$. The similarity of the complex VHS splitting and appearing of
higher order VHS in a mismatched CuO2 bilayer with VHS due to the moiré
lattice in strained twisted bilayer graphene is discussed. This makes
mismatched CuO2 bilayer quite promising for constructing quantum devices with
tuned physical characteristics.

##  I Introduction

Different van der Waals heterostructures Liu et al. (2016), such as bilayer
graphene Rozhkov et al. (2016) and transition metal dichalcogenide bilayers
Weston et al. (2020); Zhao et al. (2021) provide an actively progressing
testing ground for a variety of quantum phenomena, e.g., flat bands, electron
nematicity, etc. Indeed, if we put together two layers of van der Waals
materials, a lattice mismatch or twist angle lead to the formation of moiré
patterns with large supercells. This provides a promising means of the band
and crystal structure engineering leading, in particular, to the moiré
modulated topological order Tong et al. (2017), and paving way to such novel
field of research as twistronics Ren et al. (2020).

Bilayers of atomic cuprate oxide superconductors appear to be not less
interesting. The lattice parameter and aperiodic local lattice distortions
associated with variable microstrain are not key physical parameters either in
BCS and unconventional superconductivity theories. Nevertheless, a compelling
evidence for the key role of CuO2 lattice in the mechanism of high-temperature
superconductivity and charge density wave (CDW) phenomena in cuprate
perovskites has been reported Krockenberger et al. (2021); Park et al. (2021);
Conradson et al. (2020); Sederholm et al. (2021); Jang et al. (2017). These
findings confirm early experimental results provided by local and fast
experimental methods based on the use of synchrotron X-ray radiation Bianconi
et al. (1996a); Saini et al. (2003); Bianconi et al. (1996b, c); Di Castro et
al. (2000) pointing toward the relevant role of the micro-strain in
superconducting cuprate perovskites Bianconi et al. (2000); Agrestini et al.
(2003).

The Bi2Sr2CaCu2O8+y (BSCCO) crystal Bianconi et al. (1996b, c); Di Castro et
al. (2000) is a misfit layer compound Wiegers (1996) exhibiting an
incommensurate composite structure Wiegers (1996); Janssen et al. (2004) and
large atomic displacements from average positions in the [BiO], [SrO], [CuO2],
and [Ca] layers. These displacements form an incommensurate composite
structure, where a texture of nanoscale insulating and metallic atomic
stripes, see for example Etrillard et al. (2000); Grebille et al. (2001);
Bianconi (1994); Beskrovnyy and Jirák (2012) Bi2Sr2Ca1-xYxCu2O8+y (BSCYCO)
with six different elements in the average crystalline unit cell forming an
archetypal case of High Entropy Perovskites (HEPs). The search of new HEPs is
today a hot topic focusing on the manipulation of the complexity, where the
entropy control plays a key role in the determination of material
functionalities like it is observed in bismuth 2212 iron ferrite and ferrate
compounds Pérez et al. (1997); Sedykh et al. (2003), thermoelectric misfit
layered cobaltite oxide perovskites Lambert et al. (2001); Grebille et al.
(2007) also in the form of 2D flakes Kim et al. (2012) taking advantage of the
polymorphism of perovskite structure with competing nanoscale crystallographic
phases Gavrichkov et al. (2019); Mazza et al. (2022).

Spatially correlated incommensurate lattice stripy modulations in
Bi2Sr2CaCu2O8+y superconductor can now be visualized by the scanning submicron
X-ray diffraction focusing on the superlattice satellites of the main X-ray
reflections Kaneko et al. (2004); Poccia et al. (2011, 2020). The results show
that, while oxygen concentration and the doping level can change depending on
the sample treatments and thermal history, the incommensurate superlattice
wave vector is controlled by the misfit strain and it exhibits a remarkable
stability. The critical temperature exhibits a dome as a function of oxygen
concentration or by potassium deposition Zhang et al. (2016), but the maximum
critical temperature of the dome Bianconi and Poccia (2012) is controlled by
the lattice stripe aperiodic structure Bianconi et al. (1993) and misfit
strain Bianconi et al. (2000).

Recently the scientific interest turned toward the normal state of cuprate
perovskites exhibiting high-temperature superconductivity in the overdoped
regime in the proximity of the superconductor-to-metal transition tuned by
doping, while for many years the focus was on the underdoped regime. The
investigations of the overdoped phase Alldredge et al. (2008); Pasupathy et
al. (2008); Gomes et al. (2007); Božović et al. (2016); He et al. (2021);
Drozdov et al. (2018); Li et al. (2018); Maier et al. (2020) of Bi2212 have
found evidence that the superconducting gap appearing at temperatures well
above the critical temperature $T_{c}$, while below $T_{c}$ a large fraction
of the normal state remains uncondensed. Scanning tunneling microscopy has
shown that on the edge of the overdoped regime, the gap remains in isolated
puddles above $T_{c}$, so that superconducting puddles are intercalated by
normal metal phase. Moreover, at the same time, the Fermi surface topology
changes giving rise to an electronic topological Lifshitz transition. This
provides compelling evidence for a nanoscale phase separation giving a
granular landscape, where superconducting puddles are embedded in a normal
phase. The phase separation scenario was predicted to appear at the Lifshitz
transition several years ago M.Yu. Kagan et al. (2021); Rakhmanov et al.
(2017); Kugel et al. (2008a); Bianconi et al. (2015); Kugel et al. (2008b).

It is known that the interlayer interaction between CuO2 of BSCCO planes
produces a bilayer splitting of the 2D Cu(3$d$)O(2$p$) single band into two
bands with even and odd symmetry exhibiting a double-peaked logarithm
divergent density of states (DOS) due to standard van Hove singularities (VHS)
in a 2D metal. However, it is not known how the 1D superstructure controls the
VHS.

While many electronic band structure calculations of the average structure
BSCCO have been published in these last 33 years (see e.g., the papers
beginning from Ref. Massidda et al., 1988 up to one of the most recent papers,
Ref. Nokelainen et al., 2020), to our knowledge, the electronic band structure
calculations of BSCCO with the incommensurate superstructure with the wave
vector $0.209<q<0.215$ (see e.g., Ref. Hewat et al., 1989) known since 1989,
has never been reported. Here, we report the band structure calculation of the
strained CuO2 bilayer with the incommensurate structure focusing on the high
doping range, where correlations become negligible and the van Hove
singularity crosses the chemical potential. We focus on the variation of the
electronic topological transition from a hole Fermi surface to an electron
Fermi surface in the overdoped regime in the BSCCO double layer perovskite
taking into account the aperiodic incommensurate lattice modulation.

In Bi2Sr2CaCu2O8+y, the mobility of oxygen interstitials in the CuO2 planes is
high in the temperature range 200–380 K, and it is possible to control the
concentration $y$ of oxygen interstitials to change the hole doping $x$ in the
layered perovskite.

The insulating Bi2Sr2YCu2O8 exhibits a commensurate modulation with
$\varepsilon=0.25$, which indicates the formation of lattice stripes with
commensurate period $n=1/\varepsilon=4$ lattice units. The superconducting
Bi2Sr2CaCu2O8+y exhibits a commensurate modulation with $\varepsilon=0.209$,
which indicate the formation of lattice stripes superlattice with
incommensurate period $n=1/0.209=9/43$ lattice units. For $\varepsilon=0.2$,
the CuO2 lattice should be decorated by a commensurate stripes superlattice
with $n=1/\varepsilon=1/0.2=5$ lattice units. For intermediate values of the
superlattice wave vector $0.25<\varepsilon<0.2$, the incommensurate
superstructure Bak (1982) consists of the mixture of the $\varepsilon=1/n$ and
$\varepsilon=1/m$ order, i.e. alternating the so called $n=4$ stripes portions
and $m=5$ stripes portions. Indeed, quasi-commensurate periods
$\lambda=1/\varepsilon=(nx+my)/(n+m)$ intermediate between two main
commensurate wave vectors $1/n$ and $1/m$ with periods of $m$ and $n$ lattice
units, respectively, occur at integer numbers of lattice unit cells
$\eta=(n+m)/\varepsilon$. Each quasi-commensurate phase (QCP) corresponds to a
modulation wave locked with the underlying lattice onto a rational number. The
sequence of quasi-commensurate phases is called the devil’s staircase Bak
(1982), which have been observed in many complex materials.

The quasi-commensurate modulation period in BSCCO is given by
$1/\varepsilon=(4x+5y)/(4+5)$, where $x+y=9$. Therefore,
$\eta=9/\varepsilon=4x+5y$ varies in the range $36<\eta<45$ and at each
integer value of $\eta$ the system reaches a quasi-commensurate phase at each
step of a devil’s staircase. We have found that in the BSCCO samples, where
the chemical potential is tuned by the concentration of oxygen interstitials,
the mismatched CuO2 bilayer is tuned by the misfit strain at the quasi-
commensurate modulation with $\eta=43$ i.e. $\epsilon$=$9/43=0.209$ in
agreement with the recent scanning micro-X-ray diffraction Poccia et al.
(2020), which is the devil’s staircase for $x=2$ and $y=5$. Therefore, we have
$4\times 2+5\times 7=43$ (i.e. a quasicrystal made of 8 insulating lattice
units and 35 metallic lattice units).

Figure 1: Schematics of the CuO2 layer 1 (left panel) and 2 (middle panel).
Layer 1 is unstretched. Layer 2 is either stretched or compressed in diagonal
($x$) direction. As a result, for some commensurate stretching/compression, a
superstructure appears. The superlattice contains $2N_{1}$ elementary unit
cells of the layer 1 and $2N_{2}$ elementary unit cells of the layer 2. Right
figure shows the sketch of the reduced Brillouin zone of the superlattice. The
$\Gamma$ point is in the coordinate origin, the $M$ point has coordinates
$\mathbf{M}=(0,\,\pi/d)$, while coordinate of the $X$ point is
$\mathbf{X}=(\pi/(N_{1}d),\,\pi/d)$. The thick blue triangle shows the
contour, along which the spectra (see below) are calculated.

##  II Geometry of the system

We model the system under study as consisting of two mismatched CuO2 layers.
The layer 1 is assumed to be quadratic. For the reason described below, we
choose the reference frame as shown on the left panel in Fig. 1. For such a
choice, the vectors connecting next nearest-neighbor Cu ions have coordinates
$\mathbf{d}_{1}=d(1,\,0)$ and $\mathbf{d}_{2}=d(0,\,1)$, where $d$ is the
distance between nearest copper ions in the diagonal direction. We choose
Bianconi et al. (1996b) $d=5.37$ Å. We also introduce vectors connecting
nearest-neighbor Cu ions of the layer 1 (lattice vectors of the unit cell of
layer 1). They are expressed via vectors $\mathbf{d}_{1,2}$ according to
$\mathbf{a}_{1}=(\mathbf{d}_{1}-\mathbf{d}_{2})/2$ and
$\mathbf{a}_{2}=(\mathbf{d}_{1}+\mathbf{d}_{2})/2$. The positions of Cu ions
in the layer 1 are

$\mathbf{r}^{d}_{1\mathbf{n}}\equiv\mathbf{r}_{1\mathbf{n}}=n\mathbf{a}_{1}+m\mathbf{a}_{2}\,,$
(1)

where $\mathbf{n}=(n,\,m)$ ($n$, $m$ are integers), while the positions of
oxygen ions are

$\mathbf{r}^{p_{x}}_{1\mathbf{n}}=\mathbf{r}_{1\mathbf{n}}+\frac{1}{2}\mathbf{a}_{1}\,,\;\;\mathbf{r}^{p_{y}}_{1\mathbf{n}}=\mathbf{r}_{1\mathbf{n}}+\frac{1}{2}\mathbf{a}_{2}\,.$
(2)

Layer 2 is assumed to be stretched or compressed in the $x$ (diagonal)
direction. As a result, the vectors connecting next nearest-neighbor Cu ions
in this layer have coordinates $\mathbf{d}^{\prime}_{1}=d(1-\delta,\,0)$ and
$\mathbf{d}^{\prime}_{2}=d(0,\,1)$, where parameter $\delta$ describes the
strength of the stretching ($\delta<0$) or compression ($\delta>0$). The
vectors connecting nearest-neighbor Cu ions of the layer 2 become
$\mathbf{a}^{\prime}_{1}=(\mathbf{d}^{\prime}_{1}-\mathbf{d}^{\prime}_{2})/2$
and
$\mathbf{a}^{\prime}_{2}=(\mathbf{d}^{\prime}_{1}+\mathbf{d}^{\prime}_{2})/2$.
The positions of copper and oxygen ions in the layer 2 are

$\displaystyle\mathbf{r}^{d}_{2\mathbf{n}}$ $\displaystyle\equiv$
$\displaystyle\mathbf{r}_{2\mathbf{n}}=n\mathbf{a}^{\prime}_{1}+m\mathbf{a}^{\prime}_{2}\,,\;$
$\displaystyle\mathbf{r}^{p_{x}}_{2\mathbf{n}}$ $\displaystyle=$
$\displaystyle\mathbf{r}_{2\mathbf{n}}+\frac{1}{2}\mathbf{a}^{\prime}_{1}\,,\;\mathbf{r}^{p_{y}}_{2\mathbf{n}}=\mathbf{r}_{2\mathbf{n}}+\frac{1}{2}\mathbf{a}^{\prime}_{2}\,.$
(3)

For a rational value of $1-\delta=N_{1}/N_{2}$ ($N_{1}$ and $N_{2}$ are co-
prime positive integers), the system has a superstructure. To satisfy the
periodicity conditions, we must include two adjacent chains of copper and
oxygen ions into superlattice cell aligned in the diagonal ($x$) direction. As
a result, the superlattice cell will contain $2N_{1}$ copper ions of layer 1
and $2N_{2}$ copper ions of layer 2. It also will contain $4N_{1}$ oxygen ions
of layer 1 and $4N_{2}$ oxygen ions of layer 2. The superlattice vectors are
$\mathbf{R}_{1}=N_{1}\mathbf{d}_{1}=N_{2}(1-\delta)\mathbf{d}_{1}$ and
$\mathbf{R}_{2}=2\mathbf{d}_{2}$. The superlattice Brillouin zone (reduced
Brillouin zone, RBZ) has a shape of the rectangle with the size
$2\pi/(N_{1}d)$ in the $x$ direction and $2\pi/d$ the $y$ direction (see the
right panel of Fig. 1).

##  III Tight-binding model of mismatched CuO2 bilayer

Figure 2: Densities of states at the Fermi level for unmismatched layers as
functions of doping $x$, calculated for three different values of $t_{2p}$.
Other model parameters are: $\varepsilon_{x^{2}-y^{2}}=0$,
$\varepsilon_{p_{x}}=\varepsilon_{p_{y}}=-0.9$ eV, $\varepsilon_{a}=6.5$ eV,
$t_{1p}=1.6$ eV, $t_{0}=0.27$ eV, and $t^{aa}_{0}=0.75$ eV. In order to avoid
logarithmic singularities, the densities of states are averaged over the
energy range $\Delta E=7$ meV.

To describe electronic properties of mismatched CuO2 bilayer, we use the four-
band tight-binding model. The model includes two types of $p$ oxygen orbitals:
$p_{x}$ orbitals of oxygen ions located in positions
$\mathbf{r}^{p_{x}}_{i\mathbf{n}}$ and $p_{y}$ orbitals of oxygen ions located
in positions $\mathbf{r}^{p_{y}}_{i\mathbf{n}}$ . Model also includes two
types of $d$ orbitals of copper located in positions
$\mathbf{r}^{d}_{i\mathbf{n}}$: these are $x^{2}-y^{2}$ and ‘axial’ ($a$)
orbitals. We do not specify here the nature of the latter $a$ orbital, and
only mention its circular symmetry in $xy$ plane. It can be, e.g., $s$ or
$3z^{2}-r^{2}$ orbital of copper, or their superposition. We include this
orbital into the model because it gives the largest effect into interlayer
hybridization Andersen et al. (1995).

Figure 3: Densities of states as functions of energy calculated for nine
superstructures with $N_{1}=36$ and $N_{2}$ varying from $37$ to $45$. Model
parameters correspond to Fig. 2 with $t_{2p}=1.5$ eV. In order to avoid
logarithmic singularities, the densities of states are averaged over the
energy range $\Delta E=7$ meV. The dashed green curve in each plot corresponds
to unmismatched layers.

The Hamiltonian under study can be written as

$H=H_{1}+H_{2}+H_{12}\,.$ (4)

The terms $H_{1,2}$ correspond to each individual layer, while $H_{12}$
describes the interlayer hybridization. Let us consider first intralayer
terms. We use the simplest intralayer tight-binding Hamiltonian including only
nearest neighbor hoppings of electrons between nearest copper and oxygen ions.
The Hamiltonian $H_{i}$ ($i=1,\,2$) reads

$H_{i}=\sum_{\mathbf{n}A\sigma}\varepsilon_{A}d^{{\dagger}}_{\mathbf{n}iA\sigma}d^{\phantom{{\dagger}}}_{\mathbf{n}iA\sigma}+\\!\sum_{\langle\mathbf{nm}\rangle\atop
dp\sigma}\\!\left(t^{dp}_{i\mathbf{nm}}d^{{\dagger}}_{\mathbf{n}id\sigma}d^{\phantom{{\dagger}}}_{\mathbf{m}ip\sigma}+H.c.\right),$
(5)

where the subscript $A=\\{d,\,p\\}$ (with $d=x^{2}-y^{2},\,a$ and
$p=p_{x},\,p_{y}$) denotes the considered orbitals, and
$d^{{\dagger}}_{\mathbf{n}iA\sigma}$ and
$d^{\phantom{{\dagger}}}_{\mathbf{n}iA\sigma}$ are the creation and
annihilation operators of electron located in unit cell $\mathbf{n}$ of layer
$i$ having orbital index $A$ and spin projection $\sigma$. The first term in
(5) includes the local energies of electrons in different orbitals, while the
last term describes the nearest-neighbor hopping. The hopping amplitudes
$t^{dp}_{i\mathbf{nm}}$ are position dependent. For a given direction of
hopping $\mathbf{m}\to\mathbf{n}$, the hopping amplitudes
$t^{dp}_{i\mathbf{nm}}$ can we written in the form of $2\times 2$ matrix. The
relationships between elements of this matrix can be easily obtained from the
symmetry analysis of orbitals’ wave functions Slater and Koster (1954). Let us
introduce the Fourier transformed electron operators
$d^{\phantom{{\dagger}}}_{\mathbf{k}iA\sigma}$ and Fourier transformed hopping
amplitudes $t^{dp}_{i\mathbf{k}}$. Taking into account the symmetry of the
orbitals under study Slater and Koster (1954), one can write for layer 1

$t^{dp}_{1\mathbf{k}}=\left(\begin{array}[]{cc}-t_{1p}e^{-i\mathbf{ka}_{1}/2}&t_{1p}e^{-i\mathbf{ka}_{2}/2}\\\
-t_{2p}e^{-i\mathbf{ka}_{1}/2}&-t_{2p}e^{-i\mathbf{ka}_{2}/2}\\\
\end{array}\right),$ (6)

where we introduce two independent Slater–Koster parameters describing the
hopping between $x^{2}-y^{2}$ and $p$ orbitals ($t_{1p}$) and between $a$ and
$p$ orbitals ($t_{2p}$). The similar formula for the layer 2 can be obtained
from Eq. (6) with the change $\mathbf{a}_{1,2}\to\mathbf{a}^{\prime}_{1,2}$.
Strictly speaking, the values of parameters $t_{1p}$ and $t_{2p}$ for layer 2
are different from those for layer 1. In further analysis, however, we neglect
this difference due to the smallness of the lattice distortion. Following Ref.
Andersen et al., 1995, we use the following values of the tight-binding
parameters: $\varepsilon_{x^{2}-y^{2}}=0$,
$\varepsilon_{p_{x}}=\varepsilon_{p_{y}}=-0.9$ eV, $t_{1p}=1.6$ eV,
$t_{2p}=1.5$ eV. In our simulations, we use different values of the parameter
$\varepsilon_{a}$ describing the on-site energy of $a$ orbital, taking both
negative (almost $3z^{2}-r^{2}$ orbital) and positive (almost $s$ orbital)
values. We also performed simulations for different values of the parameter
$t_{2p}$ describing hybridization between $a$ orbitals of copper and $p$
oxygen orbitals. Independent on the values of $\varepsilon_{a}$ and $t_{2p}$,
the density of states of individual layer has only one van Hove peak
corresponding to the change of the Fermi surface topology. The position of
this peak depends on $\varepsilon_{a}$ and $t_{2p}$. In present study, we
always focus on the energy range close to the van Hove peak. The obtained
results are similar for different values of $\varepsilon_{a}$ and $t_{2p}$.

Let us consider now the interlayer hopping. According to Ref. Andersen et al.,
1995, the largest interlayer hopping amplitudes are between $a$ orbitals in
two different layers, and between $a$ and $p_{x,y}$ orbitals of different
layers. As a result, Hamiltonian $H_{12}$ takes the form

$\displaystyle H_{12}$ $\displaystyle=$
$\displaystyle\\!\\!\sum_{\mathbf{nm}p\sigma}\\!\\!\left(t^{ap}_{\bot\mathbf{nm}}d^{{\dagger}}_{\mathbf{n}2a\sigma}d^{\phantom{{\dagger}}}_{\mathbf{m}1p\sigma}+t^{pa}_{\bot\mathbf{nm}}d^{{\dagger}}_{\mathbf{n}2p\sigma}d^{\phantom{{\dagger}}}_{\mathbf{m}1a\sigma}\right)+$
(7)
$\displaystyle\\!\sum_{\mathbf{nm}\sigma}\\!t^{aa}_{\bot\mathbf{nm}}d^{{\dagger}}_{\mathbf{n}2a\sigma}d^{\phantom{{\dagger}}}_{\mathbf{m}1a\sigma}+H.c.$

The value of hopping amplitude depends on the mutual positions of two
orbitals. In our simulations, we tried different parametrizations of
interlayer hopping amplitudes depending on the possible symmetry of $a$
orbital and obtain qualitatively similar results. If $a$ orbital is assumed to
be rather of $s$ symmetry, the amplitude $t^{ap_{x,y}}_{\bot\mathbf{nm}}$ of
electron hopping from oxygen ion in layer $1$ located in the position
$\mathbf{r}^{p_{x,y}}_{1\mathbf{m}}$ to copper ion in layer $2$ located in the
position $\mathbf{r}^{d}_{2\mathbf{n}}$ is described by the following
Slater–Koster formula

$\displaystyle
t^{ap_{x,y}}_{\bot\mathbf{nm}}=\frac{(\mathbf{r}^{d}_{2\mathbf{n}}-\mathbf{r}^{p_{x,y}}_{1\mathbf{m}})\mathbf{e}_{x,y}}{\sqrt{c^{2}+(\mathbf{r}^{d}_{2\mathbf{n}}-\mathbf{r}^{p_{x,y}}_{1\mathbf{m}})^{2}}}V_{\sigma}(\mathbf{r}^{d}_{2\mathbf{n}}-\mathbf{r}^{p_{x,y}}_{1\mathbf{m}})\,,$
(8)

where $\mathbf{e}_{x,y}$ is the unit vector in the direction of
$\mathbf{a}_{1,2}$ and the function

$V_{\sigma}(r)=t_{0}\sqrt{1+8c^{2}/d^{2}}\,e^{-\frac{r-\sqrt{c^{2}+d^{2}/8}}{r_{0}}}.$
(9)

In the expression above, $t_{0}$ is the largest interlayer hopping amplitude
between $a$ and $p$ orbitals and $c$ is the interlayer distance. We choose
Bianconi et al. (1996d) $c=3.35$ Å(note that the unit cell of the bulk BSCCO
contains eight CuO2 layers). The parameter $r_{0}$ defines how fast the
function $V_{\sigma}(r)$ decays with the distance between orbitals (we choose
$r_{0}/d=0.19$). Following Ref. Andersen et al., 1995, we take $t_{0}=0.27$
eV. For the hopping amplitude between $a$ orbitals in different layers, we
have

$\displaystyle
t^{aa}_{\bot\mathbf{nm}}=V_{0}(\mathbf{r}^{d}_{2\mathbf{n}}-\mathbf{r}^{d}_{1\mathbf{m}}),\;\;V_{0}(\mathbf{r})=t^{aa}_{0}e^{-\frac{r-c}{r^{aa}_{0}}}\,,$
(10)

where we choose $t^{aa}_{0}=0.75$ eV and $r^{aa}_{0}/d=0.3$.

Following the Slater–Koster formalism Slater and Koster (1954), we tried also
a bit complicated parametrization of interlayer hopping amplitudes, which
corresponds to $a$ orbital rather of $3z^{2}-r^{2}$ symmetry. Again, the
functions similar to $V_{\sigma}(r)$ and $V_{0}(r)$, contain factors depending
on directional cosines and the factors describing exponential decay with the
distance between ions, but now the hopping amplitudes contain larger number of
such functions. Details of this parametrization can be found elsewhere. Here,
we present the results corresponding to the first parametrization, but let us
notice again that qualitatively the results are independent of the type of the
hopping amplitude parametrization.

##  IV Results

Figure 4: Densities of states at Fermi level as functions of doping $x$
calculated for nine superstructures with $N_{1}=36$ and $N_{2}$ varying form
$37$ to $45$. Model parameters correspond to Fig. 2 with $t_{2p}=1.5$ eV. In
order to avoid logarithmic singularities, the densities of states are averaged
over the energy range $\Delta E=7$ meV.

Figure 5: Energy spectra calculated for superstructures with $N_{1}=36$ and
$N_{2}=36$ (top panel, no mismatch), $N_{2}=38$ (middle panel), and $N_{2}=42$
(bottom panel). The spectra are calculated along the contour shown in Fig. 1.
The distance between $M$ and $X$ points is enhanced by factor $5$ for better
view. Model parameters correspond to Fig. 2 with $t_{2p}=1.5$ eV.

Here, we demonstrate the evolution of the density of states and energy spectra
for bilayers at different values of the lattice mismatch. We consider the
superstructures with fixed $N_{1}=36$ and different $N_{2}$ ranging from
$N_{2}=36$ (no mismatch) to $N_{2}=45$ (compressive strain). We suggest that
apical orbital $a$ is rather of $s$ orbital of copper, so its energy level
lies above $x^{2}-y^{2}$ orbital, and $\varepsilon_{a}=6.5$ eV. In this case,
the undoped bilayer has $5$ electrons inside elementary unit cell of each
layer (four electrons on $p_{x,y}$ orbitals plus one electron on
$d_{x^{2}-y^{2}}$ orbital). We are interested in the hole doped bilayers with
doping $0<x=5-n<1$, where $n$ is the number of electrons per CuO2 unit cell.
Let us note here that $x$ differs from the oxygen content $y$ in the chemical
formula of Bi2Sr2CaCu2O8+y compound. For example, according to Ref. Poccia et
al., 2011, we have $x=y-0.05$ for some $y$. In the interesting energy range,
the density of states of the individual layer has one van Hove singularity.
For unmismatched layers, the interlayer hopping terms split this van Hove
singularity into two. Their energy positions depend on the model parameters.
Thus, increasing the hopping amplitude between $a$ and $p_{x,y}$ orbitals,
$t_{2p}$, shifts these peaks to larger energies (smaller doping) in agreement
with Ref. Pavarini et al., 2001. For example, in Fig. 2, we show the densities
of states at Fermi level as functions of doping $x$ calculated for three
different values of $t_{2p}$. Below, we consider the case of $t_{2p}=1.5$ eV
as more realistic from the viewpoint of the positions of van Hove
singularities in real compounds.

Figure 6: Energy spectra calculated for superstructures with $N_{1}=4$ and
$N_{2}=4$ (top panel, no mismatch), $N_{2}=5$ (middle panel), and $N_{2}=6$
(bottom panel). The spectra are calculated along the contour shown in Fig. 1.
Model parameters correspond to Fig. 2 with $t_{2p}=1.5$ eV.

The main message of our work is that the mismatch between two CuO2 layers
generates extra van Hove singularities. For example, in Fig. 3, we show the
densities of states as function of energy calculated for nine superstructures
with $N_{1}=36$ and $N_{2}$ varying form $37$ to $45$. We clearly see extra
van Hove singularities approximately located between to peaks corresponding to
unmismatched layers. It is interesting, that the number of these extra van
Hove singularities exhibits a counterintuitive behavior with the change of
$N_{2}$: the number of peaks increases when $N_{2}$ decreases from $N_{2}=45$
to $N_{2}=37$. For better comparison, we show in a single plot, Fig. 4, the
densities of states at Fermi level as functions of doping calculated for all
nine superstructures with different $N_{2}$ and fixed $N_{1}=36$.

For better understanding of the origin of the extra van Hove singularities, we
calculated the energy spectra for several superstructures (see Figs. 5 and 6).
The spectra are calculated along the triangular contour connecting the most
symmetrical points $\Gamma$, $M$, and $X$ of the reduced Brillouin zone, see
the right panel of Fig. 1. For Fig. 5, there are a lot of bands (about $100$)
inside the energy range shown, $1.9$ eV$<E<2.25$ eV. This energy range
includes two van Hove singularities of unmismatched layers (upper panel of
Fig. 5), located at $E^{1}_{vHs}\cong 1.93$ eV and $E^{2}_{vHs}\cong 2.19$ eV.
The analysis shows that the origin of these two van Hove singularities are two
saddle points of certain bands located at $M$ point of the reduced Brillouin
zone. In the energy range approximately between $E^{1}_{vHs}$ and
$E^{2}_{vHs}$, the band crossing occurs in the region near the line connecting
$M$ and $X$ points of the RBZ, as it can be seen in top panel of Fig. 5. Such
a picture is observed for unmismatched layers. When the layer 2 becomes
stretched or compressed in comparison to layer 1, we observe a band flattening
and band splitting in the region near the line connecting $M$ and $X$ points
as it is seen in middle and bottom panels in Fig. 5. This band flattening
gives rise to extra van Hove singularities with energies lying approximately
between $E^{1}_{vHs}$ and $E^{2}_{vHs}$. Situation here is very similar to
that observed in twisted bilayer graphene Rozhkov et al. (2016), where the
band flattening is observed for large superstructures, when the twist angle
$\theta$ decreases down to critical value $\theta_{c}\sim 1^{\circ}$, while
the moiré period increases as $L\propto 1/\theta$. In our case, however, the
band flattening occurs only in the $k_{x}$-direction, since the band folding
is performed only in this direction in the momentum space. As we can see from
Fig. 5, the bands remain dispersive in the $k_{y}$-direction.

Note that the recent studies of 2D systems and moiré superlattices (e.g., in
twisted bilayer graphene) demonstrate that such systems often exhibit higher-
order VHS characterized by the power-law divergence in the density of states
in contrast to the logarithmic divergence for conventional VHS Yuan et al.
(2019); Isobe and Fu (2019); Classen et al. (2020); Guerci et al. (2022). Near
such VHS the dispersion is flatter than near conventional ones. To illustrate
such situation in our case, we show in Fig. 7 a zoomed version of the highest
DOS peak for the 36/43 structure. At a higher resolution, we can see that this
peak splits into two ones. One of them has the logarithmic divergence, whereas
another one exhibits the power-law divergence ($E^{-1/4}$) implying the
existence of a higher-order VHS.

Figure 7:  The largest DOS peak for the 36/43 structure. At a higher
resolution, we can see that this peak splits into two ones. One of them (at
$E_{1}$) has the logarithmic divergence, whereas another one (at $E_{2}$)
exhibits the power-law divergence ($E^{-1/4}$) implying the existence of a
higher-order VHS. Different prefactors of the power-law divergence and the
exponent of the power law agree with the predictions of Refs. Yuan et al.,
2019 and Isobe and Fu, 2019. The stepwise change in the DOS corresponds to the
band edge.

The spectra shown in Fig. 5 contain a lot of bands. For better illustration,
we calculated also spectra for much smaller superlattices. Namely, we
calculated spectra for superstructures with $N_{1}=4$ and $N_{2}=4,\,5$, and
$6$, see Fig. 6. Since, in these cases, the superlattice cell contains much
smaller number of atoms, the energy range of interest (approximately between
$E^{1}_{vHs}$ and $E^{2}_{vHs}$) contains much smaller number of bands. For
superstructure with $N_{1}=4$ and $N_{2}=4$ (no mismatch, top panel of Fig.
6), we see two parabolic bands with minima at $M$ points producing two van
Hove singularities at energies $E^{1}_{vHs}$ and $E^{2}_{vHs}$. By calculating
the spectrum around small circle centered at $M$ point, we have checked that
$M$ point is indeed a saddle point for corresponding bands. In the momentum
range near the line connecting $M$ and $X$ points, we see several band
crossings. When we compress layer $2$ (superstructures with $N_{2}=5$ and $6$,
middle and bottom panels in Fig. 6), this produces the band splittings, which
give rise to extra van Hove singularities.

##  V Conclusions

In conclusion, we have studied the model of double-layer BSCCO cuprate
superconductor with mismatched CuO2 bilayers in diagonal direction. The
mismatch between two layers produces the superstructure in the system. This
gives rise to several extra van Hove singularities located approximately
between energies $E^{1}_{vHs}$ and $E^{2}_{vHs}$ corresponding to two van Hove
singularities of unmismatched layers. The number of these peaks and their
energy positions depends on the quasi-commensurate superstructures of the
devil’s staircase in strained bilayers.

We have found that these extra van Hove singularities are closely related to
the flattening and splitting of the bands inside the reduced Brillouin zone of
mismatched layers. The effect of the band flattening is similar to that
observed in twisted bilayer graphene Rozhkov et al. (2016), where the large
Fermi velocity reduction occurs for superstructures with small twist angles
Lopes dos Santos et al. (2012); Sboychakov et al. (2015). In our case,
however, the band flattening appears not in the whole Brillouin zone, but only
in the direction of the band folding ($k_{x}$-direction). In the perpendicular
direction, the bands remain dispersive.

Note here that our approach relies on several approximations. For example, (i)
only nearest neighbor hoppings are considered; (ii) the rotation of the
orbitals in the stretched layer is not considered, and (iii) the hopping
parameters in the stretched layer are assumed to remain the same as those of
the non-stretched layer. Of course, these approximations make the model rather
oversimplified and the extension of this model should indeed make it more
realistic. However, just the structure of our calculations implies that the
extra van Hove singularities will not disappear and a certain partial band
flattening should remain, but maybe in some modified manner. Indeed, our
additional calculations performed using different values of the hopping
integrals in each layer lead to the qualitatively similar results.

Let us also mention that DOS at the Fermi level governs to a large extent
various thermodynamic properties such as charge compressibility, spin
susceptibility, and specific heat. The well-known effect of the large DOS at
the Fermi level on superconductivity can be here additionally enhanced since
the large DOS may result in a strong screening of the repulsive interaction.
This effect is the most clearly pronounced at the van Hove singularities,
especially at strongly divergent high-order ones.

The present results show that BSCCO is an intrinsic multi-band system, where
at high doping, the chemical potential is close to a Lifshitz transition
appearing at VHS at the metal to superconductor transition (see also Ref.
Mazziotti et al., 2021a). Therefore, superconductivity in BSCCO is not a
single-band unconventional superconducting phase but a multi-gap
superconductor at a Fano–Feshbach resonance between gaps in the BCS regime and
in the BCS–BEC crossover as described by Perali et al. in cuprates Perali et
al. (1996, 2012); Chen et al. (2012); Salasnich et al. (2019) and appearing in
diborides Agrestini et al. (2004), oxide interfaces Bianconi et al. (2014) and
in room temperature hydride superconductors Mazziotti et al. (2021b, c).

Our results are supported by the recent evidence by a high resolution STM
experiment of nanoscale phase separation in a BSCCO sample doped up to $0.27$
holes/Cu site Tromp et al. (2022), reaching the high doping range of the high
order VHS singularities, which was predicted by theory to appear where the
chemical potential approaches a high order VHS singularity near a Lifshitz
transition in strongly correlated multiband systems M.Yu. Kagan et al. (2021);
Rakhmanov et al. (2017); Kugel et al. (2008a); Bianconi et al. (2015); Kugel
et al. (2008b). The position of the Fermi level with respect to the high order
VHS peaks can be controlled by the gate voltage using modern technologies.
Therefore, this makes mismatched bilayer cuprates to be an efficient system
for fine tuning of the peaks of the electron density of states with respect to
the Fermi level. Note that this effect is the most clearly pronounced at
rather small lattice mismatch, which produces numerous additional peaks in the
density of states.

## Acknowledgments

A.O.S. acknowledge the support of the Russian Science Foundation (project No.
22-22-00464) in the part concerning numerical calculations. K.I.K. acknowledge
the support of the Russian Science Foundation (project No. 20-62-46047) in the
part concerning the data analysis. We are grateful to the Joint Supercomputer
Center of the Russian Academy of Sciences (JSCC RAS) for the provided
computational resources. A.B. acknowledges the support of the Ministry of
Science and Higher Education of the Russian Federation (Agreement No.
075-15-2021-1352).

## References

  * Liu et al. (2016) Y. Liu, N. O. Weiss, X. Duan, H.-C. Cheng, Y. Huang, and X. Duan, “Van der Waals heterostructures and devices,” Nat. Rev. Mater. 1, 16042 (2016).
  * Rozhkov et al. (2016) A. Rozhkov, A. Sboychakov, A. Rakhmanov, and F. Nori, “Electronic properties of graphene-based bilayer systems,” Phys. Rep. 648, 1 (2016).
  * Weston et al. (2020) A. Weston, Y. Zou, V. Enaldiev, A. Summerfield, N. Clark, V. Zólyomi, A. Graham, C. Yelgel, S. Magorrian, M. Zhou, et al., “Atomic reconstruction in twisted bilayers of transition metal dichalcogenides,” Nat. Nanotechnol. 15, 592 (2020).
  * Zhao et al. (2021) W.-M. Zhao, L. Zhu, Z. Nie, Q.-Y. Li, Q.-W. Wang, L.-G. Dou, J.-G. Hu, L. Xian, S. Meng, and S.-C. Li, “Moiré enhanced charge density wave state in twisted 1T-TiTe2/1T-TiSe2 heterostructures,” Nat. Mater. 21, 284 (2021).
  * Tong et al. (2017) Q. Tong, H. Yu, Q. Zhu, Y. Wang, X. Xu, and W. Yao, “Topological mosaics in moiré superlattices of van der Waals heterobilayers,” Nat. Phys. 13, 356 (2017).
  * Ren et al. (2020) Y.-N. Ren, Y. Zhang, Y.-W. Liu, and L. He, “Twistronics in graphene-based van der Waals structures,” Chin. Phys. B 29, 117303 (2020).
  * Krockenberger et al. (2021) Y. Krockenberger, A. Ikeda, and H. Yamamoto, “Atomic stripe formation in infinite-layer cuprates,” ACS Omega 6, 21884 (2021).
  * Park et al. (2021) H. Park, J. Oh, B. Song, and B. Kang, “Structural suitability of GdFeO3 as a magnetic buffer layer for GdBa2Cu3O7-x superconducting thin films,” Progr. Supercond. Cryog. 23, 14 (2021).
  * Conradson et al. (2020) S. D. Conradson, T. H. Geballe, C.-Q. Jin, L.-P. Cao, A. Gauzzi, M. Karppinen, G. Baldinozzi, W.-M. Li, E. Gilioli, J. M. Jiang, et al., “Nonadiabatic coupling of the dynamical structure to the superconductivity in YSr2Cu2.75Mo0.25O7.54 and Sr2CuO3.3,” PNAS 117, 33099 (2020).
  * Sederholm et al. (2021) L. Sederholm, S. D. Conradson, T. H. Geballe, C.-Q. Jin, A. Gauzzi, E. Gilioli, M. Karppinen, and G. Baldinozzi, “Extremely overdoped superconducting cuprates via high pressure oxygenation methods,” Condens. Matter 6 (2021).
  * Jang et al. (2017) H. Jang, S. Asano, M. Fujita, M. Hashimoto, D. H. Lu, C. A. Burns, C.-C. Kao, and J.-S. Lee, “Superconductivity-insensitive order at $q\sim 1/4$ in electron-doped cuprates,” Phys. Rev. X 7, 041066 (2017).
  * Bianconi et al. (1996a) A. Bianconi, N. L. Saini, A. Lanzara, M. Missori, T. Rossetti, H. Oyanagi, H. Yamaguchi, K. Oka, and T. Ito, “Determination of the local lattice distortions in the Cu${\mathrm{O}}_{2}$ plane of L${\mathrm{a}}_{1.85}$S${\mathrm{r}}_{0.15}$Cu${\mathrm{O}}_{4}$,” Phys. Rev. Lett. 76, 3412 (1996a).
  * Saini et al. (2003) N. Saini, H. Oyanagi, T. Ito, V. Scagnoli, M. Filippi, S. Agrestini, G. Campi, K. Oka, and A. Bianconi, “Temperature dependent local Cu-O displacements from underdoped to overdoped La-Sr-Cu-O superconductor,” Eur. Phys. J. B 36, 75 (2003).
  * Bianconi et al. (1996b) A. Bianconi, M. Lusignoli, N. L. Saini, P. Bordet, A. Kvick, and P. G. Radaelli, “Stripe structure of the ${\mathrm{CuO}}_{2}$ plane in ${\mathrm{Bi}}_{2}$${\mathrm{Sr}}_{2}$${\mathrm{CaCu}}_{2}$${\mathrm{O}}_{8+\mathit{y}}$ by anomalous x-ray diffraction,” Phys. Rev. B 54, 4310 (1996b).
  * Bianconi et al. (1996c) A. Bianconi, N. L. Saini, T. Rossetti, A. Lanzara, A. Perali, M. Missori, H. Oyanagi, H. Yamaguchi, Y. Nishihara, and D. H. Ha, “Stripe structure in the ${\mathrm{CuO}}_{2}$ plane of perovskite superconductors,” Phys. Rev. B 54, 12018 (1996c).
  * Di Castro et al. (2000) D. Di Castro, G. Bianconi, M. Colapietro, A. Pifferi, N. Saini, S. Agrestini, and A. Bianconi, “Temperature dependent local Cu-O displacements from underdoped to overdoped La-Sr-Cu-O superconductor,” Eur. Phys. J. B 18, 617 (2000).
  * Bianconi et al. (2000) A. Bianconi, G. Bianconi, S. Caprara, D. D. Castro, H. Oyanagi, and N. L. Saini, “The stripe critical point for cuprates,” J. Phys : Condens. Matter 12, 10655 (2000).
  * Agrestini et al. (2003) S. Agrestini, N. Saini, G. Bianconi, and A. Bianconi, “The strain of CuO2 lattice: the second variable for the phase diagram of cuprate perovskites,” J. Phys. A: Math. Gen. 36, 9133 (2003).
  * Wiegers (1996) G. Wiegers, “Misfit layer compounds: Structures and physical properties,” Prog. Solid State Chem. 24, 1 (1996).
  * Janssen et al. (2004) T. Janssen, A. Janner, A. Looijenga-Vos, and P. de Wolff, _International Tables for Crystallography, Mathematical, Physical and Chemical Tables, E. Prince, ed._ (Wiley, New York, 2004), vol. C of _IUCr Series_ , chap. in 9.8. Incommensurate and Commensurate Modulated Structures, 3rd ed.
  * Etrillard et al. (2000) J. Etrillard, P. Bourges, and C. T. Lin, “Incommensurate composite structure of the superconductor ${\mathrm{Bi}}_{2}{\mathrm{Sr}}_{2}{\mathrm{CaCu}}_{2}{\mathrm{O}}_{8+\delta}$,” Phys. Rev. B 62, 150 (2000).
  * Grebille et al. (2001) D. Grebille, H. Leligny, and O. Pérez, “Comment on “Incommensurate composite structure of the superconductor Bi2Sr2CaCu2O8+δ”,” Phys. Rev. B 64, 106501 (2001).
  * Bianconi (1994) A. Bianconi, “On the Fermi liquid coupled with a generalized Wigner polaronic CDW giving high Tc superconductivity,” Solid State Commun. 91, 1 (1994).
  * Beskrovnyy and Jirák (2012) A. Beskrovnyy and Z. Jirák, “Structural modulation in Bi2Sr2Ca0.4Y0.6Cu2O8+δ,” J. Phys. Conf. Ser. 340 (2012).
  * Pérez et al. (1997) O. Pérez, H. Leligny, D. Grebille, J. M. Grenêche, P. Labbé, D. Groult, and B. Raveau, “Disorder phenomena in the incommensurate compound Bi2+xSr3-xFe2O9+δ,” Phys. Rev. B 55, 1236 (1997).
  * Sedykh et al. (2003) V. Sedykh, V. Shekhtman, I. Smirnova, B. Bagautdinov, E. Suvorov, and A. Dubovitskii, “About specific features of the structure modulation in the Bi-ferrate compounds isostructural with Bi2Sr2CaCu2O8,” Physica C 390, 311 (2003).
  * Lambert et al. (2001) S. Lambert, H. Leligny, and D. Grebille, “Three forms of the misfit layered cobaltite [Ca2CoO3][CoO2]1.62. A 4D structural investigation,” J. Solid State Chem. 160, 322 (2001).
  * Grebille et al. (2007) D. Grebille, H. Muguerra, O. Pérez, E. Guilmeau, H. Rousselière, and R. Funahashi, “Superspace crystal symmetry of thermoelectric misfit cobalt oxides and predicted structural models,” Acta Cryst. B 63, 373 (2007).
  * Kim et al. (2012) J.-Y. Kim, J.-I. Kim, S.-M. Choi, Y. Soo Lim, W.-S. Seo, and H. Jin Hwang, “Nanostructured thermoelectric cobalt oxide by exfoliation/restacking route,” J. Appl. Phys. 112, 113705 (2012).
  * Gavrichkov et al. (2019) V. A. Gavrichkov, Y. Shan’ko, N. G. Zamkova, and A. Bianconi, “Is there any hidden symmetry in the stripe structure of perovskite high-temperature superconductors?,” J. Phys. Chem. Lett. 10, 1840 (2019).
  * Mazza et al. (2022) A. R. Mazza, X. Gao, D. J. Rossi, B. L. Musico, T. W. Valentine, Z. Kennedy, J. Zhang, J. Lapano, V. Keppens, R. G. Moore, et al., “Searching for superconductivity in high entropy oxide Ruddlesden–Popper cuprate films,” J. Vac. Sci. Tech. A 40, 013404 (2022).
  * Kaneko et al. (2004) S. Kaneko, Y. Shimizu, K. Akiyama, T. Ito, M. Mitsuhashi, S. Ohya, K. Saito, H. Funakubo, and M. Yoshimoto, “Modulation derived satellite peaks in x-ray reciprocal mapping on bismuth cuprate superconductor film,” Appl. Phys. Lett. 85, 2301 (2004).
  * Poccia et al. (2011) N. Poccia, G. Campi, M. Fratini, A. Ricci, N. L. Saini, and A. Bianconi, “Spatial inhomogeneity and planar symmetry breaking of the lattice incommensurate supermodulation in the high-temperature superconductor Bi2Sr2CaCu2O8+y,” Phys. Rev. B 84, 100504 (2011).
  * Poccia et al. (2020) N. Poccia, S. Y. F. Zhao, H. Yoo, X. Huang, H. Yan, Y. S. Chu, R. Zhong, G. Gu, C. Mazzoli, K. Watanabe, et al., “Spatially correlated incommensurate lattice modulations in an atomically thin high-temperature ${\mathrm{Bi}}_{2.1}{\mathrm{Sr}}_{1.9}\mathrm{Ca}{\mathrm{Cu}}_{2.0}{\mathrm{O}}_{8+y}$ superconductor,” Phys. Rev. Materials 4, 114007 (2020).
  * Zhang et al. (2016) Y. Zhang, C. Hu, Y. Hu, L. Zhao, Y. Ding, X. Sun, A. Liang, Y. Zhang, S. He, D. Liu, et al., “In situ carrier tuning in high temperature superconductor Bi2Sr2CaCu2O8+δ by potassium deposition,” Sci. Bull. 61, 1037 (2016).
  * Bianconi and Poccia (2012) A. Bianconi and N. Poccia, “Superstripes and complexity in high-temperature superconductors,” J. Supercond. Novel Magn. 25, 1403 (2012).
  * Bianconi et al. (1993) A. Bianconi, S. D. Longa, M. Missori, I. Pettiti, M. Pompa, and A. Soldatov, “Structure of the different Cu sites in the corrugated CuO2 plane in high $T_{c}$ superconductors,” Jpn. J. Appl. Phys. 32, 578 (1993).
  * Alldredge et al. (2008) J. Alldredge, J. Lee, K. McElroy, M. Wang, K. Fujita, Y. Kohsaka, C. Taylor, H. Eisaki, P. Uchida, S. Hirschfeld, and J. Davis, “Evolution of the electronic excitation spectrum with strongly diminishing hole density in superconducting Bi2Sr2CaCu2O8+δ,” Nat. Phys. 4, 319 (2008).
  * Pasupathy et al. (2008) A. N. Pasupathy, A. Pushp, K. K. Gomes, C. V. Parker, J. Wen, Z. Xu, G. Gu, S. Ono, Y. Ando, and A. Yazdani, “Electronic origin of the inhomogeneous pairing interaction in the high-$T_{c}$ superconductor Bi2Sr2CaCu2O8+δ,” Science 320, 196 (2008).
  * Gomes et al. (2007) K. K. Gomes, A. N. Pasupathy, A. Pushp, S. Ono, Y. Ando, and A. Yazdani, “Visualizing pair formation on the atomic scale in the high-$T_{c}$ superconductor Bi2Sr2CaCu2O8+δ,” Nature 447, 569 (2007).
  * Božović et al. (2016) I. Božović, X. He, J. Wu, and A. Bollinger, “Dependence of the critical temperature in overdoped copper oxides on superfluid density,” Nature 536, 309 (2016).
  * He et al. (2021) Y. He, S.-D. Chen, Z.-X. Li, D. Zhao, D. Song, Y. Yoshida, H. Eisaki, T. Wu, X.-H. Chen, D.-H. Lu, et al., “Superconducting fluctuations in overdoped Bi2Sr2CaCu2O8+δ,” Phys. Rev. X 11, 031068 (2021).
  * Drozdov et al. (2018) I. K. Drozdov, I. Pletikosić, C.-K. Kim, K. Fujita, G. Gu, J. C. S. Davis, P. Johnson, I. Božović, and T. Valla, “Phase diagram of Bi2Sr2CaCu2O8+δ revisited,” Nat. Commun. 9, 5210 (2018).
  * Li et al. (2018) X. Li, Y. Ding, C. He, W. Ruan, P. Cai, C. Ye, Z. Hao, L. Zhao, X. Zhou, Q. Wang, et al., “Quasiparticle interference and charge order in a heavily overdoped non-superconducting cuprate,” New J. Phys. 20, 063041 (2018).
  * Maier et al. (2020) T. A. Maier, S. Karakuzu, and D. J. Scalapino, “Overdoped end of the cuprate phase diagram,” Phys. Rev. Research 2, 033132 (2020).
  * M.Yu. Kagan et al. (2021) M.Yu. Kagan, K.I. Kugel, and A.L. Rakhmanov, “Electronic phase separation: Recent progress in the old problem,” Phys. Rep. 916, 1 (2021).
  * Rakhmanov et al. (2017) A. L. Rakhmanov, K. I. Kugel, M. Yu. Kagan, A. V. Rozhkov, and A. O. Sboychakov, “Inhomogeneous electron states in the systems with imperfect nesting,” JETP Letters 105, 806 (2017).
  * Kugel et al. (2008a) K. I. Kugel, A. L. Rakhmanov, A. O. Sboychakov, N. Poccia, and A. Bianconi, “Model for phase separation controlled by doping and the internal chemical pressure in different cuprate superconductors,” Phys. Rev. B 78, 165124 (2008a).
  * Bianconi et al. (2015) A. Bianconi, N. Poccia, A. O. Sboychakov, A. L. Rakhmanov, and K. I. Kugel, “Intrinsic arrested nanoscale phase separation near a topological Lifshitz transition in strongly correlated two-band metals,” Supercond. Sci. Technol. 28, 024005 (2015).
  * Kugel et al. (2008b) K. I. Kugel, A. L. Rakhmanov, A. O. Sboychakov, F. V. Kusmartsev, N. Poccia, and A. Bianconi, “A two-band model for the phase separation induced by the chemical mismatch pressure in different cuprate superconductors,” Supercond. Sci. Technol. 22, 014007 (2008b).
  * Massidda et al. (1988) S. Massidda, J. Yu, and A. Freeman, “Electronic structure and properties of Bi2Sr2CaCu2O8, the third high-$T_{c}$ superconductor,” Physica C 152, 251 (1988).
  * Nokelainen et al. (2020) J. Nokelainen, C. Lane, R. S. Markiewicz, B. Barbiellini, A. Pulkkinen, B. Singh, J. Sun, K. Pussi, and A. Bansil, “Ab initio description of the Bi2Sr2CaCu2O8+δ electronic structure,” Phys. Rev. B 101, 214523 (2020).
  * Hewat et al. (1989) E. Hewat, J. Capponi, and M. Marezio, “A model for the superstructure of Bi2Sr2CaCu2O8.2,” Physica C 157, 502 (1989).
  * Bak (1982) P. Bak, “Commensurate phases, incommensurate phases and the devil’s staircase,” Rep. Prog. Phys. 45, 587 (1982).
  * Andersen et al. (1995) O. Andersen, A. Liechtenstein, O. Jepsen, and F. Paulsen, “LDA energy bands, low-energy Hamiltonians, $t^{\prime}$, $t^{\prime\prime}$ $t_{\bot}(k)$, and $J_{\bot}$,” J. Phys. Chem. Solids 56, 1573 (1995).
  * Slater and Koster (1954) J. C. Slater and G. F. Koster, “Simplified LCAO method for the periodic potential problem,” Phys. Rev. 94, 1498 (1954).
  * Bianconi et al. (1996d) A. Bianconi, M. Lusignoli, N. L. Saini, P. Bordet, A. Kvick, and P. G. Radaelli, “Stripe structure of the CuO2 plane in Bi2Sr2CaCu2O8+y by anomalous x-ray diffraction,” Phys. Rev. B 54, 4310 (1996d).
  * Pavarini et al. (2001) E. Pavarini, I. Dasgupta, T. Saha-Dasgupta, O. Jepsen, and O. K. Andersen, “Band-Structure Trend in Hole-Doped Cuprates and Correlation with ${\mathit{T}}_{\mathit{c}\mathrm{max}}$,” Phys. Rev. Lett. 87, 047003 (2001).
  * Yuan et al. (2019) N. F. Q. Yuan, H. Isobe, and L. Fu, “Magic of high-order van Hove singularity,” Nat. Commun. 10, 5769 (2019).
  * Isobe and Fu (2019) H. Isobe and L. Fu, “Supermetal,” Phys. Rev. Research 1, 033206 (2019).
  * Classen et al. (2020) L. Classen, A. V. Chubukov, C. Honerkamp, and M. M. Scherer, “Competing orders at higher-order Van Hove points,” Phys. Rev. B 102, 125141 (2020).
  * Guerci et al. (2022) D. Guerci, P. Simon, and C. Mora, “Higher-order Van Hove singularity in magic-angle twisted trilayer graphene,” Phys. Rev. Research 4, L012013 (2022).
  * Lopes dos Santos et al. (2012) J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, “Continuum model of the twisted graphene bilayer,” Phys. Rev. B 86, 155449 (2012).
  * Sboychakov et al. (2015) A. O. Sboychakov, A. L. Rakhmanov, A. V. Rozhkov, and F. Nori, “Electronic spectrum of twisted bilayer graphene,” Phys. Rev. B 92, 075402 (2015).
  * Mazziotti et al. (2021a) M. V. Mazziotti, A. Valletta, R. Raimondi, and A. Bianconi, “Multigap superconductivity at an unconventional Lifshitz transition in a three-dimensional Rashba heterostructure at the atomic limit,” Phys. Rev. B 103, 024523 (2021a).
  * Perali et al. (1996) A. Perali, A. Bianconi, A. Lanzara, and N. L. Saini, “The gap amplification at a shape resonance in a superlattice of quantum stripes: A mechanism for high Tc,” Solid State Commun. 100, 181 (1996).
  * Perali et al. (2012) A. Perali, D. Innocenti, A. Valletta, and A. Bianconi, “Anomalous isotope effect near a 2.5 Lifshitz transition in a multi-band multi-condensate superconductor made of a superlattice of stripes,” Supercond. Sci. Technol. 25, 124002 (2012).
  * Chen et al. (2012) Y. Chen, A. A. Shanenko, A. Perali, and F. M. Peeters, “Superconducting nanofilms: molecule-like pairing induced by quantum confinement,” J. Phys.: Condens. Matter 24, 185701 (2012).
  * Salasnich et al. (2019) L. Salasnich, A. A. Shanenko, A. A. Vagov, J. A. Aguiar, and A. Perali, “Screening of pair fluctuations in superconductors with coupled shallow and deep bands: A route to higher-temperature superconductivity,” Phys. Rev. B 100, 064510 (2019).
  * Agrestini et al. (2004) S. Agrestini, C. Metallo, M. Filippi, L. Simonelli, G. Campi, C. Sanipoli, E. Liarokapis, S. De Negri, M. Giovannini, A. Saccone, et al., “Substitution of Sc for Mg in MgB2: Effects on transition temperature and Kohn anomaly,” Phys. Rev. B 70, 134514 (2004).
  * Bianconi et al. (2014) A. Bianconi, D. Innocenti, A. Valletta, and A. Perali, “Shape resonances in superconducting gaps in a 2DEG at oxide-oxide interface,” J. Phys. Conf. Ser. 529, 012007 (2014).
  * Mazziotti et al. (2021b) M. V. Mazziotti, R. Raimondi, A. Valletta, G. Campi, and A. Bianconi, “Resonant multi-gap superconductivity at room temperature near a Lifshitz topological transition in sulfur hydrides,” J. Appl. Phys. 130, 173904 (2021b).
  * Mazziotti et al. (2021c) M. V. Mazziotti, T. Jarlborg, A. Bianconi, and A. Valletta, “Room temperature superconductivity dome at a Fano resonance in superlattices of wires,” EPL(Europhysics Letters) 134, 17001 (2021c).
  * Tromp et al. (2022) W. O. Tromp, T. Benschop, J.-F. Ge, I. Battist, K. M. Bastiaans, D. Chatzopoulos, A. Vervloet, S. Smit, E. van Heumen, M. S. Golden, et al., “Puddle formation, persistent gaps, and non-mean-field breakdown of superconductivity in overdoped (Pb,Bi)2Sr2CuO6+δ,” arXiv:2205.09740 (2022), eprint https://arxiv.org/abs/2205.09740.