\section{Introduction}

The transformer architecture, originally designed for the natural language processing (NLP) tasks~\cite{transformer}, has shown promising performance on various vision problems recently~\cite{dosovitskiy2020vit, touvron2020deit, liu2021swin, setr, TransTracking,yu2021pointr,rao2021dynamicvit,cheng2021maskformer}. Different from convolutional neural networks (CNNs), vision transformer models use self-attention layers to capture long-term dependencies, which are able to learn more diverse interactions between spatial locations. The pure multi-layer perceptrons (MLP) models~\cite{tolstikhin2021mlp, touvron2021resmlp} further simplify the vision transformers by replacing the self-attention layers with MLPs that are applied across spatial locations. Since fewer inductive biases are introduced, these two kinds of models have the potential to learn more generic and flexible interactions among spatial locations from raw data.

One primary challenge of applying self-attention and pure MLP models to vision tasks is the considerable computational complexity that grows quadratically as the number of tokens increases. Therefore, typical vision transformer style models usually consider a relatively small resolution for the intermediate features (\eg $14\times14$ tokens are extracted from the input images in both ViT~\cite{dosovitskiy2020vit} and MLP-Mixer~\cite{tolstikhin2021mlp}). This design may limit the applications of downstream dense prediction tasks like detection and segmentation. A possible solution is to replace the global self-attention with several local self-attention like Swin transformer~\cite{liu2021swin}. Despite the effectiveness in practice, local self-attention brings quite a few hand-made choices (\eg, window size, padding strategy, \etc) and limits the receptive field of each layer.

In this paper, we present a new conceptually simple yet computationally efficient architecture called Global Filter Network (\emph{GFNet}), which follows the trend of removing inductive biases from vision models while enjoying the log-linear complexity in computation. The basic idea behind our architecture is to learn the interactions among spatial locations in the frequency domain. Different from the self-attention mechanism in vision transformers and the fully connected layers in MLP models, the interactions among tokens are modeled as a set of learnable \emph{global filters} that are applied to the spectrum of the input features. Since the global filters are able to cover all the frequencies, our model can capture both long-term and short-term interactions. The filters are directly learned from the raw data without introducing human priors. Our architecture is largely based on the vision transformers only with some minimal modifications. We replace the self-attention sub-layer in vision transformers with three key operations: a 2D discrete Fourier transform to convert the input spatial features to the frequency domain, an element-wise multiplication between frequency-domain features and the global filters, and a 2D inverse Fourier transform to map the features back to the spatial domain. Since the Fourier transform is used to mix the information of different tokens, the global filter is much more efficient compared to the self-attention and MLP thanks to the $\mathcal{O}(L\log L)$ complexity of the fast Fourier transform algorithm (FFT)~\cite{cooley1965algorithmfft}. Benefiting from this, the proposed global filter layer is less sensitive to the token length $L$ and thus is compatible with larger feature maps and CNN-style hierarchical architectures \emph{without modifications}. The overall architecture of {GFNet} is illustrated in Figure~\ref{fig:intro}. We also compare our global filter with prevalent operations in deep vision models in Table~\ref{tab:intro}.

Our experiments on ImageNet verify the effectiveness of {GFNet}. With a similar architecture, our model outperform the recent vision transformer and MLP models including DeiT~\cite{touvron2020deit}, ResMLP~\cite{touvron2021resmlp} and gMLP~\cite{liu2021pay}. When using the hierarchical architecture, {GFNet} can further enlarge the gap.  {GFNet} also works well on downstream transfer learning and semantic segmentation tasks.   Our results demonstrate that GFNet can be a very competitive alternative to transformer-style models and CNNs in efficiency, generalization ability and robustness.

\begin{figure}[t]
\centering
\includegraphics[width=0.925\linewidth]{figs/fig1.pdf}
\caption{\textbf{The overall architecture of the Global Filter Network}. Our architecture is based on Vision Transformer (ViT) models with some minimal modifications. We replace the self-attention sub-layer with the proposed \emph{global filter layer}, which consists of three key operations: a 2D discrete Fourier transform to convert the input spatial features to the frequency domain, an element-wise multiplication between frequency-domain features and the global filters, and a 2D inverse Fourier transform to map the features back to the spatial domain. The efficient fast Fourier transform (FFT) enables us to learn arbitrary interactions among spatial locations with log-linear complexity.  }
\label{fig:intro}
\end{figure}

\begin{table}[t]
\centering \small
\newcolumntype{g}{>{\columncolor{Gray}}X}
\begin{tabularx}{\textwidth}{g>{\Centering}p{6cm}>{\Centering}X}
\toprule
 & Complexity (FLOPs) & \# Parameters \\
\midrule
Depthwise Convolution &  $\mathcal{O}(k^2HWD)$ & $k^2D$  \\
Self-Attention & $\mathcal{O}(HWD^2 + H^2W^2D)$ & $4D^2$  \\
Spatial MLP & $\mathcal{O}(H^2W^2D)$ & $H^2W^2$ \\
\midrule
\textbf{Global Filter} & $\mathcal{O}\left(HWD \lceil \log_2(HW) \rceil + HWD\right)$  & $HWD$ \\
\bottomrule
\end{tabularx} \vspace{5pt}
\caption{Comparisons of the proposed \emph{Global Filter} with prevalent operations in deep vision models. $H$, $W$ and $D$ are the height, width and the number of channels of the feature maps. $k$ is the kernel size of the convolution operation. The proposed global filter is much more efficient than self-attention and spatial MLP. }
\vspace{-10pt}
\label{tab:intro}
\end{table}

\section{Related works}
\paragrapha{Vision transformers. } Since Dosovitskiy~\etal~\cite{dosovitskiy2020vit} introduce transformers to the image classification and achieve a competitive performance compared to CNNs, transformers begin to exhibit their potential in various vision tasks~\cite{carion2020detr,TransTracking,setr}. Recently, there are a large number of works which aim to improve the transformers~\cite{touvron2020deit,touvron2021cait,liu2021swin,wu2021cvt,jiang2021token,gong2021improve,yuan2021t2t}. These works either seek for better training strategies~\cite{touvron2020deit,gong2021improve} or design better architectures~\cite{liu2021swin,wu2021cvt,yuan2021t2t} or both~\cite{touvron2021cait,gong2021improve}. However, most of the architecture modification of the transformers~\cite{wu2021cvt,jiang2021token,liu2021swin,yuan2021t2t} introduces additional inductive biases similar to CNNs. In this work, we only focus on the standard transformer architecture~\cite{dosovitskiy2020vit,touvron2020deit} and our goal is to replace the heavy self-attention layer ($\mathcal{O}(L^2)$) to an more efficient operation which can still model the interactions among different spatial locations without introducing the inductive biases associated with CNNs.

\paragrapha{MLP-like models. } More recently, there are several works that question the importance of self-attention in the vision transformers and propose to use MLP to replace the self-attention layer in the transformers~\cite{tolstikhin2021mlp,touvron2021resmlp,liu2021pay}. The MLP-Mixer~\cite{tolstikhin2021mlp} employs  MLPs to perform token mixing and channel mixing alternatively in each block. ResMLP~\cite{touvron2021resmlp} adopts a similar idea but substitutes the Layer Normalization with an Affine transformation for acceleration. The recently proposed gMLP~\cite{liu2021pay} uses a spatial gating unit to re-weight tokens in the spatial dimension. However, all of the above models include MLPs to mix the tokens spatially, which brings two drawbacks: (1) like the self-attention in the transformers, the spatial MLP still requires computational complexity quadratic to the length of tokens. (2) unlike transformers, MLP models are hard to scale up to higher resolution since the weights of the spatial MLPs have fixed sizes. Our work follows this trend and successfully resolves the above issues in MLP-like models. The proposed GFNet{} enjoys log-linear complexity and can be easily scaled up to any resolution.

\paragrapha{Applications of Fourier transform in vision. } Fourier transform has been an important tool in digital image processing for decades~\cite{pitas2000digital,baxes1994digital}.   With the breakthroughs of CNNs in vision~\cite{he2016deep,he2017mask}, there are a variety of works that start to incorporate Fourier transform in some deep learning method~\cite{li2020falcon,yang2020fda,ding2017circnn,lee2018single,chi2020fast} for vision tasks. Some of these works employ discrete Fourier transform to convert the images to the frequency domain and leverage the frequency information to improve the performance in certain tasks~\cite{lee2018single,yang2020fda}, while others utilize the convolution theorem to accelerate the CNNs via fast Fourier transform (FFT)~\cite{li2020falcon,ding2017circnn}. FFC~\cite{chi2020fast} replaces the convolution in CNNs with an Local Fourier Unit and perform convolutions in the frequency domain.  Very recent works also try to leverage Fourier transform to develop deep learning models to solve partial differential equations~\cite{li2020fourier} and NLP tasks~\cite{lee2021fnet}. In this work, we propose to use learnable filters to interchange information globally among the tokens in the Fourier domain, inspired by the frequency filters in the digital image processing~\cite{pitas2000digital}. We also take advantage of some properties of FFT to reduce the computational costs and the number of parameters.

\section{Method}
\subsection{Preliminaries: discrete Fourier transform}\label{sec:dft}
We start by introducing the discrete Fourier transform (DFT), which plays an important role in the area of digital signal processing and is a crucial component in our GFNet. For clarity, We first consider the 1D DFT. Given a sequence of $N$ complex numbers $x[n], 0\le n\le N-1$, the 1D DFT converts the sequence into the frequency domain by:
\begin{equation}
    X[k]=\sum_{n=0}^{N-1} x[n] e^{-j(2\pi / N) k n} := \sum_{n=0}^{N-1} x[n] W_N^{kn}
    \label{equ:dft}
\end{equation}
where $j$ is the imaginary unit and $W_N=e^{-j(2\pi/N)}$. The formulation of DFT in Equation~\eqref{equ:dft}~can be derived from the Fourier transform for continuous signal by sampling in both the time domain and the frequency domain (see Appendix~\ref{appendix:DFT} for details). Since $X[k]$ repeats on intervals of length $N$, it is suffice to take the value of $X[k]$ at $N$ consecutive points $k=0,1,\ldots, N-1$. Specifically, $X[k]$ represents to the spectrum of the sequence $x[n]$ at the frequency $\omega_k = 2\pi k/N$.

It is also worth noting that DFT is a one-to-one transformation. Given the DFT $X[k]$, we can recover the original signal $x[n]$ by the inverse DFT (IDFT):
\begin{equation}
    x[n] = \frac{1}{N} \sum_{k=0}^{N - 1}X[k]e^{j(2\pi/N)kn}.\label{equ:idft}
\end{equation}

For real input $x[n]$, it can be proved that (see Appendix~\ref{appendix:DFT}) its DFT is conjugate symmetric, i.e., $X[N - k] =X^*[k]$. The reverse is true as well: if we perform IDFT to $X[k]$ which is conjugate symmetric, a real discrete signal can be recovered. This property implies that the half of the DFT $\{X[k]: 0\le k\le \lceil N / 2\rceil\}$ contains the full information about the frequency characteristics of $x[n]$.

DFT is widely used in modern signal processing algorithms for mainly two reasons: (1) the input and output of DFT are both discrete thus can be easily processed by computers; (2) there exist efficient algorithms for computing the DFT. The \emph{fast Fourier transform} (FFT) algorithms take advantage of the symmetry and periodicity properties of $W_N^{kn}$ and reduce the complexity to compute DFT from $\mathcal{O}(N^2)$ to $\mathcal{O}(N\log N)$. The inverse DFT~\eqref{equ:idft}, which has a similar form to the DFT, can also be computed efficiently using the inverse fast Fourier transform (IFFT).

The DFT described above can be extend to 2D signals. Given the 2D signal $X[m, n], 0\le m\le M-1, 0\le n\le N-1$, the 2D DFT of $x[m, n]$ is given by:
\begin{equation}
    X[u, v] = \sum_{m=0}^{M-1}\sum_{n=0}^{N-1}x[m, n]e^{-j2\pi\left(\frac{um}{M}+\frac{vn}{N}\right)}.
\end{equation}
The 2D DFT can be viewed as performing 1D DFT on the two dimensions alternatively. Similar to 1D DFT, 2D DFT of real input $x[m, n]$ satisfied the conjugate symmetry property $X[M-u, N-v] = X^*[u, v]$. The FFT algorithms can also be applied to 2D DFT to improve computational efficiency.

\subsection{Global Filter Networks}\label{sec:gknet}
\paragrapha{Overall architecture.} Recent advances in vision transformers~\cite{dosovitskiy2020vit,touvron2020deit} demonstrate that models based on self-attention can achieve competitive performance even without the inductive biases associated with the convolutions. Henceforth, there are several works~\cite{touvron2021resmlp,tolstikhin2021mlp} that exploit approaches (\eg, MLPs) other than self-attention to mix the information among the tokens. The proposed Global Filter Networks (GFNet) follows this line of work and aims to replace the heavy self-attention layer ($\mathcal{O}(N^2)$) with a simpler and more efficient one.

The overall architecture of our model is depicted in Figure~\ref{fig:intro}. Our model takes as an input $H\times W$ non-overlapping patches and projects the flattened patches into $L=HW$ tokens with dimension $D$. The basic building block of GFNet{} consists of: 1) a \emph{global filter layer} that can exchange spatial information efficiently ($\mathcal{O}(L\log L)$); 2) a feedforward network (FFN) as in~\cite{dosovitskiy2020vit,touvron2020deit}. The output tokens of the last block are fed into a global average pooling layer followed by a linear classifier.

\begin{algorithm}[t]
\caption{Pseudocode of Global Filter Layer.}
\label{alg:code}
\algcomment{\fontsize{7.2pt}{0em}\selectfont \texttt{rfft2/irfft2}: 2D FFT/IFFT for real signal
}
\definecolor{codeblue}{rgb}{0.25,0.5,0.5}
\lstset{
  backgroundcolor=\color{white},
  basicstyle=\fontsize{7.2pt}{7.2pt}\ttfamily\selectfont,
  columns=fullflexible,
  breaklines=true,
  captionpos=b,
  commentstyle=\fontsize{7.2pt}{7.2pt}\color{codeblue},
  keywordstyle=\fontsize{7.2pt}{7.2pt},
}
\begin{lstlisting}[language=python]
# x: the token features, B x H x W x D (where N = H * W)
# K: the frequency-domain filter, H x W_hat x D (where W_hat = W // 2 + 1, see Section 3.2 for details)

X = rfft2(x, dim=(1, 2))
X_tilde = X * K
x = irfft2(X_tilde, dim=(1, 2))

\end{lstlisting}
\end{algorithm}

\paragrapha{Global filter layer.} We propose global filter layer as an alternative to the self-attention layer which can mix tokens representing different spatial locations. Given the tokens $\bm{x}\in \mathbb{R}^{H\times W\times D}$, we first perform 2D FFT (see Section~\ref{sec:dft}) along the spatial dimensions to convert $\bm{x}$ to the frequency domain:
\begin{equation}
    \bm{X}=\mathcal{F}[\bm{x}]\in \mathbb{C}^{H\times W\times D},
\end{equation}
where $\mathcal{F}[\cdot]$ denotes the 2D FFT. Note that $\bm{X}$ is a complex tensor and represents the spectrum of $\bm{x}$. We can then modulate the spectrum by multiplying a learnable filter $\bm{K}\in \mathbb{C}^{H\times W\times D}$ to the $\bm{X}$:
\begin{equation}
    \tilde{\bm{X}} = \bm{K}\odot \bm{X},
\end{equation}
where $\odot$ is the element-wise multiplication (also known as the Hadamard product). The filter $\bm{K}$ is called the \emph{global filter} since it has the same dimension with $\bm{X}$, which can represent an arbitrary filter in the frequency domain. Finally, we adopt the inverse FFT to transform the modulated spectrum $\tilde{X}$ back to the spatial domain and update the tokens:
\begin{equation}
    \bm{x}\leftarrow \mathcal{F}^{-1}[\tilde{\bm{X}}].
\end{equation}
The formulation of the global filter layer is motivated by the frequency filters in the digital image processing~\cite{pitas2000digital}, where the global filter $\bm{K}$ can be regarded as a set of learnable frequency filters for different hidden dimensions. It can be proved (see Appendix~\ref{appendix:DFT}) that the global filter layer is equivalent to a depthwise \emph{global circular convolution} with the filter size $H\times W$. Therefore, the global filter layer is different from the standard convolutional layer which adopts a relatively small filter size to enforce the inductive biases of the locality. We also find although the proposed global filter can also be interpreted as a spatial domain operation, the filters learned in our networks exhibit more clear patterns in the frequency domain than the spatial domain, which indicates our models tend to capture relation in the frequency domain instead of spatial domain (see Figure~\ref{fig:vis}).  Note that the global filter implemented in the frequency domain is also much more efficient compared to the spatial domain, which enjoys a complexity of $\mathcal{O}(DL\log L)$ while the vanilla depthwise global circular convolution in the spatial domain has $\mathcal{O}(DL^2)$ complexity.  We will also show that the global filter layer is better than its local convolution counterparts in the experiments.

It is also worth noting that in the implementation, we make use of the property of DFT to reduce the redundant computation. Since $\bm{x}$ is a real tensor, its DFT $\bm{X}$ is conjugate symmetric, \ie $\bm{X}[H - u, W - v, :]=\bm{X}^*[H, W, :]$. Therefore, we can take only the half of the values in the $\bm{X}$ but preserve the full information at the same time:
\begin{equation}
    \bm{X}_r = \bm{X}[:, 0:\widehat{W}]:=\mathcal{F}_r[\bm{x}], \quad \widehat{W}=\lceil W / 2\rceil,
\end{equation}
where $\mathcal{F}_r$ denotes the 2D FFT for real input. In this way, we can implement the global filter as $\bm{K}_r\in \mathbb{C}^{H\times \widehat{W}\times D}$, which can reduce half the parameters. This can also ensure $\mathcal{F}^{-1}_r[\bm{K}_r\odot \bm{X}_r]$ is a real tensor, thus it can be added directly to the input $\bm{x}$. The global filter layer can be easily in modern deep learning frameworks (\eg, PyTorch~\cite{paszke2019pytorch}), as is shown in Algorithm~\ref{alg:code}. The FFT and IFFT are well supported by GPU and CPU thanks to the acceleration libraries like \texttt{cuFFT} and \texttt{mkl-fft}, which makes our models perform well on hardware.

\paragrapha{Relationship to other transformer-style models.} The GFNet{} follows the line of research about the exploration of approaches to mix the tokens. Compared to existing architectures like vision transformers and pure MLP models, we exhibit that GFNet{} has several favorable properties: 1) GFNet{} is more efficient. The complexity of both the vision transformers~\cite{dosovitskiy2020vit,touvron2020deit,touvron2021cait} and the MLP models~\cite{tolstikhin2021mlp,touvron2021resmlp} is $\mathcal{O}(L^2)$. Different from them, global filter layer only consists an FFT ($\mathcal{O}(L\log L)$), an element-wise multiplication ($\mathcal{O}(L)$) and an IFFT ($\mathcal{O}(L\log L)$), which means the total computational complexity is $\mathcal{O}(L\log L)$. 2) Although pure MLP models are simpler compared to transformers, it is hard to fine-tune them on higher resolution (\eg{}, from $224\times 224$ resolution to $384\times 384$ resolution) since they can only process a fixed number of tokens. As opposed to pure MLP models, we will show that our GFNet{} can be easily scaled up to higher resolution. Our model is more flexible since both the FFT and the IFFT have no learnable parameters and can process sequences with arbitrary length. We can simply interpolate the global filter $\bm{K}$ to $\bm{K}'\in\mathbb{C}^{H'\times W'\times D}$ for different inputs, where $H'\times W'$ is the target size. The interpolation is reasonable due to the property of DFT. Each element of the global filter $\bm{K}[u, v]$ corresponds to the spectrum of the filter at $\omega_u=2\pi u/H, \omega_v=2\pi v/W$ and thus, the global filter $\bm{K}$ can be viewed as a sampling of a continuous spectrum $\bm{K}(\omega_u, \omega_v)$, where $\omega_u, \omega_v\in [0, 2\pi]$. Hence, changing the resolution is equivalent to changing the sampling interval of $\bm{K}(\omega_u, \omega_v)$. Therefore, we only need to perform interpolation to shift from one resolution to another.

We also notice recently a concurrent work FNet~\cite{lee2021fnet} leverages Fourier transform to mix tokens. Our work is distinct from FNet in three aspects: (1) FNet performs FFT to the input and directly adds the real part of the spectrum to the input tokens, which blends the information from different domains (spatial/frequency) together. On the other hand, GFNet{} draws motivation from the frequency filters, which is more reasonable. (2) FNet only keeps the real part of the spectrum. Note that the spectrum of real input is conjugate symmetric, which means the real part is exactly symmetric and thus contains redundant information. Our GFNet{}, however, utilizes this property to simplify the computation. (3) FNet is designed for NLP tasks, while our GFNet{} focuses on vision tasks. In our experiments, we also implement the FNet and show that our model outperforms it.

\paragrapha{Architecture variants.} Due to the limitation from the quadratic complexity in the self-attention, vision transformers~\cite{dosovitskiy2020vit,touvron2020deit} are usually designed to process a relatively small feature map (\eg, $14\times 14$). However, our GFNet{}, which enjoys log-linear complexity, avoids that problem. Since in our GFNet{} the computational costs do not grow such significantly when the feature map size increases, we can adopt a hierarchical architecture inspired by the success of CNNs~\cite{krizhevsky2012alex,he2016deep}. Generally speaking, we can start from a large feature map (\eg, $56\times 56$) and gradually perform downsampling after a few blocks. In this paper, we mainly investigate two kinds of variants of GFNet{}: transformer-style models with a fixed number of tokens in each block and CNN-style hierarchical models with gradually downsampled tokens. For transformer-style models, we begin with a 12-layer model (\emph{GFNet-XS}) with a similar architecture with DeiT-S and ResMLP-12. Then, we obtain 3 variants of the model (\emph{GFNet-Ti}, \emph{GFNet-S} and \emph{GFNet-B}) by simply adjusting the depth and embedding dimension, which have similar computational costs with ResNet-18, 50 and 101~\cite{he2016deep}. For hierarchical models, we also design three models  (\emph{GFNet-H-Ti}, \emph{GFNet-H-S} and \emph{GFNet-H-B}) that have these three levels of complexity following the design of PVT~\cite{wang2021pyramid}. We use $4\times 4$ patch embedding to form the input tokens and use a non-overlapping convolution layer to downsample tokens  following~\cite{wang2021pyramid,liu2021swin}.  Unlike PVT~\cite{wang2021pyramid} and Swin~\cite{liu2021swin}, we directly apply our building block on different stages without any modifications.  The detailed architectures are summarized in Table~\ref{tab:arch}.

\begin{table}[t]
  \centering
\caption{Detailed configurations of different variants of GFNet. For hierarchical models, we provide the number of channels and blocks in 4 stages. The FLOPs are calculated with $224\times 224$ input.} \vspace{5pt}
    \begin{tabu}to\textwidth{l*{2}{X[c]}rr}\toprule
    Model & \#Blocks &  \#Channels & Params (M) & FLOPs (G)  \\ \midrule
    GFNet-Ti &  12 & 256  & 7 & 1.3\\
    GFNet-XS &  12 & 384  & 16 & 2.9\\
    GFNet-S &  19 & 384  & 25 & 4.5\\
    GFNet-B &  19 & 512  & 43 & 7.9\\
      \midrule
    GFNet-H-Ti &  [3, 3, 10, 3] & [64, 128, 256, 512]  & 15 & 2.1\\
    GFNet-H-S &  [3, 3, 10, 3] & [96, 192, 384, 768]  & 32 & 4.6\\
    GFNet-H-B &  [3, 3, 27, 3] & [96, 192, 384, 768]  & 54 & 8.6\\
      \bottomrule
    \end{tabu}%
  \label{tab:arch}
\end{table}%

\section{Experiments}

We conduct extensive experiments to verify the effectiveness of our GFNet. We present the main results on ImageNet~\cite{deng2009imagenet} and compare them with various architectures. We also test our models on the downstream transfer learning datasets including CIFAR-10/100~\cite{cifar}, Stanford Cars~\cite{cars} and  Flowers-102~\cite{flower}. Lastly, we investigate the efficiency and robustness of the proposed models and provide visualization to have an intuitive understanding of our method.

\begin{table}[t]
  \centering
\caption{\textbf{Comparisons with transformer-style architectures on ImageNet.} We compare different transformer-style architectures for image classification including vision transformers~\cite{touvron2020deit}, MLP-like models~\cite{touvron2021resmlp,liu2021pay} and our models that have comparable FLOPs and the number of parameters. We report the top-1 accuracy on the validation set of ImageNet as well as the number of parameters and FLOPs. All of our models are trained with $224\times 224$ images. We use ``↑384" to represent models finetuned on $384\times 384$ images for 30 epochs. } \vspace{5pt}
    \begin{tabular}{lrrrrrr}\toprule
    Model & Params (M) & FLOPs (G) & Resolution & Top-1 Acc. (\%) & Top-5  Acc.  (\%) \\ \midrule
     DeiT-Ti~\cite{touvron2020deit} & 5    & 1.2  & 224   & 72.2 &  91.1  \\
     gMLP-Ti~\cite{liu2021pay} & 6     & 1.4   & 224   & 72.0  & - \\
     \rowcolor{Gray} GFNet-Ti  & 7 & 1.3 & 224  & 74.6 & 92.2 \\ \midrule
     ResMLP-12~\cite{touvron2021resmlp} & 15    & 3.0   & 224   & 76.6  & - \\
     \rowcolor{Gray} GFNet-XS & 16 & 2.9 & 224 & 78.6 & 94.2\\ \midrule
     DeiT-S~\cite{touvron2020deit} & 22    & 4.6   & 224   & 79.8  & 95.0  \\
     gMLP-S~\cite{liu2021pay} & 20    & 4.5   & 224   & 79.4  & - \\
     \rowcolor{Gray} GFNet-S & 25 & 4.5 & 224 & 80.0 & 94.9\\ \midrule
     ResMLP-36~\cite{touvron2021resmlp} & 45    & 8.9   & 224   & 79.7  & - \\
     \rowcolor{Gray} GFNet-B & 43 & 7.9 & 224 & 80.7 & 95.1 \\
     \rowcolor{Gray} GFNet-XS↑384 & 18 & 8.4& 384 &  80.6 & 95.4\\ \midrule
     DeiT-B~\cite{touvron2020deit} & 86    & 17.5  & 224   & 81.8  & 95.6 \\
     gMLP-B~\cite{liu2021pay} & 73    & 15.8  & 224   & 81.6  & - \\
     \rowcolor{Gray} GFNet-S↑384 & 28 &  13.2 & 384 & 81.7 & 95.8\\
     \rowcolor{Gray} GFNet-B↑384  & 47 & 23.3 & 384 & 82.1 & 95.8\\ \bottomrule
    \end{tabular}%
  \label{tab:main1}
\end{table}%

\begin{table}[t]
  \centering
\caption{\textbf{Comparisons with hierarchical architectures on ImageNet.} We compare different hierarchical architectures for image classification including convolutional neural networks~\cite{he2016deep,regnet}, hierarchical vision transformers~\cite{wang2021pyramid,liu2021swin} and our hierarchical models that have comparable FLOPs and number of parameters. We report the top-1 accuracy on the validation set of ImageNet as well as the number of parameters and FLOPs. All models are trained and tested with $224\times 224$ images.} \vspace{5pt}
    \begin{tabular}{lrrrrr}\toprule
    Model & Params (M) & FLOPs (G)  & Top-1 Acc. (\%) & Top-5  Acc.  (\%) \\ \midrule
     ResNet-18~\cite{he2016deep} & 12    & 1.8    &  69.8 & 89.1  \\
     RegNetY-1.6GF~\cite{regnet} & 11    & 1.6      & 78.0  & - \\
     PVT-Ti~\cite{liu2021pay} & 13     & 1.9     & 75.1  & - \\
     \rowcolor{Gray} GFNet-H-Ti & 15 & 2.1 & 80.1 & 95.1 \\ \midrule
     ResNet-50~\cite{touvron2020deit} & 26    & 4.1    &  76.1 & 92.9  \\
     RegNetY-4.0GF~\cite{regnet} &  21    & 4.0     & 80.0  & -\\
     PVT-S~\cite{liu2021pay} & 25    & 3.8    & 79.8 & - \\
     Swin-Ti~\cite{liu2021swin} & 29 & 4.5 & 81.3 & - \\
     \rowcolor{Gray} GFNet-H-S & 32 & 4.6  &  81.5 & 95.6\\ \midrule
     ResNet-101~\cite{touvron2020deit} & 45    & 7.9      &   77.4& 93.5   \\
     RegNetY-8.0GF~\cite{regnet} & 39    & 8.0     & 81.7  & - \\
     PVT-M~\cite{liu2021pay} & 44    & 6.7     & 81.2  & - \\
     Swin-S~\cite{liu2021swin} & 50 & 8.7 & 83.0 & -\\
     \rowcolor{Gray} GFNet-H-B & 54 &  8.6 & 82.9 & 96.2 \\ \bottomrule
    \end{tabular}%
  \label{tab:main2} \vspace{-7pt}
\end{table}%

\subsection{ImageNet Classification}

\paragrapha{Setups.} We conduct our main experiments on ImageNet~\cite{deng2009imagenet}, which is a widely used large-scale benchmark for image classification. ImageNet contains roughly 1.2M images from 1,000 categories. Following common practice~\cite{he2016deep,touvron2020deit}, we train our models on the training set of ImageNet and report the single-crop top-1 accuracy on 50,000 validation images. To fairly compare with previous works~\cite{touvron2020deit,touvron2021resmlp}, we follow the most training details for our models and do not add extra regularization methods like~\cite{jiang2021token}. Different from~\cite{touvron2020deit}, we does not use EMA model~\cite{ema}, RandomEarse~\cite{zhong2020randomerasing} and repeated augmentation~\cite{repeataug}, which are important to train DeiT while sightly hurting the performance of our models. We set the gradient clipping norm to 1 for all of our models. During finetuning at the higher resolution, we use the hyper-parameters suggested by the implementation of~\cite{touvron2020deit} and train the model for 30 epochs. All of our models are trained on a single machine with 8 GPUs. More details can be found in Appendix~\ref{sec:train_details}.

\paragrapha{Comparisons with transformer-style architectures.} The results are presented in Table~\ref{tab:main1}. We compare our method with different  transformer-style architectures for image classification including vision transformers (DeiT~\cite{touvron2020deit}) and MLP-like models (ResMLP~\cite{touvron2021resmlp} and gMLP~\cite{liu2021pay}) that have similar complexity and number of parameters. We see that our method can clearly outperform recent MLP-like models such as ResMLP~\cite{touvron2021resmlp} and gMLP~\cite{liu2021pay}, and show similar performance to DeiT. Specifically, GFNet-XS outperforms ResMLP-12 by 2.0\% while having slightly fewer FLOPs. GFNet-S also achieves better top-1 accuracy compared to gMLP-S and DeiT-S. Our tiny model is significantly better compared to both DeiT-Ti (+2.4\%) and gMLP-Ti (+2.6\%) with the similar level of complexity.

\paragrapha{Comparisons with hierarchical architectures.} We compare different kinds of hierarchical models in Figure~\ref{tab:main2}. ResNet~\cite{he2016deep} is the most widely used convolutional model while RegNet~\cite{regnet} is a family of carefully designed CNN models. We also compare with recent hierarchical vision transformers PVT~\cite{wang2021pyramid} and Swin~\cite{liu2021swin}. Benefiting from the log-linear complexity, GFNet-H models show significantly better performance than ResNet, RegNet and PVT and achieve similar performance with Swin while having a much simpler and more generic design.

\paragrapha{Fine-tuning at higher resolution.} One prominent problem of MLP-like models is that the feature resolution is not adjustable. On the contrary, the proposed global filter is more flexible. We demonstrate the advantage of GFNet by finetuning the model trained at $224\times 224$ resolution to higher resolution following the practice in vision transformers~\cite{touvron2020deit}. As shown in Table~\ref{tab:main1}, our model can easily adapt to higher resolution with only 30 epoch finetuning and achieve better performance.

\begin{table}[htbp]
  \centering
  \caption{\textbf{Results on transfer learning datasets}. We report the top-1 accuracy on the four datasets as well as the number of parameters and FLOPs. } \vspace{5pt}
  \adjustbox{width=\textwidth}{
    \begin{tabu}to 1.1\textwidth{lrr*{4}{X[r]}}\toprule
    Model  & FLOPs & Params & CIFAR-10 & CIFAR-100 & Flowers-102 & Cars-196 \\\midrule
    ResNet50~\cite{he2016deep} & 4.1G   & 26M   & - & - & 96.2  & 90.0 \\
    EfficientNet-B7~\cite{tan2019efficientnet} & 37G    & 66M   & 98.9  & 91.7  & 98.8  & 94.7 \\
    ViT-B/16~\cite{dosovitskiy2020vit} & 55.4G  & 86M   & 98.1  & 87.1  & 89.5  & - \\
    ViT-L/16~\cite{dosovitskiy2020vit} & 190.7G & 307M     & 97.9  & 86.4  & 89.7  & - \\
    Deit-B/16~\cite{touvron2020deit} & 17.5G  & 86M     & 99.1  & 90.8  & 98.4  & 92.1 \\
    ResMLP-12~\cite{touvron2021resmlp} & 3.0G     & 15M    & 98.1  & 87.0    & 97.4  & 84.6 \\
    ResMLP-24~\cite{touvron2021resmlp} & 6.0G     & 30M    & 98.7  & 89.5  & 97.9  & 89.5 \\\midrule
     GFNet-XS & 2.9G   & 16M  &  98.6  &  89.1  &   98.1  &  92.8 \\
     GFNet-H-B & 8.6G   &  54M  &  99.0  &  90.3  &   98.8  & 93.2 \\\bottomrule
    \end{tabu}%
    } \vspace{-5pt}
  \label{tab:transfer
\end{table}%

\subsection{Transfer learning}

To test the generality of our architecture and the learned representation, we evaluate GFNet on a set of commonly used transfer learning benchmark datasets including CIFAR-10~\cite{cifar}, CIFAR-100~\cite{cifar}, Stanford Cars~\cite{cars} and Flowers-102~\cite{flower}. We follow the setting of previous works~\cite{tan2019efficientnet,dosovitskiy2020vit,touvron2020deit,touvron2021resmlp}, where the model is initialized by the ImageNet pre-trained weights and finetuned on the new datasets.
We evaluate the transfer learning performance of our basic model and best model. The results are presented in Table~\ref{tab:transfer}. The proposed models generally work well on downstream datasets. GFNet models outperform ResMLP models by a large margin and achieve very competitive performance with state-of-the-art EfficientNet-B7. Our models also show competitive performance compared to state-of-the-art CNNs and vision transformers.

\subsection{Analysis and visualization}

\begin{figure}[t]
\begin{subfigure}{.33\textwidth}
\includegraphics[width=\linewidth]{figs/flops.pdf}
\caption{}
\end{subfigure}
\begin{subfigure}{.33\textwidth}
\includegraphics[width=\linewidth]{figs/latency.pdf}
\caption{}
\end{subfigure}
\begin{subfigure}{.33\textwidth}
\includegraphics[width=\linewidth]{figs/memory.pdf}
\caption{}
\end{subfigure}
\caption{Comparisons among GFNet, ViT~\cite{dosovitskiy2020vit} and ResMLP~\cite{touvron2021resmlp} in \textbf{(a)} FLOPs \textbf{(b)} latency and \textbf{(c)} GPU memory with respect to the number of tokens (feature resolution). The dotted lines indicate the estimated values when the GPU memory has run out. The latency and GPU memory is measured using a single NVIDIA RTX 3090 GPU with batch size 32 and feature dimension 384. }
\label{fig:eff}
\end{figure}

\paragrapha{Efficiency of GFNet.} We demonstrate the efficiency of our GFNet in Figure~\ref{fig:eff}, where the models are compared in theoretical FLOPs, actual latency and peak memory usage on GPU. We test a single building block of each model (including one token mixing layer and one FFN) with respect to the different numbers of tokens and set the feature dimension and batch size to 384 and 32 respectively. The self-attention model quickly runs out of memory when feature resolution exceeds $56^2$, which is also the feature resolution of our hierarchical model. The advantage of the proposed architecture becomes larger as the resolution increases, which strongly shows the potential of our model in vision tasks requiring high-resolution feature maps.

\begin{figure}[t]
\centering
\begin{minipage}{0.48\textwidth}
\centering
\includegraphics[width=0.95\textwidth]{figs/tradeoff.pdf} \vspace{-5pt}
\caption{\small ImageNet acc.  \emph{vs} model complexity.}\label{fig:tradeoff}
\end{minipage}\hfill
\begin{minipage}{0.48\textwidth}
\centering
\captionsetup{type=table}
\centering
\caption{\small Comparisons among the GFNet{} and other variants based on the transformer-like architecture on ImageNet. We show that GFNet{} outperforms the ResMLP~\cite{touvron2021resmlp}, FNet~\cite{lee2021fnet} and models with local depthwise convolutions. We also report the number of parameters and theoretical complexity in FLOPs. } \vspace{5pt}
\adjustbox{width=0.95\linewidth}{
    \begin{tabular}{lccc}\toprule \small
    \multirow{2}{*}{Model} & Acc & Param & FLOPs \\
    & (\%) & (M) & (G) \\
    \midrule
    DeiT-S~\cite{touvron2020deit} & 79.8  & 22    & 4.6  \\
    \midrule
    Local  Conv ($3\times 3$) & 77.7  & 15     & 2.8  \\
    Local   Conv ($5\times 5$) & 78.1  & 15     & 2.9  \\
    Local   Conv ($7\times 7$) & 78.2  & 15     & 2.9  \\
    \midrule
    ResMLP~\cite{touvron2021resmlp} & 76.6  & 15     & 3.0  \\
    FNet~\cite{lee2021fnet}  & 71.2  & 15     & 2.9  \\
    \midrule
    GFNet-XS & 78.6  &  16   &  2.9  \\\bottomrule
    \end{tabular}%
}
\label{tab:ablation}
\end{minipage} \vspace{-10pt}
\end{figure}

\paragrapha{Complexity/accuracy trade-offs.} We show the computational complexity and accuracy trade-offs of various transformer-style architectures in Figure~\ref{fig:tradeoff}. It is clear that GFNet achieves the best trade-off among all kinds of models.

\paragrapha{Ablation study on the global filter.} To more clearly show the effectiveness of the proposed global filters, we compare GFNet-XS with several baseline models that are equipped with different token mixing operations. The results are presented in Table~\ref{tab:ablation}. All models have a similar building block ( token mixing layer + FFN ) and the same feature dimension of $D=384$. We also implement the recent FNet~\cite{lee2021fnet} for comparison, where a 1D FFT on feature dimension and a 2D FFT on spatial dimensions are used to mix tokens. As shown in Table~\ref{tab:ablation}, our method outperforms all baseline methods except DeiT-S that has 64\% higher FLOPs.

\paragrapha{Robustness \& generalization ability.} Inspired by the~\cite{mao2021rethinking}, we further conduct experiments to evaluate the robustness and the generalization ability of the GFNet.
For robustness, we consider ImageNet-A, ImageNet-C, FGSM and PGD. ImageNet-A~\cite{hendrycks2021natural_imageneta} (IN-A) is a challenging dataset that contains natural adversarial examples. ImageNet-C~\cite{hendrycks2019benchmarking_imagenetc} (IN-C) is used to validate the robustness of the model under various types of corruption. We use the mean corruption error (mCE, lower is better) on ImageNet-C as the evaluation metric. FGSM~\cite{goodfellow2014explaining_fgsm} and PGD~\cite{madry2017towards_pgd} are two widely used algorithms that are targeted to evaluate the adversarial robustness of the model by single-step attack and multi-step attack, respectively. For generalization ability, we adopt two variants of ImageNet validation set: ImageNet-V2~\cite{recht2019imagenet_v2} (IN-V2) and ImageNet-Real~\cite{beyer2020imagenet_real} (IN-Real). ImageNet-V2 is a re-collected version of ImageNet validation set following the same data collection procedure of ImageNet, while ImageNet-Real contains the same images as ImageNet validation set but has reassessed labels. We compare GFNet-S with various baselines in Table~\ref{tab:robustness} including CNNs, Transformers and MLP-like architectures and find the GFNet enjoys both favorable robustness and generalization ability.

\begin{table}
  \centering
  \caption{\textbf{Evaluation of robustness and generalization ability}. We measure the robustness from different aspects, including the adversarial robustness by adopting adversarial attack algorithms including FGSM and PGD and the performance on corrupted/out-of-distribution datasets including ImageNet-A~\cite{hendrycks2021natural_imageneta} (top-1 accuracy) and ImageNet-C~\cite{hendrycks2019benchmarking_imagenetc} (mCE, lower is better). The generalization ability is evaluated on ImageNet-V2~\cite{recht2019imagenet_v2} and ImageNet-Real~\cite{beyer2020imagenet_real}.} \vspace{5pt}
  \adjustbox{max width=\textwidth}{

    \begin{tabular}{l*{10}{c}}\toprule
    \multicolumn{1}{l}{\multirow{2}[0]{*}{Model}} & \multicolumn{1}{l}{FLOPs} & \multicolumn{1}{l}{Params} & \multicolumn{2}{c}{ImageNet} & \multicolumn{2}{c}{Generalization} & \multicolumn{4}{c}{Robustness} \\\cmidrule(lr){4-5} \cmidrule(lr){6-7} \cmidrule(lr){8-11}
          & \multicolumn{1}{c}{(G)} & \multicolumn{1}{c}{(M)} & \multicolumn{1}{c}{Top-1$\uparrow$} & \multicolumn{1}{c}{Top-5$\uparrow$} & \multicolumn{1}{c}{IN-V2$\uparrow$} & \multicolumn{1}{c}{IN-Real$\uparrow$} & \multicolumn{1}{c}{FGSM$\uparrow$} & \multicolumn{1}{c}{PGD$\uparrow$} & \multicolumn{1}{c}{IN-C$\downarrow$} & \multicolumn{1}{c}{IN-A$\uparrow$} \\\midrule
    ResNet-50~\cite{he2016deep} & 4.1   & 26  & 76.1  & 92.9   & 67.4  & 85.8  & 12.2  & 0.9   & 76.7  & 0.0 \\
    ResNeXt50-32x4d~\cite{xie2017aggregated} & 4.3   & 25    & 79.8  & 94.6  & 68.2  & 85.2  & 34.7  & 13.5  & 64.7  & 10.7 \\\midrule
    DeiT-S~\cite{touvron2020deit}  & 4.6   & 22    & 79.8  & 95.0  & 68.4  & 85.6  & 40.7  & 16.7  & 54.6  & 18.9 \\
    ResMLP-12~\cite{touvron2021resmlp} & 3.0   & 15    & 76.6  & 93.2  & 64.4  & 83.3  & 23.9  & 8.5   & 66.0  & 7.1 \\\midrule
    GFNet-S & 4.5   & 25    & 80.1  & 94.9  & 68.5  & 85.8  & 42.6  & 21.0  & 53.8  & 14.3 \\\bottomrule
    \end{tabular}%
    }
  \label{tab:robustness}%
\end{table}%

\begin{figure}[t]
\begin{subfigure}{.81\textwidth}
    \centering
    \includegraphics[height=5.1cm]{figs/freq_filter.pdf}
    \caption{Frequency domain}
\end{subfigure}%
\begin{subfigure}{.19\textwidth}
    \centering
    \includegraphics[height=5.1cm]{figs/time_filter.pdf}
    \caption{Spatial domain}
\end{subfigure}
\caption{Visualization of the learned \emph{global filters} in GFNet-XS. We visualize the original frequency domain global filters in (a) and show the corresponding spatial domain filters for the first 6 columns in (b). There are more clear patterns in the frequency domain than the spatial domain. }
\label{fig:vis} \vspace{-10pt}
\end{figure}

\paragrapha{Visualization.} The core operation in GFNet is the element-wise multiplication between frequency-domain features and the global filter. Therefore, it is easy to visualize and interpret. We visualize the frequency domain filters as well as their corresponding spatial domain filters in Figure~\ref{fig:vis}. The learned global filters have more clear patterns in the frequency domain, where different layers have different characteristics. Interestingly, the filters in the last layer particularly focus on the low-frequency component. The corresponding filters in the spatial domain are less interpretable for humans.

\section{Conclusion}
\label{sec:con}
We have presented the Global Filter Network (\emph{GFNet}), which is a conceptually simple yet computationally efficient architecture for image classification. Our model replaces the self-attention sub-layer in vision transformer with 2D FFT/IFFT and a set of learnable \emph{global filters} in the frequency domain. Benefiting from the token mixing operation with log-linear complexity, our architecture is highly efficient. Our experimental results demonstrated that {GFNet} can be a very competitive alternative to vision transformers, MLP-like models and CNNs in accuracy/complexity trade-offs.

\section*{Acknowledgment}

This work was supported in part by the National Key Research and Development Program of China under Grant 2017YFA0700802, in part by the National Natural Science Foundation of China under Grant 62125603, Grant 61822603, Grant U1813218, Grant U1713214, in part by Beijing Academy of Artificial Intelligence (BAAI), and in part by a grant from the Institute for Guo Qiang, Tsinghua University.

\bibliographystyle{plain}


\section{Introduction}

Observational surveys over the past years have uncovered the existence of supermassive black holes with masses of $\mathrm{M_{BH} \sim 10^9 M_{\odot}}$ in the first Gyr of the Universe. Today, a few hundred quasars are known at redshifts $z \gtrsim 6$, which may only represent the tip of the iceberg \citep[][and references therein]{Inayoshi_et_2020}. In addition to the two quasars known at $z \sim 7.5$ \citep{Banados_et_2018, Yang_et_2020}, the latest discovery of the most distant quasar at $z = 7.642$, with a black hole (BH) mass of $\mathrm{M_{BH} = (1.6 \pm 0.4) \times 10^9 M_{\odot}}$, is reported in \citet{Wang_et_2021}. The estimated Eddington ratios are high ($L/L_E \sim 0.7-1$), suggesting that these massive luminous quasars represent extreme objects accreting close to the Eddington limit at early times.

In contrast, the majority of the accreting black hole population in the early Universe are likely to have lower luminosities, either due to smaller masses or lower accretion rates. Such low-luminosity quasars are now starting to be detected in the Subaru High-z Exploration of Low-Luminosity Quasars (SHELLQs) survey \citep[][]{Matsuoka_et_2019}. In addition, many accreting black holes may be obscured by gas and dust, hence still awaiting to be discovered.

From a theoretical perspective, the mere existence of supermassive black holes (with $\mathrm{M_{BH} \sim 10^9 M_{\odot}}$) at the highest redshifts ($z > 6$) poses severe challenges to black hole growth models \citep[][and references therein]{Woods_et_2019}.
In the case of light seeds ($\mathrm{M_0 \sim 10^2 - 10^3 M_{\odot}}$, such as remnants of Population III stars), a continuous growth at the Eddington rate and/or several super-Eddington accretion episodes are required. Even in the case of heavy seeds ($\mathrm{M_{BH} \sim 10^4 - 10^5 M_{\odot}}$, such as direct collapse black holes or DCBH), a rapid and efficient growth is still required in order to account for the billion-solar masses observed at $z \gtrsim 6$. In any case, some form of supercritical accretion may be required at very high redshifts \citep[e.g.][]{Volonteri_Rees_2005}.

The Bondi accretion rate gives an order of magnitude estimate of how rapidly infalling gas can be supplied from large scales to the central BH, in the case of spherical symmetry and without any feedback. But, as a result of high accretion rates, a huge amount of radiation is released into the surrounding environment. The resulting radiative feedback can suppress the gas inflow, and the actual accretion rate may be limited to the Eddington rate defined by the balance between radiation pressure and gravity.

Numerical simulations suggest that the mass accretion rate can also be considerably reduced due to photoionisation feedback
\citep{Milosavljevic_et_2009, Park_Ricotti_2011, Park_Ricotti_2012}. In this picture, photons heat and ionise the surrounding gas, creating an expanding hot bubble (HII region) opposing the inflowing gas. 2D hydrodynamic simulations show that the accretion follows an oscillatory pattern with a quasi-periodic trend, due to both photoionisation heating and radiation pressure, with an average accretion rate of $\sim 0.2 \%$ of the Bondi rate \citep{Milosavljevic_et_2009}.

\citet{Park_Ricotti_2011} perform a suite of 1D and 2D hydrodynamic simulations with radiative transfer in spherical symmetry, and obtain an oscillatory behaviour for the accretion rate that is attributed to UV and X-ray photoheating. In this case, the mean accretion rate is $\sim 3 \%$ of the Bondi rate, with a steep dependence on the ambient gas temperature. A qualitatively similar picture is found in subsequent simulations, including the effects of angular momentum and radiation pressure, with a somewhat lower average accretion rate of $\sim 1\%$ of the Bondi rate \citep{Park_Ricotti_2012}. Overall, these numerical simulations indicate that the time-averaged accretion rate is mostly limited to a small fraction of the Bondi rate due to radiative feedback.

On the other hand, photon trapping can be realised for sufficiently high accretion rates. In this case, the trapping of photons within the accretion flow reduces the impact of radiative feedback, thus potentially allowing supercritical accretion \citep[][]{Wyithe_Loeb_2012}. In fact, if the emergent luminosity is limited by photon trapping ($L \lesssim L_E$), the central radiation can not much affect the gas dynamics on larger scales. Several episodes of super-Eddington accretion may then account for the rapid growth of massive black holes at high redshifts \citep{Volonteri_et_2015}. In this context, \citet{Pacucci_et_2017} discuss the most favourable conditions leading to effective BH growth at early times.

Recent radiation hydrodynamic (RHD) simulations suggest that a transition to a `hyper-Eddington' regime can occur at very high accretion rates in high-density environments. The required condition is that the size of the ionised HII region is smaller than the Bondi radius ($\mathrm{R_\mathrm{HII} < R_B}$, where $\mathrm{R_B = G M_{BH}/c_s^2}$ with the gas sound speed $c_s$); and 1D RHD simulations indicate that the resulting accretion rate can reach the Bondi value \citep{Inayoshi_et_2016}. This result seems to hold even for super-Eddington luminosities \citep{Sakurai_et_2016}.

Although most of the above mentioned numerical simulations consider primordial gas accretion, recent observations indicate the existence of metal-enriched gas at very high redshifts. Intriguingly, the host galaxies of the first quasars are found to contain large amounts of dust and metals. ALMA observations of a sample of $z > 6$ quasars indicate typical dust masses in the range $\mathrm{M_d \sim (10^7 - 10^9) M_{\odot}}$ \citep{Venemans_et_2018}. The common presence of dust grains at such high redshifts suggests early enrichment and a rapid build-up of the dust reservoir. Therefore the early growth of massive black holes is likely to occur in dusty environments.

In light of the latest observational results, namely the simultaneous detection of supermassive black holes and dust reservoirs at $z > 6$, here we explore the possibility of supercritical accretion occurring in the dusty environment of the first quasars. We generalise the concept of photon trapping to the case of dusty gas, and analyse whether dust photon trapping can be actually realised in the early Universe.

\section{Photon trapping in dusty gas}
\label{Subsection_photon_trapping}

Photon trapping occurs when the infall velocity of the accreting gas ($v_r$) exceeds the outward diffusion speed of the radiation ($c/\tau$) \citep{Begelman_1979}. The optical depth is given by $\tau = \kappa_e \rho r$, where $\kappa_\textrm{e}$ is the electron scattering opacity, and $\rho$ is the gas density related to the accretion rate by $\dot{M} = 4 \pi r^2 \rho v_r$.
The condition for the trapping of photons ($v_r > c/\tau$) is satisfied for radii smaller than a critical radius
\begin{equation}
R_\mathrm{tr,e} = \frac{\kappa_\mathrm{e} \dot{M}}{4 \pi c} .
\label{Eq_R_tre}
\end{equation}
This is the photon trapping radius, within which radiation is advected inwards. The standard photon trapping radius (equation \ref{Eq_R_tre}) is determined by electron scattering ($\kappa_e = \sigma_T/m_p\cong 0.4 \, \rm{{cm^2} g^{-1}}$, where $\sigma_T$ is the Thomson cross section) and only scales with the physical accretion rate $\dot{M}$. We note that the minimal condition for photon trapping is that the critical radius in equation \ref{Eq_R_tre} lies outside the Schwarzschild radius: $\mathrm{R_{tr,e} > R_S}$ (where $\mathrm{R_S = \frac{2 G M_{BH}}{c^2}}$), which implies a minimal accretion rate $\mathrm{\dot{M} > \dot{M}_e = \frac{8 \pi G M_{BH}}{\kappa_e c}}$. This can be compared to the Eddington accretion rate $\mathrm{\dot{M}_E = L_E/(\epsilon c^2)}$, where $\epsilon$ is the radiative efficiency. The resulting ratio is $\mathrm{\frac{\dot{M}_e}{\dot{M}_E} \sim 0.2 \left( \frac{\epsilon}{0.1} \right)}$, suggesting that the minimal condition for photon trapping due to electron scattering can be easily satisfied, even for sub-Eddington accretion rates.

Electron scattering opacity dominates in the innermost region close to the central BH, but dust opacity becomes more important on larger scales. While the inner zone is dust-free, dust grains can survive at larger radii where the temperature drops below the dust sublimation temperature. The corresponding dust sublimation radius is estimated by considering the energy balance between the absorbed and reradiated fluxes \citep[e.g.][]{Murray_et_2005} and is given by $R_\textrm{sub} = \sqrt{\frac{L}{4 \pi \sigma_\textrm{sb} T_\textrm{sub}^4}}$, where $L = f_\mathrm{E} L_\mathrm{E}$ is the luminosity expressed in units of the Eddington luminosity, $\sigma_\mathrm{sb}$ is the Stefan-Boltzmann constant, and $T_\mathrm{sub}$ is the sublimation temperature of dust grains.

By analogy with equation (\ref{Eq_R_tre}), we define a photon trapping radius for dusty gas, where the opacity is dominated by dust absorption ($\kappa_d = \sigma_d/m_p$, where $\sigma_d$ is the dust absorption cross section)
\begin{equation}
R_\mathrm{tr,d} = \frac{\kappa_d \dot{M}}{4 \pi c}  .
\label{Eq_R_trd}
\end{equation}
Since the radiation-matter coupling is enhanced in the presence of dust grains (with a typical ratio of $\kappa_d/\kappa_e \gtrsim 10^3$ in the ultraviolet band), the dust photon trapping radius is correspondingly larger: $R_\mathrm{tr,d} = (\kappa_d/\kappa_e) R_\mathrm{tr,e}$.

Dust photon trapping (DPT) can only be achieved if the actual trapping radius (equation \ref{Eq_R_trd}) exceeds the sublimation radius of dust grains
\begin{equation}
R_\mathrm{tr,d} > R_\mathrm{sub} .
\end{equation}
The above condition implies a minimal accretion rate for dusty gas: $\dot{M} > \dot{M}_\mathrm{d}$.
The critical accretion rate for DPT is given by
\begin{align}
&\dot{M}_\mathrm{d}
= \frac{4 \pi c}{\kappa_\mathrm{d}} \sqrt{ \frac{G c f_\mathrm{E} M_\mathrm{BH}}{\kappa_\mathrm{e} \sigma_\mathrm{sb} T_\mathrm{sub}^4} }   \nonumber \\
&\cong 0.3 \frac{M_{\odot}}{\mathrm{yr}} \, \left(\frac{\kappa_\mathrm{d}}{10^3 \mathrm{cm^2/g}}\right)^{-1} \left( \frac{f_\mathrm{E}}{1} \right)^\frac{1}{2}
\left( \frac{T_\mathrm{sub}}{1500 \mathrm{K}} \right)^{-2} \left( \frac{M_\mathrm{BH}}{10^5 M_{\odot}} \right)^\frac{1}{2}
\label{Eq_Mdotd}
\end{align}
In the following, we analyse under which physical conditions this critical DPT accretion rate can be achieved, and whether dust photon trapping may be actually realised.

Here we implicitly consider a quasi-spherical accretion assuming that the accreting gas has negligible angular momentum. The importance of low-angular momentum gas in promoting supercritical accretion has been emphasised by \citet{Volonteri_et_2015}.
We also note that spherical symmetry is adopted in 1D RHD simulations of hyper-Eddington BH growth \citep[e.g.][]{Inayoshi_et_2016}.

In more realistic situations, the infalling gas has non-zero angular momentum, leading to the formation of a compact accretion disc.
By conservation of specific angular momentum, the disc size may be expressed as $\mathrm{R_\textrm{disc} \sim \lambda^2 R_{\sigma} \sim \lambda^2 (G M_\mathrm{BH}/\sigma^2)}$, where $\lambda \leq 1$ is the fraction of preserved angular momentum \citep{Volonteri_et_2015}. The whole accretion disc radiation can be efficiently trapped if the disc radius is smaller than the photon trapping radius ($R_\mathrm{disc} < R_\mathrm{tr}$). We note that this condition is more easily satisfied for dust photon trapping than for electron scattering (since $\mathrm{R_{tr,d} \gg R_{tr,e}}$). In fact, for a typical value of $\lambda \sim 0.05$, the condition for efficient disc radiation trapping ($\mathrm{R_\mathrm{disc} < R_{tr,d}}$) is naturally satisfied for DPT accretion rates. Evidently, a lower $\lambda$ would imply a smaller accretion disc, further facilitating the trapping condition. A particular solution of supercritical accretion flows is the slim accretion disc model \citep[][]{Abramowicz_et_1988, Czerny_2019}. In the context of disc models, the trapping radius may be adjusted by a factor $h = H/R$, corresponding to the disc aspect ratio, which may be of order unity for geometrically thick discs puffed up by trapped radiation \citep{Ohsuga_et_2005}.

\section{Dust photon trapping: parameter space analysis}
\label{Section_DPT_parameterspace}

\begin{figure}
\begin{center}
\includegraphics[angle=0,width=0.45\textwidth]{plot_01.jpg}
\caption{Characteristic accretion rates as a function of black hole mass (with fiducial parameters, see the main text): Eddington rate (blue dotted), DPT rate (red solid), and Bondi rate (green dashed).
}
\label{Fig_dotM_MBH}
\end{center}
\end{figure}

\begin{figure*}
\begin{multicols}{2}
    \includegraphics[width=0.9\linewidth]{plot_02.jpg}\par
    \includegraphics[width=0.9\linewidth]{plot_03.jpg}\par
    \end{multicols}
\caption{DPT accretion rate normalised by the Bondi rate versus ambient gas density. DPT is realised for $\mathrm{\dot{M}_d/\dot{M}_B < 1}$, i.e. when the ratio falls below unity (marked by the horizontal yellow dotted line). Left panel: variations in the black hole mass; right panel: variations in the ambient gas temperature.
}
\label{Fig_MdotB_ninf}
\end{figure*}

The critical accretion rate for dust photon trapping (equation \ref{Eq_Mdotd}) can be compared with the Eddington accretion rate, and the ratio between the two is given by
\begin{align}
&\frac{\dot{M}_\mathrm{d}}{\dot{M}_\mathrm{E}} = \epsilon c^2 \frac{\kappa_\mathrm{e}}{\kappa_\mathrm{d}} \sqrt{ \frac{c f_\mathrm{E}}{\kappa_\mathrm{e} \sigma_\mathrm{sb} T_\mathrm{sub}^4 G M_\mathrm{BH}} } \\ \nonumber
&\cong 160 \left( \frac{\epsilon}{0.1} \right) \left( \frac{f_\mathrm{E}}{1} \right)^\frac{1}{2} \left( \frac{\kappa_\mathrm{d}}{10^3 \mathrm{cm^2/g}} \right)^{-1}  \left( \frac{T_\mathrm{sub}}{1500 \mathrm{K}} \right)^{-2} \left( \frac{M_\mathrm{BH}}{10^5 M_{\odot}} \right)^{-\frac{1}{2}} \nonumber
\label{Eq_Md_ME}
\end{align}

We see that the condition for dust photon trapping requires super-Eddington accretion for typical parameters.
Moreover, this implies that DPT accretion naturally includes photon trapping by electron scattering, while the converse is not true.

In order to realise dust photon trapping, the critical accretion rate must be smaller than the Bondi accretion rate.
The Bondi rate is given by $\dot{M}_B = (4 \pi m_p G^2 M_\mathrm{BH}^2 n_{\infty})/c_s^3$ , where $c_s = \sqrt{ \frac{k_B T_{\infty}}{\mu m_p}}$ is the sound speed of the gas, with $T_{\infty}$ and $n_{\infty}$ being the gas temperature and gas number density, respectively, and $\mu$ is the mean molecular weight ($\mu = 1.22$).
The ratio between the DPT and Bondi accretion rates is given by
\begin{align}
\frac{\dot{M}_{d}}{\dot{M}_B}
&= \frac{1}{m_p \kappa_d n_{\infty}} \left( \frac{f_E}{\kappa_e \sigma_{SB} T_{sub}^4} \right)^{1/2} \left( \frac{c k_B T_{\infty}}{ \mu m_p G M_{BH}} \right)^{3/2}  \nonumber \\
&\cong 0.03 \left( \frac{f_E}{1} \right)^{1/2} \left( \mathrm{\frac{\kappa_d}{10^3 cm^2/g}} \right)^{-1} \left( \mathrm{\frac{T_{sub}}{1500 K}} \right)^{-2} \nonumber \\
\quad \quad \quad &\left( \mathrm{\frac{T_{\infty}}{10^4 K}} \right)^{3/2} \left( \mathrm{\frac{n_{\infty}}{10^5 cm^{-3}}} \right)^{-1} \left( \frac{M_{BH}}{10^5 M_{\odot}} \right)^{-3/2}
\end{align}
We see that the ratio is smaller than unity ($\dot{M}_\textrm{d}/\dot{M}_\textrm{B} < 1$) for typical parameters, suggesting that DPT can be potentially realised.
The reference values for the ambient gas density and temperature are set to $n_{\infty} = 10^5 \mathrm{cm}^{-3}$ and $T_{\infty} = 10^4 \, \mathrm{K}$, as adopted in numerical simulations \citep{Inayoshi_et_2016, Takeo_et_2018}.

As mentioned in the Introduction, the mean accretion rate can be much lower than the Bondi rate due to photoionisation feedback, typically of the order of $\sim 1\%$ \citep{Park_Ricotti_2011, Park_Ricotti_2012}. In contrast, in sufficiently dense environments characterised by efficient photon trapping, the accretion rate can approach the Bondi value \citep{Inayoshi_et_2016}. Hence Bondi accretion rates may be achieved under certain physical conditions, without being much hampered by radiative feedback. Here we consider the simple case of gas accretion onto a static BH and consider Bondi accretion as the limiting case (for BHs moving through the surrounding medium, see e.g. \citet{Toyouchi_et_2020, Sugimura_Ricotti_2020}). We note in passing that the DPT accretion rate is a few percent of the Bondi rate for fiducial parameters (comparable to the radiation feedback-reduced Bondi rate).

In Figure $\ref{Fig_dotM_MBH}$, we compare the three critical accretion rates (Eddington rate, Bondi rate, and DPT rate) as a function of the black hole mass. We note their different scalings with BH mass: $\dot{M}_\textrm{E} \propto M_\textrm{BH}$, $\dot{M}_\textrm{B} \propto M_\textrm{BH}^2$, and $\dot{M}_\textrm{d} \propto \sqrt{M_\textrm{BH}}$, respectively.
We observe that the critical accretion rate for dust photon trapping is always larger than the Eddington rate, thus implying super-Eddington accretion. As the black hole mass increases, the difference between the two becomes less prominent (due to the shallower mass-dependence of the DPT rate, $\dot{M}_\textrm{d} \propto \sqrt{M_\textrm{BH}}$). The Bondi accretion rate exceeds the DPT rate above a certain critical mass ($\mathrm{M_\textrm{BH} \gtrsim 10^4 M_{\odot}}$); whereas for lower black hole masses, the DPT accretion rate is larger than the Bondi rate. Therefore, a certain minimal black hole mass is required in order to achieve dust photon trapping.
The condition $\dot{M}_\textrm{d}/\dot{M}_\textrm{B} < 1$ implies a critical BH mass of
\begin{align}
M_\mathrm{BH} &> \frac{c k_B T_{\infty}}{G \mu m_p} \left( \frac{1}{m_p \kappa_d n_{\infty}} \right)^{2/3} \left( \frac{f_E}{\kappa_e \sigma_{SB} T_{sub}^4} \right)^{1/3}  \nonumber \\
 &> 10^4 \, M_{\odot} \, \left( \frac{f_E}{1} \right)^{1/3} \left( \mathrm{\frac{\kappa_d}{10^3 cm^2/g}} \right)^{-2/3} \left( \mathrm{\frac{T_{sub}}{1500 K}} \right)^{-4/3} \nonumber \\&
\quad \quad \quad \quad \quad \left( \mathrm{\frac{T_{\infty}}{10^4 K}} \right) \left( \mathrm{\frac{n_{\infty}}{10^5 cm^{-3}}} \right)^{-2/3} \, ,
\label{Eq_MBH_d}
\end{align}
which corresponds to the critical mass observed in Figure $\ref{Fig_dotM_MBH}$.
We also note that the difference between the DPT and Bondi accretion rates increases with increasing BH mass (due to the steeper mass-dependence of the Bondi rate $\mathrm{\dot{M}_B \propto M_{BH}^2}$), hence the DPT condition is more easily satisfied at larger masses.

Equivalently, the DPT condition can also be expressed in terms of a critical ambient gas density as
\begin{align}
n_{\infty} > 3 \times 10^3 \mathrm{cm}^{-3} &\left( \frac{f_E}{1} \right)^{1/2}  \left( \mathrm{\frac{\kappa_d}{10^3 cm^2/g}} \right)^{-1} \left( \mathrm{\frac{T_{sub}}{1500K}} \right)^{-2} \nonumber \\
&\left( \mathrm{\frac{T_{\infty}}{10^4K}} \right)^{3/2} \left( \frac{M_{BH}}{10^5 M_{\odot}} \right)^{-3/2} \, .
\label{Eq_ninf_d}
\end{align}

In Figure \ref{Fig_MdotB_ninf} we plot the DPT accretion rate normalised by the Bondi rate ($\dot{M}_d/\dot{M}_B$) as a function of the ambient gas density. Dust photon trapping occurs when this ratio is smaller than unity. The left-hand panel of Figure \ref{Fig_MdotB_ninf} shows variations in the central BH mass. For the fiducial case with $\mathrm{M_{BH} = 10^5 M_{\odot}}$, the ratio becomes smaller than unity at gas densities larger than $\mathrm{n_{\infty} \gtrsim 3 \times 10^3 cm^{-3}}$ (cf. equation \ref{Eq_ninf_d}), where DPT can be realised.
For a smaller BH mass ($\mathrm{M_{BH} = 10^4 M_{\odot}}$), the gas density must exceed $\mathrm{n_{\infty} \gtrsim 10^5 cm^{-3}}$ to achieve DPT; while for a larger BH mass ($\mathrm{M_{BH} = 10^6 M_{\odot}}$), DPT can already be reached at lower densities of $\mathrm{n_{\infty} \gtrsim 10^2 cm^{-3}}$. Thus a larger BH mass allows the system to enter the DPT regime for lower gas densities and vice-versa. In fact, the $\dot{M}_{d}/\dot{M}_B$ ratio decreases with increasing gas density, scaling as $\propto 1/n_{\infty}$. This implies that at high gas densities, the DPT condition $\dot{M}_\textrm{d}/\dot{M}_\textrm{B} < 1$ is more easily satisfied, hence dust photon trapping can be realised for lower-mass black holes.
The right-hand panel of Figure \ref{Fig_MdotB_ninf} shows variations in the ambient gas temperature. In the fiducial case with $\mathrm{T_{\infty} = 10^4 K}$, the $\dot{M}_{d}/\dot{M}_B$ ratio becomes smaller than unity and DPT occurs for $\mathrm{n_{\infty} \gtrsim 3 \times 10^3 cm^{-3}}$.
We see that for a lower gas temperature, DPT can be reached at lower gas densities, while higher densities are required for a higher ambient gas temperature.

\begin{figure*}
\begin{multicols}{2}
    \includegraphics[width=0.9\linewidth]{plot_04.jpg}\par
    \includegraphics[width=0.9\linewidth]{plot_05.jpg}\par
    \end{multicols}
\caption{
The $\mathrm{M_{BH} - n_{\infty}}$ plane with critical boundaries (left panel): $\mathrm{R_{tr,e} = R_S}$ (blue dotted), $\mathrm{\dot{M}_{d} = \dot{M}_B}$ (red solid), and $\mathrm{R_{tr,d} = R_B}$ (green dashed).
The same $\mathrm{M_{BH} - n_{\infty}}$ plane with variations in the physical parameters (right panel): fiducial parameters (red solid), $T_{\infty} = 10^3$ K (cyan dashed), $\mathrm{T_{sub} = 2000}$ K (violet dash-dot), $\mathrm{\kappa_d = 500 \, cm^2 g^{-1}}$ (orange dotted).
}
\label{Fig_MBH_ninf_plane}
\end{figure*}

\section{The $M_\textrm{BH} - \MakeLowercase{n_{\infty}}$ plane}
\label{Section_MBH_n_plane}

Focusing on the two key parameters of the system, namely the central BH mass ($\mathrm{M_{BH}}$) and the host gas density ($n_{\infty}$), a two-dimensional plane can be defined \citep[e.g.][]{Pacucci_et_2017}.
In our case, the DPT condition $\dot{M}_\textrm{d}/\dot{M}_\textrm{B} < 1$ can be recast into the form
\begin{align}
&\left( \frac{M_{BH}}{10^5 M_{\odot}} \right)^{3/2} \left( \mathrm{\frac{n_{\infty}}{10^5 cm^{-3}}} \right) \\
\nonumber
&> 0.03 \left( \frac{f_E}{1} \right)^{1/2} \left( \mathrm{\frac{\kappa_d}{10^3 cm^2/g}} \right)^{-1}
\left( \mathrm{\frac{T_{sub}}{1500 K}} \right)^{-2}  \left( \mathrm{\frac{T_{\infty}}{10^4 K}} \right)^{3/2}
\label{Eq_MBH_ninf}
\end{align}

In Figure \ref{Fig_MBH_ninf_plane}, we show the resulting $M_\textrm{BH} - n_{\infty}$ plane. The red solid line defines the DPT condition ($\dot{M}_{d}/\dot{M}_B < 1$): dust photon trapping occurs if the source is located in the region of parameter space above this critical line. The minimal condition for photon trapping due to electron scattering ($\mathrm{R_{tr,e} > R_S}$) is indicated by the blue dotted line. At the other end, the maximal extent of the DPT zone may be estimated e.g. by assuming the Bondi accretion rate in the definition of the trapping radius (equation \ref{Eq_R_trd}). The corresponding maximal DPT radius is smaller than the Bondi radius for fiducial parameters, and here we take $\mathrm{R_{tr,d} \sim R_B}$ as the limiting condition for DPT (green dashed line). On the other side, if the dust photon trapping radius happens to be larger than the Bondi radius, the whole inner region ($r \lesssim R_\mathrm{B}$) could in principle be under trapped conditions, potentially leading to maximally efficient DPT. If we directly impose the physical condition $\mathrm{R_{tr,d} > R_B}$ for a generic mass accretion rate we have $\dot{M} > \frac{4 \pi c}{\kappa_d} \frac{G M_{BH}}{c_s^2}$. Compared to the Bondi rate, the resulting ratio would be $\dot{M}/\dot{M}_B \gtrsim 10$ for fiducial parameters, suggesting that such extreme form of accretion is rather unlikely to be realised. We note that dust photon trapping may be harder to achieve for smaller black hole masses in the early growth phase, but could be facilitated by the higher ambient gas densities at early times (see also Section  \ref{Section_Discussion}).

Figure \ref{Fig_MBH_ninf_plane} clearly illustrates that a small BH mass requires a high gas density to achieve dust photon trapping, while DPT can occur at lower gas densities for larger BH masses. We further consider how the $M_\textrm{BH} - n_{\infty}$ plane is affected by variations in the other physical parameters (right-hand panel of Figure \ref{Fig_MBH_ninf_plane}). We see that the ambient gas temperature can have a significant impact: for a lower gas temperature, the trapping condition is relaxed and the DPT area becomes larger (the critical line shifts to the lower-left). Similarly, the trapping condition is slightly relaxed for a higher dust sublimation temperature, whereas a lower dust opacity leads to a slight reduction in the DPT area.

\citet{Pacucci_et_2017} identify a `high-efficiency region' in the $M_\textrm{BH} - n_{\infty}$ plane, where a favourable combination of seed mass and gas density allows efficient accretion and thus rapid BH growth. Requiring that accretion proceeds efficiently all the way from large scales to small scales, they obtain that a black hole with mass $\mathrm{M \gtrsim 10^4 M_{\odot}}$ formed in the high-efficiency region can continuously accrete and grow. Their high-efficiency region is located towards the upper-right of the $M_\textrm{BH} - n_{\infty}$ plane, with a more massive BH being located in the high-efficiency region for a wider range of gas densities. This may roughly correspond to the DPT region represented in Figure \ref{Fig_MBH_ninf_plane}.

Assuming that photon trapping due to electron scattering limits the emerging luminosity to $\mathrm{L \lesssim L_E}$, 1D RHD simulations indicate that steady hyper-Eddington accretion can occur if the condition $(\mathrm{n_{\infty}/10^5 cm^{-3}}) > (\mathrm{M_{BH}/10^4 M_{\odot}})^{-1} (\mathrm{T_{\infty}/10^4 K})^{3/2}$ is satisfied \citep{Inayoshi_et_2016}. We note that this relation is similar to the DPT condition for $\mathrm{M_{BH} = 10^4 M_{\odot}}$, although the BH mass dependence is steeper in our case. Therefore, in broad agreement with previous studies, we obtain that more massive BHs grow more efficiently, as massive objects can accrete from a wider range of gas densities within the DPT regime.

\section{Obscured accretion and black hole mass growth}

\subsection{Obscured accretion}
\label{Subsection_Obscured_accretion}

Since photons are trapped within the accretion flow, dust photon trapping necessarily implies obscured accretion. This can have a considerable impact on the observational appearance of the source, such that the early growth of massive black holes may be hidden from view \citep{Wyithe_Loeb_2012, Volonteri_et_2015}.
A rough estimate of the column density can be obtained by assuming a fraction of the gas free-fall rate at the BH capture radius, as in \citet{Volonteri_et_2015}
\begin{equation}
\mathrm{N_H} \sim \frac{(\dot{M}/m_p) t_\mathrm{ff}}{r^2} \, ,
\end{equation}
where $t_\textrm{ff} = \sqrt{\frac{r^3}{2 G M_\mathrm{BH}}}$ is the free-fall timescale, and $r = G M_\mathrm{BH}/\sigma^2$ is the capture radius with velocity dispersion $\sigma$.
Assuming the critical DPT accretion rate (equation \ref{Eq_Mdotd}), the associated column density is given by
\begin{equation}
\mathrm{N_H}
\sim \frac{4 \pi c \sigma}{m_p \kappa_d} \sqrt{\frac{c f_E}{2 G M_{BH} \kappa_e \sigma_{SB} T_{sub}^4}} \, .
\end{equation}

Inserting the numerical values, we obtain
\begin{align}
\mathrm{N_H = 3.5 \times 10^{24} cm^{-2}} &\left( \frac{f_E}{1} \right)^{1/2} \left( \mathrm{\frac{\kappa_d}{10^3 cm^2/g}} \right)^{-1}  \left( \mathrm{\frac{\sigma}{50 km/s}} \right) \\
\nonumber
&\left( \mathrm{\frac{T_{sub}}{1500 K}} \right)^{-2} \left( \frac{M_{BH}}{10^5 M_{\odot}} \right)^{-1/2}
\label{Eq_NH}
\end{align}
For typical parameters, DPT accretion may indeed be obscured, close to the Compton thick limit ($\mathrm{N_H \sim 10^{24} cm^{-2}}$).
But we note that the column density decreases with increasing BH mass, scaling as $\mathrm{N_H \propto M_{BH}^{-1/2}}$.
For a massive BH with $\mathrm{M_{BH} \sim 10^9 M_{\odot}}$, the associated column density would be $\mathrm{N_H \sim 7 \times 10^{22} cm^{-2}}$, and the quasar could be unobscured \citep[cf.][]{Volonteri_et_2015}. Such sources may correspond to the luminous quasars detected at $z \gtrsim 6$.

\subsection{Black hole mass growth}
\label{subsection_BHmass_growth}

As dust photon trapping requires supercritical accretion rates, we expect DPT to be associated with a rapid build-up of the BH mass.
We recall that the standard Eddington accretion rate can be written as
\begin{equation}
\dot{M} = \frac{1 - \epsilon}{\epsilon} \frac{M}{t_E} \, ,
\label{Eq_dotME}
\end{equation}
where $t_E = \frac{\sigma_T c}{4 \pi G m_p}$ is the Eddington timescale.
Integrating equation (\ref{Eq_dotME}), we obtain the corresponding BH mass growth
\begin{equation}
M(t) = M_0 \exp \left( \frac{1-\epsilon}{\epsilon} \frac{t}{t_E} \right) \, ,
\end{equation}
where $M_0$ is the initial seed mass.

In a similar way, integrating the DPT accretion rate (equation \ref{Eq_Mdotd}), we obtain
\begin{equation}
M(t)
= \left( \sqrt{M_0} + \frac{2 \pi c}{\kappa_d} \sqrt{ \frac{G c f_E}{\kappa_e \sigma_{SB} T_{sub}^4}} t \right)^2 \, .
\label{Eq_MBH_DPT}
\end{equation}

Let us consider a simple numerical example of a seed BH accreting for a Salpeter time ($t \cong 5 \times 10^7$ yr). Starting from the same initial mass of $\mathrm{M_0 \sim 10^5 M_{\odot}}$, the final mass is $\mathrm{M_{BH} \sim 3 \times 10^5 M_{\odot}}$ for the Eddington accretion rate, while it is $\mathrm{M_{BH} \sim 7 \times 10^8 M_{\odot}}$ for the DPT accretion rate. In contrast to the case of Eddington-limited accretion, a BH can grow considerably under DPT conditions. Therefore DPT accretion may enable the growing BHs to reach large masses within short timescales, and may help explain the existence of supermassive BHs in the early Universe.

The DPT accretion rate (equation \ref{Eq_Mdotd}) is the minimal accretion rate required to achieve dust photon trapping; so the associated mass growth may be considered as a lower limit to the actual BH growth. However, in the presence of dust grains, one should also consider the effective Eddington limit for dusty gas \citep[e.g.][]{Fabian_2012}, which is lower than the standard Eddington luminosity. A more detailed analysis will be required to investigate the interplay between radiation pressure on dust and dust photon trapping.

\section{Discussion}
\label{Section_Discussion}

Recent observations show that the host galaxies of the first quasars at $z \gtrsim 6$ are chemically evolved systems, containing unexpectedly large amounts of dust and metals \citep[e.g.][and references therein]{Valiante_et_2017}. ALMA observations of a sample of quasars at $z \gtrsim 6$ indicate dust masses ranging from $\mathrm{M_d \sim 10^7 M_{\odot}}$ for the faintest sources to $\mathrm{M_d \sim 10^9 M_{\odot}}$ for the brightest objects \citep{Venemans_et_2018}.
A large dust mass of $\mathrm{M_d \sim 7 \times 10^7 M_{\odot}}$ is even observed in the highest-redshift quasar at $z = 7.6$ \citep{Wang_et_2021}. The presence of dust grains at such high redshifts requires rapid dust formation, possibly due to early supernovae and grain growth in the interstellar medium \citep{Lesniewska_Michalowski_2019}.
The above observational results point towards fast dust enrichment occurring in parallel to the early growth of massive black holes within the first Gyr of the Universe. This prompts us to consider the potential role of dust photon trapping in the supercritical growth of massive BHs at early times.

\subsection{Physical implications}
\label{subsection_physical_implications}

Photon trapping affects the radiative luminosity output, and thus can influence the accretion flow onto the growing BH, as well as determine the observational appearance of the source. Due to photon trapping, the emergent diffusive luminosity is usually limited to values lower than the Eddington luminosity, independently of the accretion rate. This can be seen from the simple relation $\mathrm{L \lesssim G M_\mathrm{BH} \dot{M}/R_\mathrm{tr} = L_E}$, whereby the accretion rate terms cancel out (as $R_\mathrm{{tr}} \propto \dot{M}$). Because the emergent luminosity does not exceed the Eddington limit, radiation pressure can not prevent the gas inflow, and the resulting accretion rates can in principle reach arbitrarily high values \citep{Begelman_1979}. This provides a physical mechanism potentially allowing super-Eddington accretion. Indeed, as the photons are trapped within the accretion flow, radiation pressure feedback becomes inefficient, such that the accretion rates can reach values far exceeding the Eddington rate \citep{Wyithe_Loeb_2012}.

When the emergent luminosity is limited by photon trapping ($\mathrm{L \lesssim L_E}$), the size of the ionised HII region is also likely confined inside the Bondi radius ($\mathrm{R_{HII} < R_B}$). As a consequence, photoionisation feedback can not operate efficiently and is unable to halt the inflowing gas, potentially leading to hyper-Eddington accretion \citep{Inayoshi_et_2016}. Moreover, the effect of photoionisation feedback may be further reduced for realistic accretion disc spectra, compared to the case of the single power-law spectrum usually assumed in most works \citep{Takeo_et_2019}.

As mentioned in Section \ref{Section_DPT_parameterspace}, for accretion rates satisfying the DPT condition, photon trapping due to electron scattering is also naturally realised. In the case of dust photon trapping, radiation feedback at radii larger than the dust sublimation radius is effectively suppressed, while the entire accretion disc radiation may be efficiently trapped. Under such extreme trapping conditions, most of the infalling matter may be accreted \citep{Volonteri_et_2015}. Since the accretion rate leading to dust photon trapping is much larger than the Eddington rate, it should be associated with super-Eddington accretion and a rapid build-up of the BH mass. For instance, adopting the DPT accretion rate as the minimal value, the enhancement factor could be of a few orders of magnitude compared to the Eddington-limited value.

It is known that photon trapping also modifies the emergent spectral energy distribution. Higher energy photons originating in the inner accretion flow are more efficiently trapped than lower energy photons produced from the outer regions.  As a result, the peak frequency of the spectral energy distribution shifts towards lower energies, and the AGN spectrum becomes softer with increasing accretion rate \citep{Ohsuga_et_2003}. In the case of DPT,  the whole disc radiation can be efficiently trapped, and the overall luminosity output may be considerably reduced, leading to faint AGNs. This may not be the case for photon trapping due to electron scattering, where only the inner parts of the accretion disc remain trapped (unless the accretion rate is extremely high).

Since UV/optical photons are trapped within the dusty accretion flow, the central AGN will be obscured. In this context, X-ray photons would be even more readily trapped than UV/optical photons, and thus super-Eddington accretion will be likely obscured at both X-ray and UV/optical wavelengths \citep{Wyithe_Loeb_2012}. In our picture, DPT accretion is indeed associated with significant obscuration, with large column densities expected at early times. Cosmological hydrodynamic simulations also suggest heavy obscuration on different physical scales around the growing black holes in the early Universe \citep{Trebitsch_et_2019, Ni_et_2020}. Assuming a standard radiative efficiency of $\epsilon \sim 0.1$, it has been estimated that the obscured fraction at $z > 7$ could be of the order of $\gtrsim 80\%$ \citep{Davies_et_2019}. In fact, much of the early growth of massive black holes could be obscured, and a large number of such obscured quasars may be missed in current optical surveys \citep[][and references therein]{Hickox_Alexander_2018, Inayoshi_et_2020}.

On the other hand, obscuration by dusty gas may cause higher energy photons to be reprocessed into lower energy photons, leading to an increase in the infrared (IR) emission. In high column density sources (e.g. Compton-thick AGNs), the X-ray emission will be significantly suppressed, while the emission in the IR band could be considerably enhanced. Future instruments may allow the detection of such obscured accreting BHs at very high redshifts \citep[][and references therein]{Pacucci_et_2019}.

\subsection{Comparison to other works}

Different forms of supercritical accretion have been proposed and discussed in the literature (see Introduction). While most previous works consider primordial gas composition, we discuss below a few other studies that explicitly analyse the accretion of dusty gas.

\citet{Wyithe_Loeb_2012} consider photon trapping in dusty accretion flows, and discuss the cosmological conditions allowing super-Eddington accretion. We note that their minimal accretion rate required for photon trapping is essentially equivalent to our DPT condition (equation \ref{Eq_Mdotd}). But they define the sublimation radius as a fixed multiple of the gravitational radius, while we assume a physically motivated sublimation radius, explicitly depending on the underlying parameters (Section \ref{Subsection_photon_trapping}).
The DPT condition may actually evolve with cosmic time, due to the redshift-dependence of the ambient gas density. While the number density is constrained by the halo properties in \citet{Wyithe_Loeb_2012}, here we leave $n_{\infty}$ as a free parameter and explore a wide range of ambient gas densities \citep[see also][]{Pacucci_et_2017}. Qualitatively, higher gas densities are to be expected at higher redshifts. As a consequence, when the ambient gas density was higher in the past, dust photon trapping could occur even for low-mass BHs, enabling efficient growth at early times. This suggests that the physical conditions in the early Universe were most favourable for the occurrence of DPT, allowing rapid BH growth at critical times.

 \citet{Yajima_et_2017} perform 1D RHD simulations of dusty gas accretion in the low density regime ($n_{\infty} \sim 10-10^3 \, \mathrm{cm^{-3}}$) to analyse the interplay between dust and photoionisation feedback. Compared to the accretion of primordial gas characterised by strong periodic bursts, the accretion pattern of dusty gas is found to be much smoother with smaller fluctuations. Radiation pressure on dust (neglecting the reprocessed IR emission) seems to regulate the overall accretion flow, while the effect of dust opacity (in reducing the size of the HII region) seems to be of secondary importance. In our framework, we would not expect efficient dust photon trapping to occur in such low density environments, unless the central BH is already quite massive (Section \ref{Section_MBH_n_plane}).

Extending the investigation to a higher density regime ($n_{\infty} \sim 10^5 \mathrm{cm^{-3}}$), \citet{Toyouchi_et_2019} carry out 1D RHD simulations of dusty gas to analyse the possible transition to hyper-Eddington accretion. By explicitly considering the radiative force due to reprocessed IR photons, which can penetrate far into the neutral medium beyond the HII bubble, the resulting accretion rates are reduced compared to the Bondi value. Moreover, the accretion rates are found to be lower for higher metallicities, suggesting that super-Eddington accretion is prevented above a certain critical metallicity. However, an opposite result is obtained in the newest 3D RHD simulations including radiation anisotropy and mass outflows, whereby higher metallicities lead to higher accretion rates \citep{Toyouchi_et_2020}.

In this work, we just consider a characteristic dust opacity at UV wavelengths ($\kappa_\mathrm{d,UV} \sim 10^3 \, \mathrm{cm^2} \mathrm{g^{-1}}$).  UV photons are most efficiently absorbed by dust grains, and subsequently reprocessed as IR photons (by energy conservation). The dust opacities in the IR are much lower than in the UV band, with a typical value of $\kappa_\mathrm{d,IR} \sim 5 \, \mathrm{cm^2}\mathrm{g^{-1}}$ scaling with the dust-to-gas ratio. As a consequence, the associated trapping radius is smaller by a factor $\sim (\kappa_\mathrm{d,IR}/\kappa_\mathrm{d,UV})$, and the critical accretion rate for IR photon trapping is correspondingly larger. Therefore the condition for the trapping of reprocessed IR radiation is likely more stringent, and would require accretion rates in excess of the Bondi rate. A more detailed analysis of the impact of such IR photon trapping is left for future work. In reality, the actual dust opacity also depends on grain size, structure, and composition, as well as incident radiation spectrum. Since the dust opacity grows in proportion to the dust-to-gas ratio, which in turn may scale with the metallicity, the critical DPT accretion rate is expected to be lower for higher metallicities. This suggests that dust photon trapping could be more easily realised in higher metallicity environments.

In addition, supercritical accretion can be facilitated by anisotropic feedback, either due to anisotropic emission from the accretion disc or shadowing effects on larger scales \citep{Sugimura_et_2017, Takeo_et_2018}. The exact radiation pattern is ultimately determined by the central BH spin, with low/high spins leading to prolate/oblate emission patterns, respectively \citep{Ishibashi_et_2019}. In reality, several competing mechanisms operate simultaneously, leading to opposing effects. For instance, the inclusion of the gas angular momentum and associated centrifugal barrier can hinder supercritical accretion \citep{Sugimura_et_2018}, but an accretion disc geometry may also favour the escape of IR photons in the polar directions, reducing the overall radiative feedback \citep{Toyouchi_et_2019}. Contrary to what might be expected, it has been shown that the transition to the hyper-Eddington accretion regime could even be facilitated by powerful outflows \citep{Takeo_et_2020}.

We have previously discussed how AGN radiative feedback, driven by radiation pressure on dust, may affect the surrounding environment through the development of galactic outflows \citep{Ishibashi_Fabian_2015, Ishibashi_et_2018}. We have also considered the effects of such radiation feedback-driven outflows on the early growth of massive BHs in the case of Eddington-limited accretion \citep{Ishibashi_2019}. A further step will be to investigate the dynamical coupling between radiative feedback, dust photon trapping, and supercritical BH growth in the early Universe.

\section*{Acknowledgements }

WI acknowledges support from the University of Zurich.

\section*{Data availability}

No new data were generated or analysed in support of this research.

\bibliographystyle{mn2e}