\section{Introduction}

The performance of machine learning algorithms crucially depends on the representation of data they work with. Relevant for the particular task data representation can significantly improve the performance of machine learning models on this task. As it was shown by Schmidhuber el al.~in \cite{schmidhuber}, Bengio et al.~in \cite{bengio} and I. Higgins et al.~in \cite{higgins}, disentangled (or statistically independent) representations are particularly valuable, because they are useful in a wide range of applied tasks, particularly in the computer vision field.

Advanced supervised machine learning models are known to be able to learn high-abstract relevant non-redundant representations, but they require specific preprocessing and a large amount of labeled data to achieve this. However, labeled data are usually limited and expensive to obtain. That is why the problem of learning good data representation from unlabeled data is important.

The requirements to amount of data required for learning good representation can be softened by using data augmentation technique. Initial dataset can be expanded by applying set of transformations to each training sample, where each of these transformation is known a priori to not affect desired output. However, when using data augmentation, we still need decent initial amount of labeled data, and only labeled data can be used.

A frequent situation in practice is when large amount of unlabeled samples is available, but labeling them is expensive. While data augmentation helps to significantly improve results in supervised tasks and learn good data representation during training, it usually does not help to learn more meaningful representation when applied naively in unsupervised tasks.

To see the problem, consider an example of CNN autoencoder for digit images, and assume we want it to extract an interpretable compact image representation, such that we can recognize digit from it easily (by linear classifier) regardless of rotation, scale, etc. If we provide our model an image of \textbf{3} and the same image, rotated $45^\circ$, it is really hard to infer relation between these images from training signal (pixel-wise reconstruction loss). In fact, in terms of pixel $L_2$ distance (typically used as a reconstruction loss) \textbf{3} is closer to \textbf{9} than to $45^\circ$ rotated \textbf{3} (contrary to classification task, where the rotation invariance is implicitly encoded in the training signal). Because of this, applying rotation augmentations will help the autoencoder to memorize all rotated variations of digits, but it still fails to generalize image classes over rotations. This case was experimentally investigated in \cite{kuzminykh}. This kind of generalization we want to obtain can be equivalently formulated as disentanglement of class-related part of image description from rotation related part in terms of linear subspaces, i.e. factoring by rotation group.

More generally, when we work with natural images, we want to separate image description space into two subspaces of features: the feature subspace, representing inherent properties of the depicted objects themselves and the feature subspace, representing different kinds of transformations (zooms, rotations, shifts etc).

One particular approach to such a separation was proposed by D. Kuzminykh et al.~in \cite{kuzminykh}.
By using properties of the Group Equivariant Convolutional Networks \cite{kuzminykh}, the autoencoder with specialized architecture can factor it's internal representation space by group of transformations used by group equivariant convolutions in encoder network, allowing to distill all invariant information. It is shown that such representation can successfully disentangle image classes, making them linearly separable once irrelevant transformation groups are factored out. However, that approach has number of disadvantages, arising from reliance on group equivariant convolutions: group equivariant convolution operation needs to be implemented manually for desired group, which is not always possible or require additional tricks, set of transformations need to be proper algebraic group with known structure, computation and memory costs grow with the group size.

In this paper, we introduce new, more flexible approach, inspired by observations on how biological brains do unsupervised learning. Instead of using specialized architecture to factor particular group of transformations, we change learning process to train encoder to produce more predictable image descriptions. To do so, we go away from training autoencoder on i.i.d. samples, and instead train it on conditionally connected samples, represented as tuples (image, transformation, result).

\subsection{Biological Inspiration}

Biological neural networks are known to learn good data representation in an unsupervised manner. We suppose that one of the main reasons for it is that brain always learns in context. It remembers what it has seen before and tries to predict what will happen next. Particularly, a brain of an animal always knows where the animal is moving at every moment due to the signals it exchanges with muscles.
Incredible importance of such a context was shown in paper \cite{held}. Authors described a biological experiment, which was performed on pairs of newborn kittens at early stage of their post-natal development, which is considered as crucial for their learning. Each pair was placed into an apparatus, shown on figure 1. One of the kittens (active subject, $A$) had relative freedom to move actively around axes $a$, $b$ and $c$. Another one (passive subject, $P$) was fixed inside of gondola, which repeated movements of $A$ automatically, due to mechanical transmission system. Passive subject couldn't influence its own movements. In the same time, both subjects could equally observe their environment.

\begin{figure}[htp]
    \centering
    \includegraphics[width=10cm]{Cats.png}
    \caption{Apparatus for equating motion and consequent visual feedback for an actively moving (A)
and a passive passively moved (P) S. Reprinted from \cite{held}, p.873}
    \label{fig:cats}
\end{figure}

After spending suitable amount of time inside this apparatus, kittens were released, and their responses on visual stimuli were tested. Passive kittens didn't manage to react properly on approaching objects and obstacles, as if they didn't see them properly, despite the fact, that their eyes were healthy. In the same time, active kittens didn't have such a problem. They develop normal response on visual environment. Active kittens were able to learn relationship between their own movements and the visual data they received, and this context turned out to be crucial for them to learn, how to interpret visual data in a meaningful way.

But unsupervised machine learning models usually don't have such a context. Their sources of information are limited to distinct data samples from a training set via the loss function. It doesn't have context that a biological brain has. Autoencoder models ``know'' that there is some compressed data representation and that they ``should'' find it, but found compressed representation isn't necessarily meaningful. Moreover, the task of compressing data can be performed fully without understanding the data structure. For example, JPEG algorithm is able to compress images, but is fully unable to extract useful features from it.

Having in mind the biological analogy, we think, that providing an additional context to autoencoder would let it learn more abstract and meaningful features. Particularly, we want to teach it the link between transformation (movement) and a new picture, which it gets after the movement.

\section{Related Work}

There are several different approaches to the described problem, known so far. Two especially influential of them are $\beta$-VAE \cite{held} and learning factorial codes by predictability minimization \cite{higgins}. Among more recent approaches, we want to highlight Multi-Modal Deep Clustering \cite{mmdc}, Maximizing Mutual Information Across Views \cite{mmiav} and Contrastive Learning \cite{chen}. Let us briefly highlight the insights of these works.

\subsection{Factorial Codes}

Factorial Codes are a particular variation of the lossless encoding. A factorial code of a given data piece is a vector of statistically independent features, which can be decoded into that original piece of data again. In other words, that vector should get rid from statistical dependencies inside the data piece, but save all the information, which it contained, in the same time.

A particular way for unsupervised learning of such a non-redundant representation was introduced in \cite{higgins}. There authors consider an abstract learner, which learns a representation of the data, using some hidden units. Then they assign an adaptive predictor to each hidden (representation) unit of this learner. Each predictor learn to predict the output of the corresponding unit, using the outputs of all other units of the same representational layer. In turn, each representational unit try to minimize its own predictability. This competition encourage each representational unit to learn concepts, statistically independent of those upon which the other units focused.

Authors also demonstrate a particular implementation of this idea, performing experiments with a particular neural network architecture and a learning algorithm for it.

\subsection{\texorpdfstring{$\beta$}--VAE}

Authors of $\beta$-VAE adapted well-known variational autoencoder \cite{VAE} model for learning factorized interpretable data representations. To achieve this, they specify disentangled metric and use it to enhance VAE's loss function in a way so to make it consider an ``entanglement'' as a part of the loss.

This enhanced loss function includes a hyperparameter $\beta$, which balances latent channel capacity and independence constraints with reconstruction accuracy (and which gives the model it's name). The edge case of $\beta = 1$ gives the classical VAE model. The further $\beta$ parameter is from 1, the more the ``entanglement'' part of the loss matters.

\subsection{Multi-Modal Deep Clustering}

In \cite{mmdc} a new algorithm for training a CNN-based neural network for natural images clustering is proposed. This algorithm combines main, unsupervised task of "noise as target" encoding with auxiliary task of rotation prediction, in order to obtain a representation with meaningful, clearly defined clusters.

"Noise as target" method aligns embeddings of images with target points sampled from a Gaussian Mixture. For this purpose it utilizes loss function based on the distance between image embedding and a closest target point. The additional task of rotation prediction encourages rotations of the same image to map into close vectors in the representation layer.

\subsection{Maximizing Mutual Information Across Views}

Authors of \cite{mmiav} encourage their model to maximize mutual information between features extracted from independently-augmented copies of the same image. The purpose of using such augmentations in this particular work is to create multiple views of a shared context, and then to force a neural network to extract the most informative features, shared by all these views.

\subsection{Simple Framework for Contrastive Learning}

Authors of SimCLR \cite{chen} also aim to extract transformation-independent representations of visual data. They achieve this goal by combining original architecture, composed of encoder and an auxiliary network, with a suitable dataset augmentation. It's easy to see that our approaches share much in common; however, they have deep differences as well. We will point out these differences below, right after giving brief outline of the Contrastive Learning framework itself.

In SimCLR, convolutional enconder takes two different transformations of the same image as independent inputs and learns representation of each one inside its latent space. In turn auxiliary fully-connected network called the ``projection head'' learns to map (project) representations corresponding to transformations of the same images into vectors, close in terms of the cosine similarity. Both networks train together by minimizing normalized temperature-scaled cross entropy loss through gradient descent (see the original paper \cite{chen} for more details).

In our work auxiliary network plays another role: it learns transformations instead of projections. Besides, \cite{chen} operates with representations, which drop information about transformations, while our model aims to save information about transformations in a disentangled manner. It is achieved by basing our model on an autoencoder instead of just encoder.

\bigskip

All these results have shown the way in which the resulting data representation is useful in a variety of tasks, especially in the computer vision field. However, first two papers attempted to approach the problem as solely a part of a compression task. Recognizing the huge importance of their results, we nevertheless admit that those techniques are not sufficient to learn good representation of more complicated datasets. In turn, the last three papers revolve around extraction of transformation invariants, but without keeping information about transformations themselves in the latent layer. So, despite the fact that such representations are proved to be useful, they are not decodable.

Our model, on the other hand, uses these invariants for disentangling learned representation, keeping it decodable. To achieve this, we augment the underlying VAE model in such a way, that it not only solves the compression task, but also tries to predict certain objects transformations. We see it as a part of an important learning context discussed above, which $\beta$-VAE and factorial codes models lack by themselves. We also provide formal mathematical description of our idea, using the Category Theory language as well as the Bayesian inference, justifying our architectural solutions.

\section{Model Description}

\subsection{Categorical Background}

\begin{figure}
    \centering
    \nodefontsize

    \begin{tikzpicture}
        \foreach \s in {1,...,\n}
        {
           \node[draw, rectangle, rotate=(\s*360/\n)] at ({1.5*360/\n * (\s - 1)}:\radius - \s*1.5 mm + 5mm) {$3$};
        }
        \node at (0, 0) (Set) {$Set$};
    \end{tikzpicture}
    \hspace{\setspace}
    \begin{tikzpicture}
        \foreach \s in {1,...,\n}
        {
           \node[draw, rectangle, rotate=(\s*360/\n)] at ({360/\n * (\s - 1)}:\radius) {$3$};
           \draw[->, >=latex] ({360/\n * (\s - 1)+\margin}:\radius)
             arc ({360/\n * (\s - 1)+\margin}:{360/\n * (\s)-\margin}:\radius);
        }
        \node at (0, 0) (Cat) {$Category$};
    \end{tikzpicture}

    \caption{Set~vs.~Category}
    \label{set_vs_category}
    \normalsize
\end{figure}

Since we want to take into account the relationships between samples in our dataset, we suggest describing our dataset in the terms of the Category Theory. To begin, we recall the base definition of a category.

\smallskip

\begin{definition}
    A category $\mathcal{C}$ is an aggregate of two classes --- the class of objects $\Ob(\mathcal{C})$ and the class of morphisms $\Hom(\mathcal{C})$ with the following properties:
    \begin{enumerate}
    \item Each pair of objects $a, b \in \Ob(\mathcal{C})$ corresponds to a set $H_\mathcal{C}(a, b) \subset \Hom(\mathcal{C})$ of morphisms.

    \item Each morphism $\alpha \in \Hom(\mathcal{C})$ belongs to one and only one set $H_\mathcal{C}(x, y)$, where $x$ and $y$ are some objects from $\Ob(\mathcal{C})$.

    \item There is a composition operation defined in the class $\Hom(\mathcal{C})$. The composition of $\alpha \in H_\mathcal{C}(a, b)$ and $\beta \in H_\mathcal{C}(b, c)$ gives a morphism $\beta\alpha \in H_\mathcal{C}(a, c)$. This operation is associative, i.e. $\gamma(\beta\alpha) = (\gamma\beta)\alpha$ for any three morphisms $\alpha \in H_\mathcal{C}(a, b),\, \beta \in H_\mathcal{C}(b, c),\, \gamma \in H_\mathcal{C}(c, d)$.

    \item Each set $H_\mathcal{C}(a, a)$ contains a special morphism $1_a$, such that for any morphisms $\alpha \in H_\mathcal{C}(x, a)$ and $\beta \in H_\mathcal{C}(a, y)$ $1_a \alpha = \alpha$ and $\beta 1_a = \beta$. This special morphism is called the identity morphism of object $a$.
    \end{enumerate}
\end{definition}

\smallskip

For more detailed explanation, see \cite{category}.

\smallskip

Returning to our situation, let's consider the set $\mathbf{X}$ of observations and the set $\mathbf{A}$ of allowed transformations between these observations. For example, $\mathbf{X}$ can be the set of some natural images, and $\mathbf{A}$ --- the set of possible rotations, shifts and rescalings of each object. Now, we introduce a category $\mathcal{X}$ in a following way:

\begin{itemize}
\item we use the objects from $\mathbf{X}$ as objects of $\mathcal{X}$: $\Ob(\mathcal{X}) = \mathbf{X}$;
\item the pair of objects $x_1$ and $x_2$ from $\Ob(\mathcal{X})$ is connected by a morphism $\alpha \in H_\mathcal{X}(x_1, x_2)$, if and only if there exists a transformation $a \in \mathbf{A}$, such that $x_2 = a(x_1)$;
\item the composition of morphisms in $\mathcal{X}$ is consistent with the composition of transformations from $\mathbf{A}$.
\end{itemize}

Our main object of interest is, however, the latent space $\mathbf{Z}$ of our model, since we suppose that each object from $\mathbf{X}$ is generated by some hidden variable $z \in \mathbf{Z}$.
Let's introduce a category $\mathcal{Z}$ in a following way:
\begin{itemize}
\item we use the objects from $\mathbf{Z}$ as objects of $\mathcal{Z}$;

\item we say, that there is a morphism $a': Z_1 \to Z_2$, if there is a morphism $a: X_1 \to X_2$. We denote this new set of morphisms $Hom(\mathcal{Z}) = \mathbf{A}'$.
\end{itemize}

Let's also denote the embedding, learned by our model, by $P : \mathbf{X} \to \mathbf{Z}$.

Now, let's consider the diagram on figure \ref{commutativity}. The goal of our model in the language of category theory is to make this diagram commutative. I.e.~for each $x_1 \in \mathbf{X}, a \in \mathbf{A}$, we want the existence of $a' \in \mathbf{A}'$, such that the equality $a'(P(x_1)) = P(a(x_1))$ holds. To achieve this we are using graph model, depicted on figure \ref{graph}.

\begin{figure}
    \centering

    \begin{tikzpicture}[node distance=2cm, auto]
    \node (X1) {$X$};
    \node (X2) [right of=X1] {$X$};
    \node (Z1) [below of=X1] {$Z$};
    \node (Z2) [below of=X2] {$Z$};

    \draw [decoration={markings,mark=at position 1 with
    {\arrow[scale=\arrowheadsizediag,>=stealth]{>}}},
    postaction={decorate}] (X1) -- node [above=0.5mm] {$a$} (X2);
    \draw [decoration={markings,mark=at position 1 with
    {\arrow[scale=\arrowheadsizediag,>=stealth]{>}}},
    postaction={decorate}] (Z1) -- node [above=0.5mm] {$a'$} (Z2);

    \draw [decoration={markings,mark=at position 1 with
    {\arrow[scale=\arrowheadsizediag,>=stealth]{>}}},
    postaction={decorate}] (X1) -- node [left=0.5mm] {$P$} (Z1);
    \draw [decoration={markings,mark=at position 1 with
    {\arrow[scale=\arrowheadsizediag,>=stealth]{>}}},
    postaction={decorate}] (X2) -- node [right=0.5mm] {$P$} (Z2);

    \draw [dashed, decoration={markings,mark=at position 1 with
    {\arrow[scale=\arrowheadsizediag,>=stealth]{>}}},
    postaction={decorate}] (X1) -- node [left=0.5mm] {} (Z2);
    \end{tikzpicture}

    \caption{Diagram we want to be commutative}
    \label{commutativity}
\end{figure}

\begin{figure}
    \centering
    \nodefontsize

    \begin{tikzpicture}[
    node distance=\nodedistance mm,
    nonterminal/.style=
    	circle,minimum size=\nodesize mm,
    	thick,draw=black,font=\ttfamily},
    terminal/.style=
    	circle,minimum size=\nodesize mm,
    	thick,draw=black,font=\ttfamily}]

    \node (Z1) [nonterminal] {$Z_1$};
    \node (X1) [nonterminal,below=of Z1] {$X_1$};
    \node (Z2) [nonterminal, right=of Z1] {$Z_2$};
    \node (X2) [nonterminal,below=of Z2] {$X_2$};

    \node (a)  [terminal,above=of Z1] {$a'$};

    \draw [decoration={markings,mark=at position 1 with
    {\arrow[scale=\arrowheadsize,>=stealth]{>}}},
    postaction={decorate}] node [above=0.8cm] {$p(z_1)$} (Z1) -- (X1) node [below=0.7cm] {$p(x_1|z_1)$} ;

    \draw [decoration={markings,mark=at position 1 with
    {\arrow[scale=\arrowheadsize,>=stealth]{>}}},
    postaction={decorate}] (Z2) -- (X2) node [below=0.7cm] {$p(x_2|z_2)$};

    \draw [decoration={markings,mark=at position 1 with
    {\arrow[scale=\arrowheadsize,>=stealth]{>}}},
    postaction={decorate}] (Z1) -- node [above=0.1cm] {} (Z2);

    \draw [decoration={markings,mark=at position 1 with
    {\arrow[scale=\arrowheadsize,>=stealth]{>}}},
    postaction={decorate}] (a) --  (Z2) node [above=0.8cm] {$p(z_2|z_1,a')$};

    \draw [dashed, decoration={markings,mark=at position 1 with
    {\arrow[scale=\arrowheadsize,>=stealth]{>}}},
    postaction={decorate}] (X1.west) -- (-\nodesize/10, -\nodesize/10-\nodedistance/10) -- node [left=0.1cm] {$q(z_1|x_1)$} (-\nodesize/10, 0) --  (Z1.west);

    \draw [dashed, decoration={markings,mark=at position 1 with
    {\arrow[scale=\arrowheadsize,>=stealth]{>}}},
    postaction={decorate}] (X2.east) -- (\nodedistance/10 + \nodesize * 2 / 10, -\nodesize/10 -\nodedistance/10) -- node [right=0.1cm] {$q(z_2|x_2)$} (\nodedistance/10 + \nodesize * 2 / 10, 0) -- (Z2.east);
    \end{tikzpicture}

    \caption{Graph model}
    \label{graph}
    \normalsize
\end{figure}

\subsection{Variational Inference}

Let us first remind the main steps of the variational inference for the classical Variational Autoencoder \cite{VAE} and then derive it for our modified model.

Recall the identity (see, for instance, \cite{tutorial_on_vae})
\begin{equation}
\label{elbo_identity}
     \mathbb{E}_{z\sim q(z|x)} [\log p(x|z)] - \KL{q(z|x)}{p(z)} = \mathcal{L}(x),
\end{equation}
which holds for any distribution $q(z|x)$ over the same domain $Z$ as $p(z|x)$. Here
\begin{equation}
\label{elbo_definition}
    \mathcal{L}(x) = \log p(x) - \KL{q(z|x)}{p(z|x)}
\end{equation}
is the evidence lower bound (ELBO) and
\begin{equation}
    \KL{p_1(z)}{p_2(z)}= \mathbb{E}_{z\sim p_1(z)} \log\frac{p_1(z)}{p_2(z)}
\end{equation}
is the Kullback–Leibler divergence between two distributions. Since the latter is always non-negative one can obtain the following lower bound of the evidence from definition (\ref{elbo_definition}):
\begin{equation}
    \log p(x) \geq \mathcal{L}(x).
\end{equation}
Moreover, according to (\ref{elbo_identity}), in order to maximize the data likelihood (with respect to some trainable parameters $\theta$ hidden in $p(x|z)$ and $q(z|x)$) one has to maximize objective function $\mathcal{L}(x)$ along with minimizing the KL-divergence between $q(z|x)$ and $p(z|x)$. As usual, we assume that our encoding model $q(z|x)$ is complex enough to nullify the divergence term without affecting the log-likelihood term, so it suffices just to maximize $\mathcal{L}(x)$.

In our case the space $X$ of all observable variables consists of triples $(x_1,x_2,a) \in \mathbf{X}\times\mathbf{X}\times\mathbf{A}$ of two images and a transformation. We want our generative process to produce triples for which $x_2=a(x_1)$ more likely than other, thus we want to maximize $p(x_1,x_2,a)$ for these points. The space $Z$ of all latent variables consists of triples $(z_1,z_2,a') \in \mathbf{Z}\times\mathbf{Z}\times\mathbf{A}$. Here $a'$ is a formal duplicate of observable variable $a$ in the latent space involved in the generative process. Recall that the objective of auxiliary encoding function $q(z|x)$ is to give us distribution over $z$ values that are likely to produce $x$, but its form is completely a matter of our choice. We use the following one:
\begin{equation}
\label{q_form}
    q(z_1,z_2,a'|x_1,x_2,a) = q(z_1|x_1) q(z_2|x_2) p(a'|a).
\end{equation}
The conditional distribution $p(a'|a)$ can be written in its explicit form $p(a'|a) = p(a|a') = \delta(a - a')$.

The equality that follows from our graph model depicted on figure \ref{graph} will also be useful:
\begin{equation}
\label{p_form}
    p(x_1,x_2,a|z_1,z_2,a') = p(x_1|z_1) p(x_2|z_2) p(a|a').
\end{equation}
Thus the whole generating process in our model is the following:
\begin{multline}
    p(x_1,x_2,a) = \mathbb{E}_{z_1,z_2,a'\sim p(z_1,z_2,a')} p(x_1,x_2,a|z_1,z_2,a')\\
    = \mathbb{E}_{z_1\sim p(z_1)} \mathbb{E}_{z_2\sim p(z_2|z_1,a)} p(x_1|z_1) p(x_2|z_2) p(a).
\end{multline}

In order to obtain the computable form of the evidence lower bound to train the model one has to substitute (\ref{q_form}) and (\ref{p_form}) into the identity (\ref{elbo_identity}) (written for specified above observable and hidden variables):
\begin{multline}
\label{th_objective}
    \mathcal{L}(x_1,x_2,a) = \mathbb{E}_{z_1,z_2,a'}[\log p(x_1,x_2,a|z_1,z_2,a')] - \KL{q(z_1,z_2,a'|x_1,x_2,a)}{p(z_1,z_2,a')}\\
    = \mathbb{E}_{z_1}[\log p(x_1|z_1)] +  \mathbb{E}_{z_2}[\log p(x_2|z_2)] +\log p(a) -\\
    - \KL{q(z_1|x_1)}{p(z_1)} - \mathbb{E}_{z_1}\KL{q(z_2|x_2)}{p(z_2|z_1,a)}.
\end{multline}
Here expectation $\mathbb{E}_{z_1,z_2,a'}$ is taken over distribution $q(z_1,z_2,a'|x_1,x_2,a)$ and expectations $\mathbb{E}_{z_i}$ over distributions $q(z_i|x_i)$ respectively.

In our experiments we did not use naturally augmented dataset $\bar{\mathcal{D}} = \{(x_1^i,x_2^i,a^i)|i\in I\}$, but did the augmentation by picking a transformation $a$ from finite predefined set $A = \{a^i|i\in I'\}$. Thus, given sampled $z_i\sim q(z_i|x_i)$ and $a\sim p(a)=\mathcal{U}_A(a)$ for every data point $x \in \mathcal{D}$, the objective (\ref{th_objective}) to maximize takes the form
\begin{multline}
\label{objective}
    \mathcal{L}(x,a(x),a|z_1,z_2) = \log p(x|z_1) + \log p(a(x)|z_2) + \log p(a) -\\
    - \KL{q(z_1|x)}{p(z_1)} - \KL{q(z_2|a(x))}{p(z_2|z_1,a)}.
\end{multline}
Here constant (with respect to trainable parameters) term $\log p(a)$ does not affect the maximization and can be omitted. However, one has to keep it if the distribution over possible transformations $p(a)$ is trainable itself.

\subsection{Details of Implementation}

More specifically, our model consists of two CNN encoders and two CNN decoders with the common parameter set and one full-connected network, which learns transformations (i.e.~morphisms $a'$ from $A'$) in their latent space $\mathbf{Z}$. All these networks are learning together via gradient descent with Adam optimizer, maximizing variational lower estimate of data likelihood (\ref{objective}).

In our particular implementation we used the following set  $A$ of allowed transformations:
\begin{itemize}
\item{up, down, left and right shifts by $3$px and $6$px;}
\item{scaling with coefficients $1.15$ and $1.32$;}
\item{rotations by $12^\circ$ and $24 ^\circ$.}
\end{itemize}

We assigned unique number for each transformation from this list. After it, at each training step we pick an object $x_1$ from our dataset with a particular transformation $a$ from $\mathbf{A}$ and obtain object $x_2$ by applying $a$ to $x_1$. Next, we compute the output of variational autoencoder on both objects --- $x_1$ and $x_2$.
Simultaneously we estimate a probability $p(z_2|z_1, a)$ of the object $z_2$ given $z_1$ and its desirable transformation $a$ with fully-connected network, mentioned above.

\section{Experiment Results}

In this section we compare the proposed model with several standard unsupervised algorithms that allow one to obtain a compressed representation of data in the latent space. We use the achievable accuracy of a linear classifier in the latent space as a measure of disentanglement of the obtained representation. We perform all the experiments on the MNIST dataset.

\subsection{2-dimensional Latent Space}

To demonstrate the distinctive features of the used algorithms we display the MNIST dataset embedded in 2-dimensional latent space by each of them (figures \ref{pca_2dim}--\ref{our_2dim}). To explore stability of such representation we also encode a sample of each of the ten handwritten digits from the dataset rotated $360^\circ$ with the step $3^\circ$. We call the obtained trajectories the orbits of these digits in the latent space.

\begin{figure}
    \centering
    \subfloat[The MNIST dataset in the latent space]{{
        \includegraphics[width=.45\textwidth]{experiments/2dim/pca_full_2dim.png}
    }}
    \qquad
    \subfloat[Orbit of each sample digit in the latent space]{{
        \includegraphics[width=.45\textwidth]{experiments/2dim/pca_orbits_2dim.png}
    }}
    \caption{PCA embedding (first two components)}
    \label{pca_2dim}
\end{figure}

\begin{figure}
    \centering
    \subfloat[The MNIST dataset in the latent space]{{
        \includegraphics[width=.45\textwidth]{experiments/2dim/umap_full_2dim.png}
    }}
    \qquad
    \subfloat[Orbit of each sample digit in the latent space]{{
        \includegraphics[width=.45\textwidth]{experiments/2dim/umap_orbits_2dim.png}
    }}
    \caption{UMAP embedding (with default hyperparameters, see \cite{umap})}
    \label{umap_2dim}
\end{figure}

\begin{figure}
    \centering
    \subfloat[The MNIST dataset in the latent space]{{
        \includegraphics[width=.45\textwidth]{experiments/2dim/vae_full_2dim.png}
    }}
    \qquad
    \subfloat[Orbit of each sample digit in the latent space]{{
        \includegraphics[width=.45\textwidth]{experiments/2dim/vae_orbits_2dim.png}
    }}
    \caption{VAE embedding}
    \label{vae_2dim}
\end{figure}

\begin{figure}
    \centering
    \subfloat[The MNIST dataset in the latent space]{{
        \includegraphics[width=.45\textwidth]{experiments/2dim/our_full_2dim.png}
    }}
    \qquad
    \subfloat[Orbit of each sample digit in the latent space]{{
        \includegraphics[width=.45\textwidth]{experiments/2dim/our_orbits_2dim.png}
    }}
    \caption{Our embedding (with hyperparameters $\beta = 0.01,\, \gamma = 5.0$)}
    \label{our_2dim}
\end{figure}

\begin{figure}[H]
    \centering
    \begin{tabular}{||c|c||}  \hline
         Algorithm & Linear separability \\  \hline \hline
         PCA & 44.7\% \\ \hline
         UMAP & 95.5\% \\ \hline
         VAE & 63.0\% \\ \hline
         Our & 65.7\% \\  \hline
    \end{tabular}
    \caption{Accuracy of a linear classifier on the 2-dimensional latent space of different algorithms}
    \label{results_2dim}
\end{figure}

One can see measured accuracy on the testing set of a linear classifier trained on the training set embedded in the latent representation by each algorithm in the figure \ref{results_2dim}. Our algorithm along with classical VAE is not well-suited for data visualization task, thus it is outperformed by UMAP, which is.

\subsection{64-dimensional Latent Space}

\begin{figure}
    \centering
    \subfloat[The MNIST dataset]{{
        \includegraphics[width=.45\textwidth]{experiments/64dim/vae_full_64dim.png}
    }}
    \qquad
    \subfloat[Orbit of each sample digit]{{
        \includegraphics[width=.45\textwidth]{experiments/64dim/vae_orbits_64dim.png}
    }}
    \caption{First two principal components projection of the 64-dimensional VAE embedding}
    \label{vae_64dim}
\end{figure}

\begin{figure}
    \centering
    \subfloat[The MNIST dataset]{{
        \includegraphics[width=.45\textwidth]{experiments/64dim/our_full_64dim.png}
    }}
    \qquad
    \subfloat[Orbit of each sample digit]{{
        \includegraphics[width=.45\textwidth]{experiments/64dim/our_orbits_64dim.png}
    }}
    \caption{First two principal components projection of our 64-dimensional embedding}
    \label{our_64dim}
\end{figure}

\begin{figure}[H]
    \centering
    \begin{tabular}{||c|c||}  \hline
         Algorithm & Linear separability \\  \hline \hline
         PCA & 91.7\% \\ \hline
         UMAP & 96.0\% \\ \hline
         VAE & 96.2\% \\ \hline
         Our & 98.4\% \\  \hline
    \end{tabular}
    \caption{Accuracy of a linear classifier on the 64-dimensional latent space of different algorithms}
    \label{results_64dim}
\end{figure}

Once the latent space is large enough, our algorithm outperforms all the baselines and allows us to obtain better visualization (by projecting onto the first two principal components in the latent space). The smoothness and stability of the orbits indicates that the latent code changes well-predictable under geometric transformations.

\section{Conclusion and Future Work}

In this paper we have described a particular approach that allows a model to learn transformation-aware representation of real world objects without using labeled data. We have also shown that our model outperforms other well-known models in high-dimensional latent space in some aspects. We have provided mathematical generalization of the proposed principle, using the language of the Category theory, and described Bayesian inference for our model.

We believe that our ideas can be developed further in several directions and can be useful in other applications:

\begin{itemize}
    \item Being generalized to the wider class of transformations (beyond Gaussian), our model can learn to perform stochastic multi-modal transformations properly. Now if new objects appear on the scene after applying particular transformation, our model tends to predict many inconsistent intermediate states that look like partially appeared objects. One of the possibilities here is to adapt Generative Adversarial Network for our purposes, or use variations of VAE-like models with more flexible distribution families.

    \item Learning geometry of observed world and unsupervised scene understanding. Trained on a video sequence labeled with observer's movements, our model will require deep understanding of the scene in order to predict the next frame of the video. Such labeled observation sequences may be automatically generated and extracted from 3D simulations. In this class of tasks multi-modal probability models described above can be especially useful.
\end{itemize}


\section{Introduction}

The discovery of Dirac and Weyl semimetals (DWSMs) has provided a
rich arena for probing novel gapless phases of matter with unique
transport properties and topological features\,\citep{Armitage2018}.
Several types of gapless systems featuring Dirac or Weyl points in
three-dimensional (3D) momentum space have been realized\,\citep{Borisenko2014,Han-Jin2017,Xu613}.
The simplest DWSMs exhibit twofold or fourfold degenerate linear-band
touching points at the Fermi level with isotropic velocities and a
possible replication into disjoint momentum-space valleys. Their pointlike
Fermi surface is protected against band gap opening due to either
topological constraints---in Weyl systems with broken time-reversal
({\small{}$\mathcal{T}$}) or inversion symmetries ({\small{}$\mathcal{P}$})\,\citep{Armitage2018}
--- or crystal symmetries in {\small{}$\mathcal{TP}$}-symmetric
Dirac systems\,\citep{Young2012,Yang2014}. Thus, any clean DWSM
is an incompressible electron gas with a quadratically vanishing density
of states (DoS). \textcolor{black}{Despite the inefficient charge
screening at the node, this paradigm is believed to survive electron-electron
Coulomb interactions, giving way to a marginal Fermi liquid behavior\,\citep{PhysRevB.90.121107,PhysRevB.92.045104}.
In addition, weakly interacting DWSMs can display strongly renormalized
Fermi velocities, but the nodes' integrity and topology are expected
to remain robust\,\citep{Interact_PhysRevB.94.241102,Interact_PhysRevB.97.161102,Interact_PhysRevB.98.241102,Interact_PhysRevLett.113.136402,Roy16,Roy17,Roy18,maciejko2014}. }

An outstanding question is whether random on-site potentials ubiquitous
in realistic systems (e.g., due to impurities in the crystal lattice)
can give way to a compressible diffusive metallic phase with a finite
nodal DoS\,\citep{Fradkin1986_2,Shindou2009,Goswami2011,Yang2014,Nandkishore2014,Syzranov2015,Sbierski2016,Pixley2016,Pixley2016_2,pixley2016_3,Syzranov2016,Buchhold2018,Buchhold2018_2,Ziegler2018,Goncalves2020,Wilson2020,Kobayashi2020}.
An early result by Fradkin predicted that Dirac nodes are stable in
{\small{}$d\!=\!2\!+\!\epsilon$} dimensions, below some critical
disorder strength\,\citep{Fradkin1986_2}. The robustness of DWSMs
against weak random perturbations is best visualized by considering
a massless particle moving through a short-ranged random potential
of strength{\small{} $W$}. Since, near a node, the de Broglie wavelength,
{\small{}$\lambda\!=\!\hbar v/E$}, largely exceeds the disorder correlation
length, the central limit theorem applies and the fluctuations around
the average potential inside a volume {\small{}$\lambda^{d}$ }must
scale as {\small{}$\delta V\!\!\propto\!W\!\lambda^{\!{\scriptscriptstyle -d/2}}$}.
In{\small{} $d\!=\!3$}, the fluctuations vanish as{\small{} $E^{{\scriptscriptstyle 3/2}}$,}
i.e., faster than the band energy near a node, rendering the semimetal
phase stable.

The early field-theoretical point of view has been recently questioned
by nonperturbative calculations\,\citep{Nandkishore2014, Pixley2016_2},
hinting that 3D gapless phases can become unstable due to the emergence
of\emph{ zero-energy states} bound to statistically rare regions of
the disorder potential landscape. According to this picture, the nodal
DoS remains nonzero for arbitrarily weak disorder without any signature
of singular behavior. Evidence for avoided quantum criticality (AQC)
facilitated by localized nodal eigenstates has been provided by lattice
simulations of a 3D Dirac model with uncorrelated on-site disorder\,\citep{Pixley2016_2}.
Challenging these findings, Buchhold \textit{et al.} noted that rare
events are preceded by scattering resonances which always carry zero
spectral weight at a node. Furthermore, as they are only possible
for fine-tuned (\textit{``magical}'') impurity configurations, this
would imply that the nodal DoS cannot be lifted at variance with the
AQC scenario\,\citep{Buchhold2018,Buchhold2018_2}. Their claim is
backed by the exact solution of a Weyl node with a spherical impurity.\textcolor{blue}{{}
}\textcolor{black}{These paradoxical findings have attracted significant
attention recently\,\citep{Pixley2016_2,roy_universal_2016,Ziegler2018,Wilson2020}
because they question the phase stability of incompressible 3D gapless
phases in realistic conditions. Moreover,}\textit{\textcolor{black}{{}
}}\textcolor{black}{the driving mechanism for semimetal-to-metal transitions
in the phase diagram of dirty 3D DWSMs remains elusive.}

\textcolor{black}{In this paper, we resolve this conundrum by tackling
the spherical impurity problem using two complementary theoretical
approaches. First, within a continuum model, we argue that an unforeseen
non-analytic behavior of scattering phase shifts at the node obstructs
a direct use of Friedel's sum rule (FSR)} \,\textcolor{black}{\citep{Buchhold2018,Buchhold2018_2}.
This difficulty can be overcome by keeping track of the level statistics
in systems of increasingly large volume at fixed impurity concentration.
The puzzling behavior of the phase shifts is explained by the emergence
and sudden disappearance of bound states}\textit{\textcolor{black}{\emph{
upon tuning the impurity potential across a fine-tuned}}}\textit{\textcolor{black}{{}
``magical value}}\textcolor{black}{''. Crucially, our level statistics
analysis reveals that rare bound states are accompanied by a }\textcolor{black}{\emph{continuum
of low-energy resonances}} surrounding the nod\textcolor{black}{e
in realistic material systems with a }\textcolor{black}{\emph{diversity}}\textcolor{black}{{}
of random} short-range impurities.\textcolor{red}{{} }Such near-critical
impurity configurations give effective statistical weight to \textit{\textcolor{black}{magical}}
impurities and ultimately endow the node with a finite average DoS.
Second, we carry out high-precision numerical calculations in a lattice
version of the problem, hosting one or more impurities. Remarkably,
our real-space simulations not only unambiguously demonstrate the
destabilization of a 3D DWSM by a diversity of random near-critical
impurities, but also quantitatively agree with the continuum theory
predictions\textcolor{black}{{} in the dilute impurity regime. These
findings provide the missing link between continuum and lattice approaches
to the DWSM theoretical puzzle and unambiguously pinpoint the newly
unveiled }\textit{\textcolor{black}{statistical significance of near-critical
impurities}}\textcolor{black}{{} as the driving mechanism for AQC. Hence,
dilute impurities can only destabilize a DWSM provided their random
parameters are drawn with a probability density which is nonzero on
(at least) one}\textit{\textcolor{black}{{} magical}}\textcolor{black}{{}
value.} Lastly, we note that subtleties in disordered Dirac systems
have a long history\,\citep{Atkinson2000,Hirschfeld2002, Pepin2001, Adagideli2002, Altland2002}.\,For
instance, in 2D $d$-wave superconductors, there are four low-energy
quasiparticle Dirac valleys, and scalar impurities are pair breaking.
The latter induce resonances that, in the strong scattering limit,
turn into sharp DoS peaks at $E\!\!=\!\!0$ (Majorana zero modes)\,\citep{Balatsky95, Adagideli99}.

\textcolor{black}{\vspace{-0.1cm}
 }

\textcolor{black}{The remainder of the paper is organized as follows:\,In
Sec.\,\ref{sec:Continuum-theory}, we present the theoretical tools
for calculating the DoS correction induced by dilute spherical scalar
impurities hosted within a single-node continuum model of a Dirac
semimetal. We further argue our conclusions to remain valid in 3D
Weyl semimetals. In Sec.\,\ref{sec:Impurity-induced-change-in},
we highlight the main caveats implied by a direct use of FSR and show
how to obtain consistently the thermodynamic limit DoS change due
to a finite (albeit small) concentration of impurities.\,In Sec.\,\ref{sec:Near-critical-impurities-lift}
this theory is used to predict the conditions on which AQC holds in
a Dirac semimetal. Our predictions are validated by high-resolution
real-space calculations in Sec.\,\ref{sec:Lattice-simulations}.
Finally, in Sec.\,\ref{sec:Outlook} we summarize our main findings
and highlight future directions for further study.}

\vspace{-0.3cm}

\section{\label{sec:Continuum-theory}Continuum theory}

\vspace{-0.2cm}

We start by reviewing the low-energy description o\textcolor{black}{f
a noninteracting single-valley 3D DWSM. The Hamiltonian can be written
as $\mathcal{H}_{0}\!\!=\!v\,\mathbf{\boldsymbol{\alpha}}\!\cdot\mathbf{p}$,
with $\hbar\!\equiv\!1$, $\alpha^{i}\!=\!\sigma^{x}\!\!\otimes\!\!\sigma^{i}$,
$v$ being the Fermi velocity, and $\mathbf{p}\!=\!-\imath\boldsymbol{\nabla}$
being the momentum operator. Here, $\sigma_{i}$ ($i\!\!=\!\!x,\!y,\!z$)
denote Pauli matrices acting on internal spin space. Introducing a
scalar impurity potential in the Hamiltonian breaks translation invariance,
but if $\mathcal{U}\!\left(\mathbf{r}\right)\!\!=\!\mathcal{U}\!\left(\abs{\mathbf{r}}\right)$,
rotational symmetry around the impurity center is preserved. For concreteness,
we consider a spherical well or plateau potential, $\mathcal{U}\!\left(\mathbf{r}\right)\!\!=\!\!\lambda\Theta\left(b\!-\!\left|\mathbf{r}\right|\right)$\,\citep{Nandkishore2014}.
We note that this model is suitable to describe realistic multiple-valley
Dirac or Weyl semimetals insofar as the impurity radius $b$ is much
larger than the lattice spacing (thus effectively suppressing intervalley
scattering). The coupling of distinct Weyl sectors at each valley
is also absent due to the scalar structure of the impurity potential.
The eigenstates of $\mathcal{H}\!=\!\mathcal{H}_{0}\!+\!\mathcal{U}\!\left(\abs{\mathbf{r}}\right)$
can be written as}

\begin{equation}
\Psi_{j,j_{z}}^{\kappa}\!\left(\mathbf{r}\right)\!=\!\left[\frac{f_{j}^{\kappa}\!\left(r\right)}{r}\Theta_{j,j_{z}}^{-\kappa}\!\left(\hat{\mathbf{r}}\right),\frac{\iota g_{j}^{\kappa}\!\left(r\right)}{r}\Theta_{j,j_{z}}^{\kappa}\!\!\left(\hat{\mathbf{r}}\right)\right]^{\text{T}},\label{eq:SpinorAngularStructure-1}
\end{equation}
where $j\in\{\nicefrac{1}{2},\nicefrac{3}{2},\cdots\}$ and $j_{z}\in\left\{ -j,-j+1,\cdots,j\right\} $
are the total angular momentum quantum numbers, while $\kappa=\pm1$
labels the eigenvalues of $\mathcal{K}\!=\!\gamma^{0}\!\cdot\!(2\mathbf{L}\cdot\mathbf{S}\!-\!1)$,
i.e., $\kappa(j\!+\!\nicefrac{1}{2})$. Furthermore, $\Theta_{j,j_{z}}^{\pm}\!\!\left(\hat{\mathbf{r}}\right)$
are orthonormal spin-$\nicefrac{1}{2}$ spherical harmonics and $f_{j}^{\kappa}\left(r\right)$/$g_{j}^{\kappa}\left(r\right)$
are radial functions. For nonzero energies, the latter are radial
spherical waves with phase shifts introduced by the central potential
(see Appendix\,\eqref{sec:AppendixA} for additional details). In
each $j$ sector, the scattering phase shifts $\delta_{j}$ induced
by a spherical well or plateau are obtained by constraining the spinor
to be continuous at $r\!=\!b$. One obtains, after a somewhat lengthy
calculation,
\begin{widetext}
\vspace{-0.4cm}

\begin{equation}
\tan\delta_{j}\left(\varepsilon,u\right)=\frac{\text{sgn}\left(\varepsilon\!-\!u\right)J_{j+1}\left(\abs{\varepsilon}\right)J_{j}\left(\abs{\varepsilon-u}\right)-\text{sgn}\left(\varepsilon\right)J_{j}\left(\abs{\varepsilon}\right)J_{j+1}\left(\abs{\varepsilon-u}\right)}{\text{sgn}\left(\varepsilon\right)\text{sgn}\left(\varepsilon-u\right)Y_{j+1}\left(\abs{\varepsilon}\right)J_{j}\left(\abs{\varepsilon-u}\right)-Y_{j}\left(\abs{\varepsilon}\right)J_{j+1}\left(\abs{\varepsilon\!-\!u}\right)},\label{eq:B/A-1-1-1}
\end{equation}

\vspace{-0.35cm}

\end{widetext}

\noindent where $u\!\equiv\!\lambda b/v_{\text{F}}$, $\varepsilon\!\equiv\!Eb/v_{\text{F}}\!\neq\!(0,u)$,
and $J_{n}(x)\,[Y_{n}(x)]$ are Bessel functions of the first (second)
kind. We underline that Eq.\,(\ref{eq:B/A-1-1-1}) is equivalent
to that obtained in Refs.\,\citep{Nandkishore2014,Buchhold2018_2}
for the Weyl equation. This is unsurprising because our $4\times4$
Dirac model is gapless and the impurity potential has a scalar structure
(i.e. distinct Weyl sectors at each valley remain decoupled). However,
Eq.\,\eqref{eq:B/A-1-1-1} only defines $\delta_{j}\!\left(\varepsilon,u\right)$\,mod$\pi$.
The ambiguity corresponds, at most, to a global change in the sign
of the wavefunction. In order to obtain a unique definition of $\delta_{j}\left(\varepsilon,u\right)$,
one needs to choose a reference point, i.e., as the potential is switched
off ($u\!\to\!0$), the require that the phase shifts vanish across
the entire spectrum\textcolor{black}{. A way to guarantee this is
to enforce that $\delta_{j}\left(\varepsilon\!\to\!\pm\infty,u\right)\!=\!-u$\,\citep{Ma85,Ma_2006}\,}\footnote{\textcolor{black}{As $\varepsilon\!\rightarrow\!\pm\infty$, one has
$\tan\delta_{j}\left(\varepsilon,u\right)\!\rightarrow\!-\tan u$.
This differs from the limiting behavior introduced in Ref.~\citep{Ma85},
$\lim\!{}_{\varepsilon\to\pm\infty}\!\delta_{j}\left(\varepsilon,u\right)\!\!=\!\!\mp u$.
The difference comes from the $\text{sgn}(\varepsilon)$ factors introduced
when defining the radial scattering wavefunctions {[}see Eqs.\,\eqref{eq:g+FreeSolution-3-2}---\eqref{eq:g-FreeSolutionOutside}
of Appendix\,\ref{sec:AppendixA}{]}}}\textcolor{black}{, which is achieved by a trick explained in Appendix\,\ref{sec:AppendixC}.
Since the analysis is qualitatively similar in all $j$ sectors, in
what follows we focus on the $\delta_{1/2}\!\left(\varepsilon\right)$
phase shift (see Fig.\,\ref{fig:PhaseShift - -Pi to Pi}) using the
previous convention (for completeness, plots for other $j$'s and
around different magical $u$'s are provided in Appendix\,\ref{sec:AppendixC}).
For $u\!\!=\!\!u_{c}^{j}$, the phase shift is shown to have a physical
$\pi$ discontinuity at $\varepsilon\!\!=\!\!0$, which marks the
occurrence of zero-energy bound states}\,\citep{Ma_2006}\textcolor{black}{.
For} the 3D massless Dirac equation, bound states at $\varepsilon\!=\!0$
can appear, for particular wells or plateaus, whenever a decoupling
of the radial equations for $f_{j}^{\pm}\left(r\!>\!b\right)/g^{\pm}\left(r\!>\!b\right)$
occurs\,\citep{Nandkishore2014}. In this case, the admissible (asymptotically
decreasing) solutions are simple power laws, $g_{j,k}^{+}\left(r\!>\!b\right)/f_{j,k}^{-}\left(r\!>\!b\right)\!=\!\mathcal{B}^{\pm}r^{-j-1/2}$
and $f_{j,k}^{+}\left(r\!>\!b\right)/g_{j,k}^{-}\left(r\!>\!b\right)\!\!=\!\!0$.
As shown in Appendix\,\eqref{sec:AppendixA}, such spinors are only
continuous at $r\!=\!b$, if the potential satisfies $J_{j}\left(\abs u\right)\!=\!0$.
Hence zero-energy states are allowed in a single-impurity Dirac problem
provided the parameters are fine-tuned, i.e., $\abs{\lambda b/v}\!=\!u_{c}^{j}$
is a root of $J_{j}\left(x\right)$. The critical parameters $\{u_{c}^{j}\}$
are dubbed \emph{magical values}, as they would correspond to rare
regions in a disordered landscape where nonperturbative zero-energy
modes are possible\,\citep{Nandkishore2014,Pixley2016_2}. Note that
these are true (squared-normalizable) impurity bound states within
the Dirac continuum and they have a degeneracy of $2\!\left(\!2j\!+\!1\!\right)$.
These bound-states manifest themselves as a $\pi$ discontinuity at
$\varepsilon\!=\!0$ in the phase shifts when the parameter $u$ crosses
a critical value of that angular momentum channel (inset of Fig.\,\ref{fig:PhaseShift - -Pi to Pi}).
This is in accord with Levinson's theorem for Dirac particles\,\citep{Ma85,Ma_2006}
which states that given an appropriate convention, the number of bound
states is encapsulated in discontinuous $\pi$ jumps of the phase
shifts at zero momentum.

\begin{figure}[t]
\vspace{-0.2cm}

\begin{centering}
\hspace{-0.1cm}\includegraphics[scale=0.165]{Fig1.pdf}
\par\end{centering}
\vspace{-0.35cm}

\caption{\label{fig:PhaseShift - -Pi to Pi}Plot of the energy-dependent phase
shift $\delta_{{\scriptscriptstyle 1/2}}\!\left(\varepsilon,u\right)$
in accordance with the prescription $\delta_{j}\!\left(\varepsilon\!\to\!\pm\infty,u\right)\!=\!-u$.
Several values of $u$ are plotted around $u_{c}=\pi$\textcolor{black}{.\,The
inset depicts a close-up of the curves near $\varepsilon\!=\!0$.}\vspace{-0.8cm}
 }
\end{figure}

\vspace{-0.7cm}

\section{\label{sec:Impurity-induced-change-in}Impurity-induced change in
the DoS\textup{ }}

\vspace{-0.2cm}

\noindent The change in the DoS induced by an isolated impurity is
conventionally calculated using FSR,

\vspace{-0.55cm}

\begin{equation}
\Delta\nu(\varepsilon,u)=\frac{2}{\pi}\sum_{j=1/2}^{\infty}\left(2j+1\right)\frac{\partial\delta_{j}(\varepsilon,u)}{\partial\varepsilon},\label{eq:FSR}
\end{equation}

\textcolor{black}{\vspace{-0.2cm}
 }

\noindent \textcolor{black}{which measures the variation in the number
of states per unit energy (an extensive quantity, not to be confused
with $\Delta\nu$ per unit volume hereafter denoted by $\Delta\rho$).
Strik}ingly, the phase-shift discontinuity caused by the impurity
bound state in the 3D DWSM problem precludes the direct use of FSR,
a peculiar effect that has gone unnoticed in earlier studies. Therefore,
to determine $\Delta\nu$, we adopt a strategy based on counting states
within a finite energy window adapted from Friedel's original reasoning\,\citep{Friedel1952}.
First, we restrict the Dirac fermions to lie inside a finite sphere
of radius $R$, $S_{R}$. The Hermiticity of $\mathcal{H}$ is guaranteed
if the Hilbert space is restricted to states with vanishing current
across the spherical surface, $\partial S_{R}$, that is, to a subspace
where any two spinors $\Psi$ and $\text{\ensuremath{\Phi}}$, satisfy
$\varoiint_{\partial S_{R}}\!d\mathbf{S}\!\cdot\!\left[\!\Psi_{\mu}^{\dagger}\left(\mathbf{r}\right)\!\boldsymbol{\alpha}_{\mu\nu}\Phi_{\nu}\left(\mathbf{r}\right)\!\right]\!=\!0.$
In Eq.\,\eqref{eq:SpinorAngularStructure-1}, this is true if $\cos\!\delta_{j}\left(\varepsilon,u\right)J_{j}\left(\abs{\varepsilon}\!R\right)\!-\!\sin\!\delta_{j}\left(\varepsilon,u\right)Y_{j}\left(\abs{\varepsilon}\!R\right)\!=\!0.$
For each angular momentum sector, this condition quantizes the allowed
energy levels, which we denote by $\varepsilon_{n}^{j}$. The number
of ($j$-sector) levels inside the energy window $[\varepsilon_{0}-\Delta\varepsilon/2,\varepsilon_{0}+\Delta\varepsilon/2]$
is changed by the impurity due to an inwards or outwards migration
of levels from regions of width $\simeq\delta_{j}\left(\varepsilon\pm\nicefrac{\Delta\varepsilon}{2},u\right)/R$
{[}up to $\mathcal{O}\left(R^{-2}\right)${]} near the respective
boundaries. This mechanism is illustrated in Fig. \ref{fig:LowEnergyDoS_Altland}(a).
The variation in the number of $j$ states inside the probing window
is

\vspace{-0.45cm}

\begin{align}
\Delta N_{j}\!\left({\scriptstyle \varepsilon_{0},\Delta\varepsilon,u}\!\right) & \!=\!\frac{2\left(2j\!+\!1\right)}{\pi}\!\left[\delta_{j}\!\left({\scriptstyle \varepsilon_{0}+\frac{\Delta\varepsilon}{2},u}\right)\!-\!\delta_{j}\left({\scriptstyle \varepsilon_{0}-\frac{\Delta\varepsilon}{2},u}\right)\right]\!.\label{eq:Number of States Variation}
\end{align}

\vspace{-0.2cm}

\noindent For a finite $\Delta\varepsilon$, Eq.\,\eqref{eq:Number of States Variation}
is accurate in the asymptotic limit $R\gg1$ and for $\varepsilon_{0}\pm\nicefrac{\Delta\varepsilon}{2}\neq0$
as explained and illustrated in Appendix\,\eqref{sec:AppendixB}.

Next, we consider the intensi\textcolor{black}{ve DoS ($\Delta\rho$)
indu}ced by a finite\textcolor{black}{{} density} ($c$) of impurities
in a volume $V$ focusing on \textcolor{black}{single scattering events.
The neglect of quantum-coherent scattering by multiple impurities
is justified in the dilute regime, where quantum interference corrections
are suppressed by a factor of $1/(k_{F}l)$\,\citep{altland_condensed_2010},
where $l\!\propto\!c^{-1}$ is the mean free path and $k_{F}=E/v\hbar$
is the Fermi wavevector. This is an important assumption confirmed
precisely by our numerical simulations below.}\textcolor{red}{{} }Formally,
the DoS is obtained by the limiting procedure, $\Delta\rho\left(\varepsilon_{0}\right)\!=\!\lim_{\Delta\varepsilon\to0^{+}}\lim_{V\rightarrow\infty}\left[\Delta N\left(\varepsilon_{0},\Delta\varepsilon,\left\{ u_{i}\right\} ,V\right)/V\Delta\varepsilon\right],$
where $i$ indexes the impurity and $\Delta N\left(\varepsilon_{0},\Delta\varepsilon,\left\{ u_{i}\right\} ,V\right)$
is the variation in the total number of states. Assuming that $\left\{ u_{i}\right\} $
are drawn from a probability density function $p\left(u\right)$,
the thermodynamic limit then reads

\vspace{-0.45cm}

\begin{align}
\sum_{i}\frac{\Delta N_{j}(\varepsilon_{0},\Delta\varepsilon,u_{i})}{V}\!\!\underset{{\scriptscriptstyle V\to\infty}}{\longrightarrow}\!c\!\int\!du\,p(u)\Delta N_{j}(\varepsilon_{0},\Delta\varepsilon,u),\label{eq:ThermodynamicLimitDiversity}
\end{align}
\vspace{-0.45cm}

\noindent where $\Delta N_{j}\!\left(\varepsilon_{0},\Delta\varepsilon,u\right)$
is given by Eq.~(\ref{eq:Number of States Variation}). The final
expression for the DoS variation due to a dilute set of random impurities,
$\Delta\rho\left(\varepsilon_{0}\right)\!=\!\!\sum_{j}\Delta\rho_{j}\left(\varepsilon_{0}\right)$,
is obtained from

\vspace{-0.45cm}

\begin{multline}
\Delta\rho_{j}(\varepsilon_{0})\!=\!c\!\lim_{{\scriptscriptstyle \Delta\varepsilon\to0^{+}}}\!\!\int du\,p\left(u\right)\frac{\Delta N_{j}\left(\varepsilon_{0},\Delta\varepsilon,u\right)}{\Delta\varepsilon}.\label{eq:Number of States Variation-2}
\end{multline}

\vspace{-0.3cm}

\noindent

\noindent
\begin{figure}[t]
\vspace{-0.3cm}

\begin{centering}
\hspace{-0.25cm}\includegraphics[scale=0.235]{Fig2.pdf}
\par\end{centering}
\vspace{-0.35cm}

\caption{\label{fig:LowEnergyDoS_Altland}(a) Motion of energy levels triggered
by the central impurity. (b) Plots of $\Delta\rho_{{\scriptscriptstyle 1/2}}\!\left(\varepsilon,u\right)$
for selected values of $u$ around $u_{\text{c}}\!=\!\pi$. The curves
were obtained using Eq.(\ref{eq:FSR}) with the phase shifts represented
in Fig.\ref{fig:PhaseShift - -Pi to Pi}. As one approaches $u_{\text{c}}$,
$\Delta\rho_{1/2}\!\left(\varepsilon\right)$ takes the form of a
low-lying peak, which gets narrower and closer to $\varepsilon\!=\!0$
(roughly conserving the area between zeros). The inset shows that
$\Delta\rho_{1/2}\!\left(\varepsilon,u\right)$ is always equal to
zero at $\varepsilon\!=\!0$ for any non-critical $u$. (c) Depiction
of $\Delta N_{1/2}\!\left(\varepsilon\!=\!0,\pi\right)/\Delta\varepsilon$
converging towards a distribution $4\delta\!\left(u\!-\!\pi\right)$
as $\Delta\varepsilon\!\to\!0^{+}$. (d) Theoretical prediction for
$\Delta\rho_{1/2}\!\left(\varepsilon\right)$ due to dilute spherical
impurities with a Gaussian diversity of width $\sigma$ around $u_{c}^{{\scriptscriptstyle \!1\!/\!2}}\!\!=\!\pi$.
The inset depicts the total DoS for $10^{-6}b^{-1}$ impurities per
volume. \vspace{-0.8cm}
}
\end{figure}

\noindent The order of limits in Eqs. (\ref{eq:ThermodynamicLimitDiversity})
and (\ref{eq:Number of States Variation-2}) is essential. The integration
over $u$ must be done prior to taking the $\Delta\varepsilon\!\to\!0^{{\scriptscriptstyle +}}$
limit\textcolor{black}{. This is reminiscent of lattice simulations,
where the resolution parameter must be sent to zero only after the
thermodynamic limit has been taken (see below). }If $\delta_{j}\!\left(\varepsilon,u\right)$
is differentiable at $\varepsilon\!=\!\varepsilon_{0}$, the $\Delta\varepsilon\!\to\!0^{{\scriptscriptstyle +}}$
limit can be safely brought inside the integral, and one obtains $\Delta\rho_{j}\left(\varepsilon_{0}\right)\!=\!c\left(4j\!+\!2\right)\!\av{\partial\delta_{j}\left(\varepsilon,u\right)\!/\!\partial\varepsilon|_{\varepsilon=\varepsilon_{0}}\!}_{\!u}\!/\pi,$
where $\langle f\rangle_{u}\!=\!\int\!du\,p\left(u\right)\!f(u)$,
i.e., the familiar FSR. A direct application of FSR\textcolor{black}{\,{[}Eq.\,\eqref{eq:FSR}{]}}
was employed in Refs.\,\citep{Buchhold2018,Buchhold2018_2} to determine
the DoS in Weyl systems with statistical fluctuations of $u$ around
a critical value $u_{c}$, leading to the conclusion that $\Delta\rho_{j}\!\!\left(0\right)\!=\!0$.
This was inferred from the fact that $\partial\delta_{j}\left(\varepsilon,u\right)/\partial\varepsilon|_{\varepsilon=0}\!=\!0$
for any $u\!\neq\!u_{c}^{j}$; see Fig.\,\ref{fig:LowEnergyDoS_Altland}(b).
\textcolor{black}{Since critical configurations ($u\!=\!u_{c}^{j})$
have zero statistical measure, FSR would seemingly imply a vanishing
average DoS at $\varepsilon\!=\!0$. In the remainder of this paper,
we show that the difficulty arising from the discontinuous $\delta_{j}\!\left(\!\varepsilon,\!u\!=\!u_{\text{c}}\!\right)$
can be overcome by carefully accounting for the level statistics in
the infinite volume limit of a DWSM with random impurities. Such a
procedure does not alter the fate of $\rho\left(\varepsilon\!=\!0\right)$
in the presence of a finite concentration of impurities with a fixed
$u$\,\citep{Buchhold2018_2,Nandkishore2014}. However, the conclusions
on the stability of a DWSM are changed dramatically if a }\textit{\textcolor{black}{continuous
statistical diversity of ``near-critical'' impurities}}\textcolor{black}{{}
exists. }

\vspace{-0.6cm}

\section{\label{sec:Near-critical-impurities-lift}Near-critical impurities
lift the nodal density of states}

\vspace{-0.2cm}

\noindent \textcolor{black}{A finite concentration of exactly critical
wells or plateaus would introduce a macroscopic number of nodal bound
states. However, for a diversity of random impurities with potential
strengths drawn from a probability distribution function $p(u)$,
such fine-tuned configurations appear with zero probability and cannot
yield statistically significant contributions to the bulk nodal DoS\,\citep{Buchhold2018,Buchhold2018_2}.
Nevertheless, we find that low-energy resonances due to }\textit{\textcolor{black}{near-critical}}\textcolor{black}{{}
configurations ($u\!\approx\!u_{c}^{j}$) provide such a contribution.~The
phase shifts of such impurities signal the emergence of the zero-energy
bound states by a sharp resonance, namely, a quick $\pi$ variation
of $\delta_{j}\!\left(\varepsilon\right)$ as $u\!\to\!u_{c}^{j}$
originated in the valence band, which moves towards $\varepsilon\!=\!0$
and becomes sharper while always keeping $\delta_{j}\!\left(0,u\right)\!=\!0$.\,At
$u\!=\!u_{c}$, the situation is delicate because $\delta_{j}\!\left(\varepsilon\right)$
is no longer differentiable at $\varepsilon\!=\!0$.\,In that case,
one must work with Eq.\,\eqref{eq:Number of States Variation-2}
directly, and since there is a zero-energy $\pi$ discontinuity in
the phase shifts, a Dirac-$\delta$ distribution around the $u_{\text{c},n}^{j}$
emerges as the limit of plateau functions with a conserved integral
equal to $4$\,(the degeneracy of the $j\!=\!1/2$ single-impurity
bound states).\,This limit is depicted in Fig.\,\ref{fig:LowEnergyDoS_Altland}(a).\,An
immediate implication is that a DWSM is unstable to dilute random
impurities provided $p(u_{c}^{j})\!\neq\!0$ for at least one critical
$u_{c}^{j}$.\,Such a condition implies that the stability of a 3D
semimetallic phase ultimately depends on the type of impurity model
and the disorder statistics, i.e., whether it supports the resonant
mechanism driven by a continuous distribution of near-critical impurities.
These findings, supported below by accurate lattice simulations, show
that dirty 3D DWSMs with near-critical impurities are inherently unstable,
which sheds light on the previously reported AQC\,\citep{Pixley2016_2,pixley2016_3,Syzranov_2018,Ziegler2018,Goncalves2020,Wilson2020}.}

\textcolor{black}{In Fig.\,\ref{fig:LowEnergyDoS_Altland}(d), we
plot the change in the $j\!\!=\!\!1/2$ DoS due to a dilute diversity
of}\textit{\textcolor{black}{{} ``near-critical impurities}}\textcolor{black}{''.
The diversity is characterized by a Gaussian distribution $p(u)\!=\!\exp\left[-\left(u\!-\!\pi\right)/2\sigma^{2}\right]/\sqrt{2\pi}\sigma$
around $u_{\text{c}}\!=\!\pi$. The DoS is clearly lifted around $\varepsilon\!=\!0$,
forming a sharp symmetrical bump. For this diversity model, the peak
is Gaussian shaped near its center, and the corresponding area is
conserved as $\sigma\!\to\!0$. In this limit, a $4\delta\!\left(\varepsilon\right)$
distribution forms, i.e., all impurities are critical, each having
a fourfold degenerate zero-energy bound state.}

\vspace{-0.5cm}

\section{\label{sec:Lattice-simulations}Lattice simulations}

\vspace{-0.2cm}

Our prediction for the lifting of the DoS due to near-critical impurities
has been based on a continuum model for a single-node Dirac semimetal.
However, real Dirac materials and numerical simulations live in the
realm of lattice models, featuring several nodes and warped band structures.
To validate our previous conclusions, we perform real-space simulations
on a simple cubic lattice ($\mathcal{L}_{C}$ of parameter $a$ and
linear size $L$) with a four-orbital Hamiltonian derived from the
continuum Hamiltonian $\mathcal{H}$\,\citep{Pixley2016,Pixley2016_2,pixley2016_3},
namely,

\begin{figure}[t]
\begin{centering}
\hspace{-0.3cm}\includegraphics[scale=0.24]{Fig3.pdf}
\par\end{centering}
\vspace{-0.3cm}

\caption{\label{fig:LowEnergyDoS}(a) DoS change due to an impurity of critical
strength $u\!=\!\pi\,v/a$ and radius $16\,a$ inside a supercell
of $256^{3}$ sites. Vertical widths are $95\%$ statistical error
bars, and dashed lines are the continuum theory predictions. (b) Plot
of $\protect\av{\rho\!\left(E\!=\!0\right)}_{u}$ with several impurities
of radius $16\,a$ inside the simulated supercell of $512^{3}$ sites
for different resolutions $\eta$. The gray line is the dilute regime
prediction. The inset shows converged $\protect\av{\rho\!\left(E\right)}_{u}$
for three concentrations against the predictions of Fig.\,\ref{fig:LowEnergyDoS_Altland}(d)
(black lines).\vspace{-0.7cm}
}
\end{figure}

\vspace{-0.5cm}

\begin{align}
\negthickspace\!H & =\!\negthickspace\sum_{{\scriptscriptstyle \mathbf{R}\in\mathcal{L_{\text{C}}}}}\negthickspace\left[\frac{\imath v}{2a}\Psi_{\mathbf{{\scriptscriptstyle R}}}^{\dagger}\!\cdot\!\alpha^{j}\!\cdot\!\Psi_{{\scriptscriptstyle \mathbf{R}+a\hat{e}_{j}}}\!+\!\frac{\mathcal{U}\left(\mathbf{R}\right)}{2}\Psi_{\mathbf{{\scriptscriptstyle R}}}^{\dagger}\!\cdot\!\Psi_{\mathbf{{\scriptscriptstyle R}}}\!+\!\text{H.c.}\right]\!.\label{eq:LatticeModel}
\end{align}

\vspace{-0.2cm}

\noindent The DoS is calculated by means of accurate Chebyshev polynomial
expansions of the resolvent operator~$\delta(E\!-\!H)$\,\footnote{\textcolor{black}{In all numerical calculations, the variable $E$
is to be understood as a dimensionless energy, measured in units of
$a/\hbar v$.}} in very large systems as implemented in the open-source quantum transport
code $\mathtt{QUANTUM}\;\mathtt{KITE}$\,\citep{Joao2020}. The energy
resolution reads as $\eta\!=\!\pi\Delta E/2M$, where $\Delta E$
is the bandwidth of the Hamiltonian matrix and $M$ is the truncation
order of the polynomial expansion\,\citep{Weise2006,Ferreira2015}.
The DoS is obtained from $\rho(E)\!=\!\lim_{\eta\rightarrow0}\lim_{D\rightarrow\infty}D^{-1}\textrm{Tr}[\langle\delta_{\eta}(E\!-\!H)\rangle]$,
where $\langle...\rangle$ denotes disorder averaging and $D\!=\!4L^{3}$
is the Hilbert space dimension. \textcolor{black}{In order to simulate
systems with a vanishing mean-level spacing, thereby performing calculations
bounded only by $\eta$, we randomly sample over twisted boundaries\,\citep{Goncalves2020,Wilson2020}.
}This approach allowed the DoS to be calculated with unprecedented
working spectral resolutions as low as $\eta\!\simeq\!4\!\times\!10^{-4}v/a$\textcolor{black}{,
whose full convergence requires $M\!\approx\!30\,000$ polynomials.
}More technical details are provided in Appendix\,\ref{sec:AppendixD}.

\textcolor{black}{Figure\,\ref{fig:LowEnergyDoS}(a) shows the average
DoS induced by critical impurities, $\Delta\rho\left(E\right)\!=\!\rho_{\text{imp}}\left(E,u\!=\!\pi\right)\!-\!2E^{2}/\pi^{2}$,
in the dilute regime. The numerical data are compared with our analytical
results {[}Eq.\,\eqref{eq:Number of States Variation-2} and discussion
thereafter{]}, including the eightfold valley degeneracy, and properly
convoluted with Gaussian functions of width $\eta$ to mimic the finite
numerical spectral resolution. The lifting of the DoS at the node
and the underlying near-critical impurity mechanism are borne out
by the spectral calculations, which show excellent quantitative agreement
with the continuum theory, provided the impurity radius is large enough
{[}see Fig.\,\ref{fig:LowEnergyDoS}(a) and additional numerical
evidence in Appendix\,\ref{sec:AppendixD}{]}. In Fig.\,\ref{fig:LowEnergyDoS}(b),
we present an analogous calculation for a system having several impurities
inside the supercell. The impurities are placed randomly without superpositions,
and their strengths drawn from a Gaussian distribution $\text{\ensuremath{\mathcal{N}}}\!\!\left(\mu\!=\!\pi\,v/a\!\!,\sigma\!=\!0.3v/b\right)$.
The continuum theory prediction for the low-energy bump in the DoS
is reproduced in the diluted limit {[}see inset of Fig.\,\ref{fig:LowEnergyDoS}(b){]},
and the law $\rho\!\left(E\!=\!0\right)\!\propto\!c$ remains accurate
up to $10^{-6}$ impurities per unit cell. The overshooting for higher
concentrations is due to multi-impurity effects, which become more
effective as impurities are pushed closer together. }

\vspace{-0.3cm}

\section{\label{sec:Outlook}Conclusions and Outlook}

\vspace{-0.2cm}

\textcolor{black}{Here, we have shown that AQC must occur in 3D Dirac
semimetals having dilute short-range scalar impurities, if their random
parameters have a nonzero probability density at so-called }\textit{\textcolor{black}{magical}}\textcolor{black}{{}
values, where nodal bound states appear. These results were based
on a continuum formulation of the problem treated at the single-impurity
level and quantitatively confirmed by high-resolution lattice simulations
in a gapless }\textcolor{black}{\emph{multivalley}}\textcolor{black}{{}
Dirac model hosting $\approx10^{-9}\!\!\!\!-\!10^{-6}$ random scalar
impurities impurities per unit cell. The perfect agreement between
theory and numerical simulations gives confidence that the newly unveiled
resonant mechanism stemming from diverse near-critical impurities
is a crucial piece in the DWSM quantum criticality puzzle. Moreover,
disparities with previous work\,\citep{Buchhold2018,Buchhold2018_2}
are explained by the presence of physical $\pi$ jumps in the scattering
phase shifts that prevent a direct use of FSR for diverse impurities
around the aforementioned }\textit{\textcolor{black}{magical}}\textcolor{black}{{}
parameters. Similar conclusions are expected to hold for Weyl semimetals,
as scalar impurities do not couple the different Weyl sectors in the
infinite volume limit of our Dirac model. Meanwhile, our lattice nodal
DoS calculations show a crossover from a dilute regime at very low
impurity concentration (with DoS scaling linearly with $c$) to an
intermediate impurity density regime ($c\gtrsim10^{-6}$ impurities
per unit cell), where the DoS diverges from the analytical prediction.
This behavior can be traced to quantum-coherent multiple-impurity
scattering events, which are neglected in our continuum theory.}

\textcolor{black}{A related, but nontrivial, question concerns the
validity of these conclusions when dealing with lattice models having
uncorrelated on-site disorder. In light of our theory, as well as
earlier work\,\citep{Pixley2016_2,Ziegler2018}, one reasonably expects
the semimetallic phase to be unstable for unbounded distributions.
However, such systems with highly concentrated and atomic-sized (on-site)
impurities are exactly in the regime where the lattice results deviate
from continuum predictions, hinting at yet further subtleties when
relating the fate of disordered DWSM phases with rare-event bound
states.}

\vspace{-0.3cm}

\section{\textup{Acknowledgments}}

\vspace{-0.2cm}

We acknowledge support by the Portuguese Foundation for Science and
Technology through Strategic Funding No. UIDB/04650/2020, Projects
No. POCI-01-0145-FEDER-028887 (J.P.S.P. and J.M.V.P.L) and No. CEECIND/02936/2017
(B.A.), and Grant No. PD/BD/142774/2018 (J.P.S.P.). A.F. acknowledges
financial support from the Royal Society through a Royal Society University
Research Fellowship. The numerical calculations were performed on
the Viking Cluster, a high performance computing\,(HPC) facility
provided by the University of York. We are grateful for support from
the University of York's HPC service and their team. The authors are
grateful to A. Altland, M. Buchhold, M. Gonçalves, S.\,M. João, A.
Antunes, and J.\,M.\,B. Lopes dos Santos for valuable discussions.
The authors finally thank the anonymous referees for useful comments
which helped to improve the clarity of this work.

\vspace{-0.25cm}