\section{Introduction}

Graph neural network (GNN) \citep{scarselli2008graph,bruna2013spectral,defferrard2016convolutional,kipf2016semi,velivckovic2017graph,hamilton2017inductive,gilmer2017neural,xu2018powerful} has become the prominent approach to learn representations for graphs, such as social networks \citep{li2019encoding}, biomedical information networks \citep{yan2019groupinn}, communication networks \citep{suarez2021graph}, n-body systems \citep{kipf2018neural}, etc. Most GNNs rely on the message passing mechanism \citep{gilmer2017neural} to implement the interactions between neighbouring nodes. Despite its huge success, message passing GNNs still faces two fundamental but fatal drawbacks: it fails to generalize to heterophily where neighboring nodes share dissimilar features or labels, and a simple multilayer perceptron can outperform many GNNs \citep{zhu2020beyond}, this limit GNNs extending to many real-world networks with heterophily;
it is also observed the node representations become indistinguishable when stacking multiple layers, and suffers sharply performance drop, resulting in the so-called ``over-smoothing'' problem \citep{li2018deeper}, which prevent GNNs to utilize high-order neighborhood information effectively.

To address these two drawbacks, numerous approaches have been proposed. Most of them concentrate on the aggregation stage of message passing. Some design signed messages to distinguish neighbors belong to different classes \citep{yang2021diverse,bo2021beyond,luan2021heterophily,yan2021two}, allowing GNNs to capture high-frequency signals; \citet{min2020scattering} design specific filters to capture band-pass signals; some apply personalized aggregation with reinforcement learning \citep{lai2020policy} or neural architecture search \citep{wang2022graph}; others attempt to aggregate messages not only from the direct neighbors, but also from the embedding space \citep{pei2020geom} or higher-order neighbors \citep{zhu2020beyond}. These aggregator designs have achieved good performance, however, they primarily focus on the single-round message passing process and ignore the integration of messages from multiple hops. Another line of works focus on the effective utilization of multiple hops information, which is mainly accomplished by designing various skip-connections. \citet{klicpera2018predict,chen2020simple} propose initial connection to prevent ego or local information from being ``washed out'' by stacking multiple GNN layers; inspired by ResNet \citep{he2016deep}, some works \citep{li2019deepgcns,chen2020simple,cong2021provable} explores the application of residual connection on GNNs to improve the gradients; others combined the output of intermediate GNN layers with well-designed components, such as concat \citep{xu2018representation,zhu2020beyond}, learnable weights \citep{zhu2020beyond,abu2019mixhop,liu2020towards}, signed weights \citep{chien2020adaptive}, or RNN-like architectures \citep{xu2018representation,sun2019adagcn}. These works are simple yet effective, however, they can only model the information \textit{within few hops}, but unable to model the information exactly \textit{at some orders}, this lead to a mixing of features at different orders;
besides, many of these approaches \citep{li2019deepgcns,abu2019mixhop,chen2020simple,zhu2020beyond,chien2020adaptive,cong2021provable} are unable to make personalized decisions for each node. These deficiencies result in suboptimal performance. In addition to caring about the model side, other approaches focus on how to modify the graph structure. These methods are called ``graph rewiring'', including randomly removing edges \citep{rong2019dropedge} or nodes \citep{feng2020graph}, or computing a new graph with heuristic algorithms \citep{suresh2021breaking,zeng2021decoupling}. In general, these algorithms are not learnable and thus only applicable to certain graphs.

\begin{wrapfigure}{r}{0.41\textwidth}
\vspace{-0.4cm}
\centering
\includegraphics[width=1.0\linewidth]{nested_tree.pdf}
\caption{Aligning the hierarchy of a rooted-tree $\mathcal{T}_{v}^{(k)}$ underlying the graph with the node embedding of the root node $v$. Neighboring nodes within $k$ hops of edges to $v$ naturally form a depth $k$ subtree. Messages passed to $v$ from nodes within this subtree are restricted to the first $P_v^{(k)}$ neurons in the node embedding of $v$.}
\label{fig.nested_tree}
\vspace{-0.45cm}
\end{wrapfigure}

Unlike the previous works, we address both problems by designing the combine stage of message passing and emphasize the importance of it.
\textit{The key idea is to integrate an inductive bias from rooted-tree hierarchy, let GNN encode the neighborhood information exactly at some orders and avoid feature mixing within hops.}
The combine stage has been rarely focused before, most works simply implement it as a self-loop. This would result in an unreasonable mixing of node features \citep{zhu2020beyond}.
To avoid this ``mixing'', \citet{hamilton2017inductive,xu2018representation,zhu2020beyond} concat the node representation and the aggregated message, which has been identified as an effective design to deal with heterophily \citep{zhu2020beyond}. However, keeping the embedding dimension constant across layers, the local information will be squeezed at an exponential rate. The most related work to ours is Gated GNN \citep{li2015gated}, it applys a GRU in the combine stage and strengthen the expressivity, but fails to prevent feature mixing, limiting the performance.

In this paper, we present the message passing mechanism in an ordered form. That is, the neurons in the node embedding of a certain node is aligned with the hierarchies in the rooted-tree of the node. Here, by rooted-tree of a node we refer to the tree with the node itself as the root and its neighbors as its children. Recursively, for each child, its children are again the neighboring nodes of the child. (c.f. Figure \ref{fig.nested_tree}). We achieve the alignment by proposing a novel ordered gating mechanism, which controls the assignment of neurons to encode subtrees with different depth. Experimental results on an extensive set of datasets show that our model could alleviate the heterophily and over-smoothing problem at the same time. Our model provides the following practical benefits:

\begin{itemize}
\item We design the combine stage guided by the rooted-tree hierarchy, a very general topological inductive bias with least assumption about the neighbors' distribution, this allows a flexible integration of information at different orders, and leading superior performance in both heterophily and homophily.

\item The ordered gating mechanism prevent the mixing of node features within hops, enable us to model information at different orders. This open a door to extract similar neighborhood patterns for each node under heterophily; this also make it easy to preserve ego and local information, then effectively alleviating over-smoothing.

\item Our model aligns neighboring structures with blocks in node embeddings through an explicit gating mechanism, thus the gating mechanism could provide visualizations to reveal the connecting type of the data and offer explainability.
\end{itemize}

\section{Our Approach}
In this section, we will first present a general formalization of GNNs. Then we show how to align the hierarchy of the rooted-tree to the node embeddings, the intuitions behind it, as well as the two necessary components to imeplement it. Finally, our model is presented in Section \ref{puttogether}, with few cases to explain it's practical benefits.

\subsection{Graph Neural Networks}
Consider an unweighted and self-loop-free graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$, with $\mathcal{V}$ and $\mathcal{E}$ denoting its nodes and edges, respectively. Let $v \in \mathcal{V}$ being one of the nodes with its label $y_{v} \in \mathcal{Y}$, and assume that there are $N$ nodes in the graph. The adjacency matrix for graph $\mathcal{G}$ is denoted as $A \in \mathbb{R}^{N \times N}$.
Assume that the input features for each node form a matrix $X \in \mathbb{R}^{N \times F}$, where $F$ is the dimension of input features, with $X_v \in \mathbb{R}^N$ being the features of node $v$. Before feeding it into the GNN, an input projection function $f_{\theta}(\cdot)$ is optionally incorporated to nonlinearly project the features into a space mathing the node embedding dimensions, i.e.,
\begin{equation}
Z_v=f_{\theta}(X_v)
\end{equation}
where $D$ is the dimension of node embeddings, and $Z_v \in \mathbb{R}^{D}$ is feed to the first layer of GNN, i.e., $h_{v}^{(0)}=Z_{v}$. In our work, we simply parameterize $f_{\theta}(\cdot)$ with an MLP.

Generally, a GNN layer with message passing can be divided into two stages: the \textit{aggregate} stage and the \textit{combine} stage. For the node $v$ at the $k$-th layer, it first aggregates information from its neighbors, forming the message vector $m_{v}^{(k)}$. Then in the second stage $m_{v}^{(k)}$ is combined with $v$'s ego representation $h_{v}^{(k-1)}$ from the last layer. Formally,
\begin{equation}
m_{v}^{(k)}=\operatorname{AGGREGATE}^{(k)}\left(\left\{h_{u}^{(k-1)}: u \in \mathcal{N}(v)\right\}\right), h_{v}^{(k)}=\operatorname{COMBINE}^{(k)}\left(h_{v}^{(k-1)}, m_{v}^{(k)}\right)
\end{equation}
where $\mathcal{N}(v)$ is the set of node $v$'s direct neighbors. The $\operatorname{AGGREGATE}(\cdot)$ operator denotes a permutation-invariant function on all the neighboring nodes of $v$, such as addition, mean pooling or max pooling, optionally followed by a transformation function. The transformation function could be an MLP \citep{xu2018powerful} or transformation matrix \citep{kipf2016semi}. The $\operatorname{COMBINE}(\cdot)$ operator denotes a combine function that fuses the ego information and the aggrated message.

Our model only modifies the \textit{combine} stage, which means it could be seamlessly extended to many popular GNNs, such as GCN \citep{kipf2016semi}, GAT \citep{velivckovic2017graph}, GIN \citep{xu2018powerful}, etc., by changing the \textit{aggregate} stage correspondingly. In this paper, we simply implement the \textit{aggregate} stage as computing the mean value of neighboring vectors.

\subsection{Aligning Rooted-Tree with Node Embedding}
\label{align_rooted_tree}
From the viewpoint of a node, the neighboring nodes within multiple hops around it can be naturally presented as a rooted-tree \citep{xu2018powerful,xu2018representation,liu2020towards}. For node $v$, its $k$-th order rooted-tree $\mathcal{T}_{v}^{(k)}$ is composed by putting its first order neighbors as its children, and then recursively for each child, putting the neighboring nodes of the child as its children (except the child that is also the ancester of the node), unitl its $k$-th order neighbors (Figure \ref{fig.nested_tree}).

Since there is only one round of message passing in each GNN layer, stacking multiple layers of GNNs corresponds to doing multiple rounds of message passing. And since for each single round of message passing, only direct neighbors exchange their information, the node representation of $v$ in the $k$-th layer $h_v^{(k)}$ is only depending on the nodes within its $k$-th order rooted-tree $\mathcal{T}_{v}^{(k)}$. Our goal is to properly encode the information within the $k$-th order rooted-tree into the fixed size vector $h_v^{(k)}$.

Note that there is a nested structure between rooted-trees of $v$, and we are going to utilize this structure as anchors to align the rooted-trees with the node embedding of $v$. It is obvious that the $k-1$-th rooted-tree $\mathcal{T}_{v}^{(k-1)}$ is a subtree of the $k$-th rooted-tree  $\mathcal{T}_{v}^{(k)}$, with the same root $v$ and the same tree structure up to depth $k-1$. Thus we have (assuming the model has $K$ layers)
\begin{equation} \label{relativemagnitude}
    \mathcal{T}_{v}^{(0)} \subseteq \mathcal{T}_{v}^{(1)} \subseteq \cdots \subseteq \mathcal{T}_{v}^{(k)} \subseteq \cdots \subseteq \mathcal{T}_{v}^{(K)}.
\end{equation}
where $\mathcal{T}_v^{(0)}$ corresponds to the neurons used to represent the node's ego information.
As $k$ grows, $\mathcal{T}_{v}^{(k)}$ becomes larger and more complex, requiring more neurons to encode its information. As a result, it is natural to assume that the neurons used to represent message passing within $\mathcal{T}_{v}^{(k-1)}$ should be a subset of those for $\mathcal{T}_{v}^{(k)}$. Correspondingly, we are going to implement the same nested structure in the node embedding of $v$. We first order all the neurons in $h_v$ into a sequence, thus we can index them with numbers. For information within $\mathcal{T}_{v}^{(k-1)}$, we allow it to be encoded within the first $P_v^{(k-1)}$ neurons. As for the next level rooted-tree $\mathcal{T}_{v}^{(k)}$, more neurons are involved, corresponding to the first $P_v^{(k)}$ neurons. Apparently here we have $P_v^{(k-1)} \leq P_v^{(k)}$.

There is a critical property that comes with this way of ordering the neurons. Consider the neurons between the two split points $P_v^{(k-1)}$ and $P_v^{(k)}$. They reflect the delta between $k-1$ and $k$ rounds of message passing, encoding the neighbor information exactly \textit{at $k$-th order}. Intuitively, it organizes the message passing in an \textit{ordered} form, avoid the mixing of neighbor features at different orders within a hop. We implement the ordered neurons with the gating mechanism (which will be introduced in the following sections in detail), thus it's natural to visualize the gates as a way to probe the message passing of the model.

Now we have the node embedding of $v$ split by a set of spliting points
\begin{equation} \label{relativemagnitude2}
    P_v^{(0)} \leq P_v^{(1)} \leq \cdots \leq P_v^{(k)} \leq \cdots \leq P_v^{(K)}.
\end{equation}
Each of them is predicted by a GNN layer, with $P_v^{(K)}=D$ pointing to the end of the node embedding. If we can learn the positions of these spliting points and ensure their relative magnitudes, then the alignment could be realized. We are going to show how we predict the spliting points through a gating mechanism, and how we lock their relative magnitudes with the differentiable OR operator. Our model naturally emerges by putting everything together in Section \ref{puttogether}.

\subsection{The Split Points}
Consider the split point $P_{v}^{(k)}$ that split up the ordered node embedding into two blocks, with the left block containing neurons indexed between $[0, P_{v}^{(k)}-1]$ and the right block $[P_{v}^{(k)}, D]$. We can represent the split by a $D$-dimensional gating vector $g_v^{(k)}$, with its left $P_{v}^{(k)}$ entries are $1$s and the rest are $0$s. This gating vector could then be used in the \textit{combine} stage to control the integration between the ego representation as input to the $k$-th layer $h_{v}^{(k-1)}$ and the aggregated context $m_{v}^{(k)}$.
\begin{equation}
h_{v}^{(k)}=g_{v}^{(k)} \circ h_{v}^{(k-1)} + \left(1-g_{v}^{(k)}\right) \circ m_{v}^{(k)}
\end{equation}
Where $\circ$ denotes element-wise multiplication. Intuitively, for the left $P_{v}^{(k)}$ neurons in $h_{v}^{(k)}$, since their corresponding gating vectors in $g_{v}^{(k)}$ are $1$s, they suppress the newly aggregated context $m_{v}^{(k)}$ and only inherit the last layer's output $h_{v}^{(k-1)}$. And vice versa for the right block of neurons with gating vectors of $0$.

Ideally, we would prefer a crispy split, with a clear boundary at $P_{v}^{(k)}$ and binary gating vectors in $g_{v}^{(k)}$. However, this will result in discretized operations, rendering the model non-differentiable. In this work, we resort to ``soften'' the gates by predicting the expectations of them, thus keeping the whole model differentiable. The expectation vector $\hat{g}_{v}^{(k)}$ is presented as the cumulative sum of the probability of each position in the node embedding being the split point $P_{v}^{(k)}$.\footnote{We refer the readers to Appendix \ref{proof} for a proof of $\hat{g}_{v}^{(k)}$ being the expectation of the gating vector $g_{v}^{(k)}$.} Concretely, $\hat{g}_{v}^{(k)}$ is parameterized as
\begin{equation}
\label{eq.ordered_gating}
\hat{g}_{v}^{(k)} = \operatorname{cumax}_{\leftarrow}\left( f_{\xi}^{(k)} \left( h_{v}^{(k-1)}, m_{v}^{(k)} \right) \right) = \operatorname{cumax}_{\leftarrow}\left( W^{(k)} \left[ h_{v}^{(k-1)}; m_{v}^{(k)} \right]+b^{(k)} \right)
\end{equation}
where $\operatorname{cumax}_{\leftarrow}(\cdot)=\operatorname{cumsum}_{\leftarrow}(\operatorname{softmax}(\ldots))$ is the composite operator of cumulative sum and softmax, with the subscript $\leftarrow$ indicating that the direction of cumulation is from right to left. $f_{\xi}^{(k)}(\cdot)$ denotes a function fusing the two vectors $h_{v}^{(k-1)}$ and $m_{v}^{(k)}$. Here we simply put it as a linear projection with bias $b^{(k)}$. The $\left[ h_{v}^{(k-1)}; m_{v}^{(k)} \right]$ indicates the concatenation of the two vectors $h_{v}^{(k-1)}$ and $m_{v}^{(k)}$.

\subsection{The Differentiable $\operatorname{OR}$ Operator}
In the above section we provide a differentiable way to predict the position of the split points, as well as corresponding gates that implement the split operation. However, there's no guarantee that the predicted split points will comply with the restrictions in Equation \ref{relativemagnitude2}. Not ensuring Equation \ref{relativemagnitude2} would break the alignment between the rooted-trees and node embeddings.

To ensure such a relative magnitude, we conduct a bitwise $\operatorname{OR}$ operator between the newly computed gating vector $\hat{g}_{v}^{(k)}$ and its predecessor. We change the notation to $\tilde{g}_{v}^{(k)}$ to discriminate it from the originally computed gating vector $\hat{g}_{v}^{(k)}$. Again, the $\operatorname{OR}$ operator is not differentiable in nature. We instead implement the softened version of it by
\begin{equation}
\label{softor}
\tilde{g}_{v}^{(k)} = \operatorname{SOFTOR}(\tilde{g}_{v}^{(k-1)}, \hat{g}_{v}^{(k)}) = \tilde{g}_{v}^{(k-1)} + (1-\tilde{g}_{v}^{(k-1)}) \circ \hat{g}_{v}^{(k)}
\end{equation}
We call it the $\operatorname{SOFTOR}(\cdot)$ operator.

\begin{figure
\centering
\includegraphics[width=\linewidth]{rooted_tree_cases.pdf}
\caption{Some rooted-tree cases, each rooted-tree is constructed from root node by message passing, node's color indicate it's label. Case 1 and Case 2 show two cases under heterophily; Case 3 is under homophily, node 1 \& 2 are at the center and the border of a community, respectively.}
\label{fig.showing_cases}
\end{figure}

\subsection{Putting it All Together}
\label{puttogether}
With the above discussion, we can get the updating rule for the $k$-th GNN layer by putting everything together, leading to the \textit{Ordered GNN} model:

\begin{equation}
\begin{aligned}
m_{v}^{(k)}&=\operatorname{MEAN}\left(\left\{h_{u}^{(k-1)}: u \in \mathcal{N}(v)\right\}\right) \\
\hat{g}_{v}^{(k)}&=\operatorname{cumax}_{\leftarrow}\left( f_{\xi}^{(k)} \left( h_{v}^{(k-1)}, m_{v}^{(k)} \right) \right) \\
\tilde{g}_{v}^{(k)}&=\operatorname{SOFTOR}\left(\tilde{g}_{v}^{(k-1)}, \hat{g}_{v}^{(k)}\right)\\
h_{v}^{(k)}&=\tilde{g}_{v}^{(k)} \circ h_{v}^{(k-1)} + \left(1-\tilde{g}_{v}^{(k)}\right) \circ m_{v}^{(k)}
\end{aligned}
\end{equation}

The first equation corresponds to the \textit{aggregate} stage, while the latter three correspond to the \textit{combine} stage. In practice, we can decrease the parameters in the gating network (parameters in the second equation) by applying a chunking trick, i.e., we set $D_{m}=D/C$ as the size of the gating vector, where $D$ is the dimension of the GNN's hidden state, and $C$ is the chunking size factor. In this setting, each gate controls $C$ neurons instead, thus drastically reducing the number of gates to be predicted.

In order to explain the practical benefits our model provided for heterophily and over-smoothing problems, we illustrate three cases in the Figure \ref{fig.showing_cases} and analyse them with the insights from related theoretical works. \citet{luan2021heterophily,ma2021homophily} notice that heterophily is not always harmful for GNNs, and \citet{ma2021homophily} points out that a simple GCN can work well if similar neighbor patterns are shared across nodes with the same label, they called it ``good heterophily''.
In Case 1, two root nodes share the same first-order neighborhood pattern, one GNN layer can capture this ``good heterophily'' and obtain good performance, however, a multi-layer GNN will mix in the second-order neighbor features that have different pattern, the captured ``good heterophily'' would be destroyed. Our model can protect the captured ``good heterophily'' and reject that ``bad heteropihly'' by learning same split points for layer 1 \& 2.
In Case 2, two root nodes' first-order neighbor patterns differ, but they share the same second-order neighborhood pattern. Since an ordinary GNN can only capture information \textit{within $k$-hops}, they can not capture this ``good heterophily'' \textit{at second-order}, whereas our model could do so with the help of split points prediction. Furthermore, we can filter out the first-order ``bad heterophily''.

Some works \citep{li2018deeper,chen2020measuring,ma2021homophily} show that neighbor feature smoothing is the key reason why GNN works, but too much smoothing lead to the mixing of nodes from different classes. Case 3 shows that while a GNN allows node 1 to smooth neighbor features within it's class, it also mixing node features belong to different classes for node 2. With the fine-grained ordered gating mechanism, our model can allow node 1 to aggregate similar features while filtering out information beyond node 2's own, thus helpful for homophily and over-smoothing.

The analysis of these cases shows that it is critical to predict a set of \textit{flexible} split points. By analyzing the convergence rate of the gates, we prove that our model only constrain the relative position of the split points, the delta between split points is very flexible, i.e. the number of neurons assigned to encode the information at different orders is flexible. We provide the proof in Appendix \ref{proof_conv_speed}.

\section{Experiments}
\label{part.Experiments}
In this section, we test our model on 11 datasets. We mainly evaluate our model in 3 aspects: performances on homophily data, on heterophily data, and the robustness towards over-smoothing. We further use 2 larger datasets to see if our model could scale well w.r.t. the size of the dataset, these two datasets are from large-scale homophily and heterophily benchmarks, respectively.

\begin{wraptable}{r}{0.6\textwidth}
\vspace{-0.5cm}
\caption{Dataset statistics.}
\label{table.dataset_stats}
\centering
\vspace{0.1cm}
\resizebox{0.6\textwidth}{!}{
\begin{tabular}{c|ccccc}
\toprule
Dataset &  Classes &  Nodes &   Edges &  Features & $\mathcal{H}(G)$ \\
\midrule
Cora &        7 &   2,708 &    5,429 &      1,433 & 0.81\\
PubMed &        3 &  19,717 &   44,338 &       500 & 0.80\\
CiteSeer &        6 &   3,327 &    4,732 &      3,703 & 0.74\\
Cornell &        5 &    183 &     280 &      1,703 & 0.30\\
Chameleon &        5 &   2,277 &   31,421 &       2,325 & 0.23\\
Squirrel &        5 &   5,201 &   198,493 &       2,089 & 0.22\\
Actor &        5 &   7,600 &   26,752 &       931 & 0.22\\
Texas &        5 &    183 &     295 &      1,703 & 0.21\\
Wisconsin &        5 &    251 &     466 &      1,703 & 0.11\\
ogbn-arxiv &       40 & 169,343 & 1,166,243 &       128 & 0.66\\
arXiv-year &        5 & 169,343 & 1,166,243 &       128 & 0.22\\
\bottomrule
\end{tabular}
}
\vspace{-0.2cm}
\end{wraptable}

\paragraph{Datasets \& Experiment Setup}
\label{exp_setup}
For homophily data, we use three citation network datasets, i.e., Cora, CiteSeer, and PubMed \citep{sen2008collective}. While six heterophilous web network datasets \citep{pei2020geom,musae}, i.e., Actor, Texas, Cornell, Wisconsin, Squirrel and Chameleon, are used to evaluate our model on the heterophily setting. In citation networks, nodes and edges correspond to documents and citations, respectively. Node features are bag-of-words representation of the documents. In web networks, nodes and edges represent web pages and hyper-links, with node features being bag-of-words representation of web pages. For the over-smoothing problem, following \citet{rong2019dropedge}, we use Cora, CiteSeer, and PubMed as testbeds. For the two larger datasets, we select ogbn-arxiv from \citet{hu2020open} and arXiv-year from \citet{lim2021new}. These two dataset have same nodes and edges that represent papers and citations. In ogbn-arxiv, node labels represent the papers' academic topics, which tend to be a homophily setting, while in arXiv-year, the node labels represent the publishing year, thus exhibiting more heterophily label patterns. We include the edge homophily score $\mathcal{H}(G)$ \citep{yan2021two} for each dataset, which represents the homophily level, together with the statistics of the datasets summarized in Table \ref{table.dataset_stats}.
We use the Adam optimizer \citep{kingma2014adam}, and apply dropout and $L_2$ regularization for $f_{\theta}$ and $f_{\xi}^{(k)}$. We also insert a LayerNorm \citep{ba2016layer} layer after every other layers. We set the hidden dimension as 256 and chunk size factor as 4. For each dataset, unless otherwise noted, a 8 layer model is used (which is already deeper than most popular GNN models).
We perform a grid search to tune the hyper-parameters for all the models.
More details are listed in Appendix \ref{exp_detail}.

\subsection{Homophily and Heterophily}
\label{setting_homo_hete}
We evaluate our model's performance on aforementioned 9 datasets, following the same data splits as \citet{pei2020geom}. Results are shown in Table \ref{table.results.hete_homo}. We report the mean classification accuracy with the standard deviation on the test nodes over 10 random data splits. For baselines, we include classic GNN models such as vanilla GCN \citep{kipf2016semi}, GAT\citep{velivckovic2017graph}, and GraphSAGE\citep{hamilton2017inductive}; state-of-the-art methods that are specifically designed for heterophily: H2GCN \citep{zhu2020beyond}, GPRGNN \citep{chien2020adaptive}, and GGCN \citep{yan2021two}; state-of-the-art models on three homophily citation network datasets: GCNII \citep{chen2020simple}, Geom-GCN \citep{pei2020geom}. We also include models that utilize the information from different hops: MixHop \citep{abu2019mixhop}, JK-Net \citep{xu2018representation}; and finally a 2-layer MLP. For the baselines that have multiple architectural variants, we choose the best one for each dataset and denote the model with an asteriod ``*''. All the baseline results are from \citet{yan2021two}, except for those on H2GCN, JK-Net and MixHop are from \citet{zhu2020beyond}. All of them were evaluated on the same set of data splits as ours.

\begin{table
\caption{Experimental results on homophily and heterophily datasets, for each dataset, we bold the model with the best performance.}
\label{table.results.hete_homo}
\centering
\resizebox{\textwidth}{!}
{
\begin{tabular}{c|cccccc|ccc}
\toprule
& Texas & Wisconsin & Actor & Squirrel & Chameleon & Cornell & CiteSeer & PubMed & Cora \\
\midrule
\textbf{Ordered GNN} & $\textbf{86.22}{\scriptstyle\pm4.12}$ & $\textbf{88.04}{\scriptstyle\pm3.63}$ & $\textbf{37.99}{\scriptstyle\pm1.00}$ & $\textbf{62.44}{\scriptstyle\pm1.96}$ & $\textbf{72.28}{\scriptstyle\pm2.29}$ & $\textbf{87.03}{\scriptstyle\pm4.73}$ & $77.31{\scriptstyle\pm1.73}$ & $\textbf{90.15}{\scriptstyle\pm0.38}$ & $\textbf{88.37}{\scriptstyle\pm0.75}$ \\
\midrule
H2GCN* & $84.86{\scriptstyle\pm7.23}$ & $87.65{\scriptstyle\pm4.98}$ & $35.70{\scriptstyle\pm1.00}$ & $36.48{\scriptstyle\pm1.86}$ & $60.11{\scriptstyle\pm2.15}$ & $82.70{\scriptstyle\pm5.28}$ &  $77.11{\scriptstyle\pm1.57}$ & $89.49{\scriptstyle\pm0.38}$ & $87.87{\scriptstyle\pm1.20}$ \\
GPRGNN & $78.38{\scriptstyle\pm4.36}$ & $82.94{\scriptstyle\pm4.21}$ & $34.63{\scriptstyle\pm1.22}$ & $31.61{\scriptstyle\pm1.24}$ & $46.58{\scriptstyle\pm1.71}$ & $80.27{\scriptstyle\pm8.11}$ &  $77.13{\scriptstyle\pm1.67}$ & $87.54{\scriptstyle\pm0.38}$ & $87.95{\scriptstyle\pm1.18}$ \\
GGCN & $84.86{\scriptstyle\pm4.55}$ & $86.86{\scriptstyle\pm3.29}$ &
$37.54{\scriptstyle\pm1.56}$ & $55.17{\scriptstyle\pm1.58}$ & $71.14{\scriptstyle\pm1.84}$ & $85.68{\scriptstyle\pm6.63}$ &  $77.14{\scriptstyle\pm1.45}$ & $89.15{\scriptstyle\pm0.37}$ & $87.95{\scriptstyle\pm1.05}$ \\
\midrule
GCNII* & $77.57{\scriptstyle\pm3.83}$ & $80.39{\scriptstyle\pm3.4}$ & $37.44{\scriptstyle\pm1.30}$ & $38.47{\scriptstyle\pm1.58}$ & $63.86{\scriptstyle\pm3.04}$ & $77.86{\scriptstyle\pm3.79}$ & $77.33{\scriptstyle\pm1.48}$ &
$\textbf{90.15}{\scriptstyle\pm0.43}$ &
$\textbf{88.37}{\scriptstyle\pm1.25}$ \\
Geom-GCN* & $66.76{\scriptstyle\pm2.72}$ & $64.51{\scriptstyle\pm3.66}$ & $31.59{\scriptstyle\pm1.15}$ & $38.15{\scriptstyle\pm0.92}$ & $60.00{\scriptstyle\pm2.81}$ & $60.54{\scriptstyle\pm3.67}$ &  $\textbf{78.02}{\scriptstyle\pm1.15}$ & $89.95{\scriptstyle\pm0.47}$ & $85.35{\scriptstyle\pm1.57}$ \\
\midrule
MixHop & $77.84{\scriptstyle\pm7.73}$ & $75.88{\scriptstyle\pm4.90}$ & $32.22{\scriptstyle\pm2.34}$ & $43.80{\scriptstyle\pm1.48}$ & $60.50{\scriptstyle\pm2.53}$ & $73.51{\scriptstyle\pm6.34}$ & $76.26{\scriptstyle\pm1.33}$ & $85.31{\scriptstyle\pm0.61}$ & $87.61{\scriptstyle\pm0.85}$ \\
JK-Net* & $83.78{\scriptstyle\pm2.21}$ & $82.55{\scriptstyle\pm4.57}$ & $35.14{\scriptstyle\pm1.37}$ & $45.03{\scriptstyle\pm1.73}$ & $63.79{\scriptstyle\pm2.27}$ & $75.68{\scriptstyle\pm4.03}$ & $76.05{\scriptstyle\pm1.37}$ & $88.41{\scriptstyle\pm0.45}$ & $85.96{\scriptstyle\pm0.83}$ \\
\midrule
GCN & $55.14{\scriptstyle\pm5.16}$ & $51.76{\scriptstyle\pm3.06}$ &  $27.32{\scriptstyle\pm1.10}$ & $53.43{\scriptstyle\pm2.01}$ & $64.82{\scriptstyle\pm2.24}$ &
$60.54{\scriptstyle\pm5.30}$ &  $76.50{\scriptstyle\pm1.36}$ & $88.42{\scriptstyle\pm0.50}$ & $86.98{\scriptstyle\pm1.27}$ \\
GAT & $52.16{\scriptstyle\pm6.63}$ & $49.41{\scriptstyle\pm4.09}$ & $27.44{\scriptstyle\pm0.89}$ & $40.72{\scriptstyle\pm1.55}$ & $60.26{\scriptstyle\pm2.5}$ & $61.89{\scriptstyle\pm5.05}$ &  $76.55{\scriptstyle\pm1.23}$ & $86.33{\scriptstyle\pm0.48}$ & $87.30{\scriptstyle\pm1.10}$ \\
GraphSAGE & $82.43{\scriptstyle\pm6.14}$ & $81.18{\scriptstyle\pm5.56}$ & $34.23{\scriptstyle\pm0.99}$ & $41.61{\scriptstyle\pm0.74}$ & $58.73{\scriptstyle\pm1.68}$ & $75.95{\scriptstyle\pm5.01}$ &  $76.04{\scriptstyle\pm1.30}$ & $88.45{\scriptstyle\pm0.50}$ & $86.90{\scriptstyle\pm1.04}$ \\
\midrule
MLP & $80.81{\scriptstyle\pm4.75}$ & $85.29{\scriptstyle\pm3.31}$ & $36.53{\scriptstyle\pm0.70}$ & $28.77{\scriptstyle\pm1.56}$ & $46.21{\scriptstyle\pm2.99}$ & $81.89{\scriptstyle\pm6.40}$ &  $74.02{\scriptstyle\pm1.90}$ & $87.16{\scriptstyle\pm0.37}$ & $75.69{\scriptstyle\pm2.00}$ \\
\bottomrule
\end{tabular}
}
\end{table}

Our Ordered GNN model consistently achieves SOTA both on homophilous and heterophilous datasets. On homophilous datasets, we match the performance of the SOTA model GCNII. However, GCNII performs worse on heterophilous datasets, leaving a gap of more than 20 percents on Squirrel dataset. The amazing performance on this dataset might related to the ``over-squashing'' problem \citep{alon2020bottleneck,topping2021understanding} (see Appendix \ref{Impact}). We also achieve SOTA on the heterophilous datasets and consistently outperform baselines across all datasets. The superior performances on both homophily and heterophily data show the generality of our model.

\begin{figure
\centering
\includegraphics[width=0.9\linewidth]{results_over_smoothing}
\caption{Node classification accuracy results for over-smoothing problem with various depths.}
\label{fig.results_over_smoothing}
\end{figure}

\subsection{Over-smoothing}
\label{setting_smooth}
We test our model's performance with depth 2, 4, 8, 16, 32 to evaluate the ability in alleviating the over-smoothing problem, we include following models for comparison: classical GNN models GCN \citep{kipf2016semi}, GAT \citep{velivckovic2017graph} and GraphSAGE \citep{hamilton2017inductive}, models that modify the architecture to alleviate over-smoothing: ResGCN \citep{li2019deepgcns}, JK-Net \citep{xu2018representation}, and IncepGCN \citep{rong2019dropedge}. We also include three state-of-the-art methods ResGCN+Drop, JK-Net+Drop, and IncepGCN+Drop that is equipped with the advanced training trick DropEdge \citep{rong2019dropedge} which is developed for alleviating over-smoothing. We reuse the results from \citet{rong2019dropedge} for all baselines except GAT. For GAT, we conduct a separate grid search for each depth to tune it's hyperparameters, with the search space the same to ours for fair comparison. Detailed configurations can be found in Appendix \ref{exp_detail}. We follow the same data split from \citet{rong2019dropedge} as used in baselines.

Performances are shown in Figure \ref{fig.results_over_smoothing}. Compared with the baseline models, our model consistently achieve the best performance across all datasets and all layers. Furthermore, for each dataset, our model's performance doesn't drop significantly as some baselines. Note that this is achieved without using training tricks such as DropEdge \citep{rong2019dropedge}. This shows that our model can well alleviate the over-smoothing problem.

\begin{figure}
\centering
\subfigure[Homophily]{
\label{fig.gate_values.Homophily}
\begin{minipage}[b]{0.26\linewidth}
\label{fig.gate_values.CiteSeer_geom_8}
\includegraphics[width=1\linewidth]{CiteSeer_geom_8.png}\vspace{4pt}
\label{fig.gate_values.PubMed_8}
\includegraphics[width=1\linewidth]{PubMed_geom_8.png}\vspace{4pt}
\end{minipage}
}
\subfigure[Heterophily]{
\label{fig.gate_values.Heterophily}
\begin{minipage}[b]{0.26\linewidth}
\label{fig.gate_values.Actor_8}
\includegraphics[width=1\linewidth]{Actor_8.png}\vspace{4pt}
\label{fig.gate_values.Texas_8}
\includegraphics[width=1\linewidth]{Texas_8.png}\vspace{4pt}
\end{minipage}
}
\subfigure[Over-smoothing]{
\label{fig.gate_values.Over-smoothing}
  \begin{minipage}[b]{0.26\linewidth}
\label{fig.gate_values.Cora_32}
\includegraphics[width=1\linewidth]{Cora_full_32.png}\vspace{4pt}
\end{minipage}
}
\caption{Visualization of gating vectors across all layers in our models on (a) homophily, (b) heterophily and (c) over-smoothing settings. Vertical and horizontal axes are gating vectors and chunk channels, respectively. Each box in the boxplot refect the distribution of the value of the corresponding gate across the whole dataset.}
\label{fig.gate_values}
\end{figure}

\subsection{Visualizing the Gating Vectors}
We checkout the gating vectors across all layers in our model, plotted in Figure \ref{fig.gate_values}. We can find that the gates automatically learns to behave differently for homophilous and heterophilous data, and learns to suppress message passing when needed to prevent over-smoothing.

On homophilous data the split points distributes more uniformly at first layer, with a gradually increasing trend of gating vectors. As the layer gets deeper, gates start to saturate. After layer 4, most of the positions are saturated, leaving the last very few positions open for update. This coinsides with the fact that 4 layers are the depth when many GNNs starts to show the over-smoothing phenomenon.

On heterophilous data, a large block of gates become saturated even in the first layer, emphasizes the importance of a node's ego information by suppressing most of the incoming messages. This observation is consistent with the amazing performance of a simple MLP model when deal with heterophily (see Table \ref{table.results.hete_homo}). Moreover, on Texas dataset the saturated gates keep the same for layer 1 and 2, however enlarge on layer 3, reflect the flexibility of our gating mechanism and the ability to capture ``good heterophily'' for high-order neighbors as we described in Section \ref{puttogether}.

We also plot the gating vectors of our model on Cora dataset in Figure \ref{fig.gate_values.Over-smoothing}. We can see that almost all the gates are closed when layers become very deep. This phenomenon reveals Ordered GNN's superior performance on deep models: it manages to keep node representations distinguishable by learning to reject incoming messages from farther-away neighbors.

\subsection{Scaling to Larger Datasets}
\label{setting_large}
The larger ogbn-arxiv and arXiv-year datasets contain much more nodes, and present a complex linking pattern. We compared our model with these models: baseline graph neural models SGC \citep{wu2019simplifying}, GCN \citep{kipf2016semi}, label propagation based methods LP \citep{zhou2003learning}, C\&S \citep{huang2020combining}, strong baselines: JK-Net \citep{xu2018representation}, APPNP \citep{klicpera2018predict}, and state-of-the art model UniMP \citep{shi2020masked},
GCNII \citep{chen2020simple} for ogbn-arxiv, and LINK \citep{lim2021new}, H2GCN \citep{zhu2020beyond} for arXiv-year. Features, links and labels are inherited from \citet{hu2020open}, all the models follow the same data split from \citet{hu2020open} and \citet{lim2021new}. The baselines results are from \citet{lim2021new}, \citet{li2020deepergcn} and \citet{huang2020combining}.

For homophily dataset ogbn-arxiv, our model beats all the baseline methods except UniMP. This is understandable since UniMP based on the graph Transformer much larger model size, and uses the extra node label information. For heterophily dataset, our method outperform the SOTA models. The performance on these two large dataset demonstrate the scalability and expressivity of our model to capture the complex patterns in large networks.

\begin{table}
\caption{Summary of classification accuracy results on larger benchmarks, we bold the model with the best performance.}
\label{table.results.larger_benchmark}
\centering
\resizebox{\textwidth}{!}
{
\begin{tabular}{l|ll|ll|ll|ll|l}
\toprule
& GCN & SGC & LP & C\&S & JK-Net & APPNP & GCNII & UniMP & \textbf{Ordered GNN} \\
ogbn-arxiv & $71.74{\scriptstyle\pm0.29}$ & $69.39$ & $68.50$ & $72.62$ & $72.19{\scriptstyle\pm0.21}$ & $66.38$ & $72.74{\scriptstyle\pm0.16}$ & $\textbf{73.11}{\scriptstyle\pm0.21}$ & $72.81{\scriptstyle\pm0.21}$ \\
\midrule
& GCN & SGC & LP & C\&S & JK-Net & APPNP & LINK & H2GCN & \textbf{Ordered GNN} \\
arXiv-year & $46.02{\scriptstyle\pm0.26}$ & $32.83{\scriptstyle\pm0.13}$ & $46.07{\scriptstyle\pm0.15}$ & $42.17{\scriptstyle\pm0.27}$ & $46.28{\scriptstyle\pm0.29}$ & $38.15{\scriptstyle\pm0.26}$ & $53.97{\scriptstyle\pm0.18}$ & $49.09{\scriptstyle\pm0.10}$ & $\textbf{54.49}{\scriptstyle\pm0.29}$ \\
\bottomrule
\end{tabular}
}
\end{table}

\subsection{Ablation Study}

We perform an ablation study to evaluate our proposed gating mechanism, Table \ref{table.results.ablation} shows the results. ``Bare GNN'' denotes a simple GNN model with a mean aggregator; ``+simple gating'' denotes adding a gating module where each gate is predicted independently by the gating network, this is equivalent to the Gated GNN \citep{li2015gated} with a mean aggregator; ``++ordered gating'' denotes an ordered gating module as we introduced in Equation \ref{eq.ordered_gating}; ``+++SOFTOR'' indicates our proposed model Ordered GNN. We test on three settings: homophily, heterophily, and over-smoothing. We use CiteSeer for both homophily and over-smoothing settings and Wisconsin for heterophily setting. We set model layers as 8 for the first two settings and 32 for the last one.

\begin{wraptable}{r}{0.55\textwidth}
\vspace{-0.55cm}
\caption{Ablation study}
\label{table.results.ablation}
\centering
\vspace{0.1cm}
\resizebox{0.55\textwidth}{!}{
\begin{tabular}{l|ccc}
\toprule
Model & Homophily & Heterophily & Over-smoothing \\
\midrule
Bare GNN & $73.79{\scriptstyle\pm1.57}$ & $53.92{\scriptstyle\pm5.79}$ & $74.70$ \\
\quad +simple gating & $75.78{\scriptstyle\pm1.49}$ & $86.27{\scriptstyle\pm5.23}$ & $79.30$ \\
\quad \quad ++ordered gating & $76.45{\scriptstyle\pm1.35}$ & $85.69{\scriptstyle\pm4.90}$ & $79.40$ \\
\quad \quad \quad +++SOFTOR & $77.31{\scriptstyle\pm1.73}$ & $88.04{\scriptstyle\pm3.63}$ & $80.80$ \\
\bottomrule
\end{tabular}
}
\vspace{-0.25cm}
\end{wraptable}

As we can see, the ``Bare GNN'' model performs poorly in all three settings, especially in the heterophily setting, revealing the challenge of these settings; the ``+simple gating'' model and the ``++ordered gating'' model perform closely, both providing performance improvements, and let the model comparable to strong baselines; with the help of our proposed differentiable OR module, the ``+++SOFTOR'' model obtains SOTA performance, demonstrating the effectiveness of our proposed methods. Overall, the ablation study verified each component we designed for combine stage of message passing is useful and necessary.

\section{Conclusion}
We propose Ordered GNN, a novel GNN model that aligns the hierarchy of the rooted-tree of a central node with the ordered neurons in its node representation. This alignment implements a more organized message passing mechanism, by ordering the passed information into the node embedding w.r.t. their distance towards the node. This model results to be both more explainable and stronger in performance. It achieves the state-of-the-art on various datasets, ranging across both homophily and heterophily settings. Moreover, we show that the model could effectively prevent over-smoothing, when the model gets very deep. We discuss the connections with related works, the limitation of our model and broader impact on building more powerful GNNs in Appendix \ref{Discussion} due to the space limit.

\section{Ethics Statement}
In this work, we propose a novel GNN model to deal heterophily and over-smoothing problem. Our model integrates inductive bias from graph topology and does not contain explicit assumptions about the labels of nodes on graph, thus avoiding the occurrence of (social) bias. However, the gating mechanism of our model needs to learn from samples with labels, and it is possible to be misleading by targeting corrupted graph structure or features, which can result in negative social impact, because both internet data and social networks can be expressed with graphs. We can develop corresponding ethical guidelines, and constraint the usage of our model to avoid such problems.

\section{Reproducibility Statement}
We provide the introduction for our used datasets and model configuration in Section \ref{exp_setup}, the space of grid search to tune our model and the hardware \& software experimental environment are introduced in Appendix \ref{exp_detail}, we also post every set of hyper-parameters to reproduce our results in Table \ref{table.hyper.hete_homo}, Table \ref{table.hyper.over_smoothing} and Table \ref{table.hyper.larger}, including the baseline model's. The download, load, and pre-process steps are introduced in Appendix \ref{exp_detail}. We've test our model on 3 settings, the setting for heterophily \& homophily are introduced in Section \ref{setting_homo_hete}, the setting for over-smoothing is introduced in Section \ref{setting_smooth}, and the setting for larger benchmarks are introduced in Section \ref{setting_large}. In our released source code, we list the steps needed to reproduce results, and also list all the hyper-parameters with yaml files.

\section*{Acknowledgement}
This work was sponsored by the National Natural Science Foundation of China (NSFC) grant (No.
62106143), and Shanghai Pujiang Program (No. 21PJ1405700).


\section{\label{sec:Introduction } Introduction }

Hyperuniform (HU) systems are characterized by the suppression of large-wavelength density fluctuations owing to a hidden order which is not apparent on short lengthscales~\cite{Torquato2003,TORQUATO20181}. Systems possessing long-range translational and orientational order, e.g., crystalline structures and periodic fields, represent trivial cases of HU structures. More interestingly, a HU character may be present in disordered systems with partial or full suppression of density fluctuations for small wavevectors (similarly to crystalline solids) and isotropic behaviors with respect to spatial directions (as for liquids). These features are well described for point patterns by the structure factor $S(\mathbf{k})$ with $\mathbf{k}$ the wavevector. Ideal HU characters corresponds to have a scaling $S(\mathbf{k})\sim |\mathbf{k}|^\alpha$ for $|\mathbf{k}|\rightarrow 0$ and $\alpha > 0$. Generalizations have been proposed to account for such a property in systems as, e.g., heterogeneous media as well as scalar and vector fields \cite{Torquato2016,Ma2017,TORQUATO20181}. In experimental settings, quantification of the HU character (how close a system is to the ideal HU case) is typically considered through the metrics $H=S(|\mathbf{k}|\rightarrow 0)/\max(S(|\mathbf{k}|))$, with condition $H\lesssim 10^{-2}$ and $H\lesssim 10^{-4}$ commonly denoting nearly and effective HU character, respectively \cite{Torquato2006,Kim2018}.

\begin{figure*}[ht!]

   \includegraphics[width=1\textwidth]{fig1.png}
   \caption{
    \textit{Influence of FIB single scanning mode}. a) High-resolution scanning electron microscopy (SEM) images of a Ge layer on SOI processed by FIB with a single scan strategy using low, medium, and high resolution (respectively 256 × 221 pixels, 512 × 442 pixels and 1024 × 884 pixels) with ion beam fluency in the 10$^{15}$ ions cm$^{-2}$ range. b) SEM image of a Ge layer on SOI processed by FIB with a multiple scan strategy. In this condition, high and low resolution exhibit the same final morphology with a nearly-HU character. All these patterns above were fabricated on a stack as follows: the structure is composed of a top amorphous germanium layer on a SOI substrate with a top thin silicon layer of 7~nm and a 2000~nm thick SiO$_{2}$ for the BOX. (c) Spectral density $\widetilde{\chi}(\mathbf{k})$ for the samples illustrated in panels (a) and (b). The white-to-black color scale is saturated at $0.75\max(\widetilde{\chi}(\mathbf{k}))$. $K$ corresponds to the wavevector length at which the $\max(\widetilde{\chi}(\mathbf{k}))$ is obtained.}
    \label{fig:figure1}
\end{figure*}

Materials and architectures possessing a nearly, effective or ideal HU character have gained a lot of interest in recent years. In addition to fundamental physical aspects~\cite{GAB2002,TORQUATO20181}, their importance
relies in the novel possibilities in light~\cite{FRO2017,ZHO2016,yu2020,FLO2009,MAN2013a,FRO2016,wu2017} and carrier management with topologically-protected electronic states~\cite{MIT2018} as well as vortices in superconductors \cite{rumi2019}. A crucial aspect consists of finding a proper fabrication method that allows the realization of a HU character, e.g., small enough $H$, at wavelengths of interest for the targeted applications. Design protocols for the realization of HU patterns have been developed \cite{FLO2009}. They have usually been combined with top-down approaches. Although powerful, see e.g. \cite{CAS2017,granchi2022near}, such approaches usually show poor scalability. This becomes crucial for patterns meant for large-area applications with features on the nanometer scale, such as metals and nanostructured semiconductor films for photonics applications.

Devising scalable approaches for the fabrication of HU materials featuring nanostructured architectures represents then an important goal. Self-assembly paradigms or hybrid bottom-up/top-down approaches have been exploited for decades to obtain nano- and micro-structures from film-like settings \cite{Medeiros-Ribeiro1998,Ye2011,Naffouti2017}. Self-assembly based on spinodal-like, solid-state dewetting has also recently proved successful for the fabrication of effectively HU, SiGe-based nano-architectures with a footprint of about 100~nm~\cite{salvalaglio2020}. Such an approach can be scaled-up to large surfaces~\cite{benali2020}. Moreover, soft nano-imprint lithography of sol-gel materials (such as silica and titania) was also exploited to fabricate nearly-HU metasurfaces using dewetted micro-and nano-architectures as hard masters \cite{chehadi2021}.

A remarkable example of self-assembly in crystalline and amorphous Ge semiconductor thin-films consists of the impact of ion beams, widely investigated both from experimental~\cite{bischoff2011,kolibal2011low,rudawski2013,bottger2013} and theoretical point of views~\cite{bradley1988,nord2002,hartmann2009}.
Formation of periodic, ordered~\cite{wei2009,fritzsche2012} and disordered nano-architectures with this approach has been largely addressed so far, and we thus refer the reader to the existing literature~\cite{bradley1988,cuerno1995,castro2005,kolibal2011low}. Central aspects of this approach are its versatility, as the features can be tuned by varying the ion beam parameters (e.g., ion species, energy, and dose ~\cite{bottger2013,hooda2016}, incidence angle~\cite{mollick2014}, sample temperature~\cite{bischoff2011} and stochiometry~\cite{alkhaldi2016}) and processable surface that can reach 200 mm wafers. These structures can also extend over several tens of micrometers in depth, below the surface, resulting in actual three-dimensional, porous nanoarchitectures with an average cavity radius as small as a few tens of nm~\cite{romano2012,alkhaldi2016}.

Here we report on Ge-based nano-structures obtained by the impact of a high-energy, Ga$^{+}$ ion beam on Ge layers resulting in ordered and disordered nano-architectures and characterize their morphological features. In Sect.~\ref{sec:experiments} we show the formation of strongly anisotropic and wavy structures as well as isotropic disordered ones, obtained by adjusting ion current, supplied dose, and scanning parameters. We focus then on the isotropic disordered pattern, being the main target of this study, and we discuss its features in Sect.~\ref{sec:experiments}. In particular, the nearly-HU character of the disordered structures is assessed through the analysis of their spectral density that accounts for the formation of Ge-based nearly-HU nano-architectures having a record footprint of about 40~nm. The analysis via Minkowski functionals further assesses the deviations from a (Gaussian) random field and supports the presence of correlations in the patterns.
A proof of concept for potential applications is then given in Sect.~\ref{sec:application}, showing peculiar reflection properties of the resulting nanostructured film. For the sake of readibility, we report in the main text the main evidence and results, while including additional materials in the Supporting Information.
Conclusions are summarized in Sect.~\ref{sec:conclusions}

\section{Experiments}
\label{sec:experiments}

Experiments were performed on a 50-nm thick amorphous germanium layer deposited on a silicon-on-insulator (SOI) substrate composed of a top silicon layer of 7~nm on a 2000~nm thick SiO$_{2}$ layer. After chemical cleaning (5 seconds in 10$\%$ HF solution in H$_{2}$O in N$_{2}$ atmosphere), amorphous Ge is deposited by molecular beam epitaxy (Riber R32, in ultra-high vacuum of $\sim$10$^{-10}$ torr, at room temperature and growth rate of 10~\AA/min).

Although very similar structures can be obtained on crystalline bulk Ge (Supplementary Information Figure 5), we chose to implement them in ultra-thin Ge layers on commercial ultra-thin SOI. This substrate is more suitable for electronic and photonic device integration and is thus more appealing than the bulk counterpart. For the sake of thoroughness, we mention that several other experimental conditions were explored, such as a-Ge 80~nm thick deposited atop SOI (Supplementary Information Figure 5) or custom-made germanium on insulator wafers with 200~nm a-Ge atop 300~nm of SiO$_2$ (not shown). They always provided similar results to those shown here, accounting for the possibility of fabricating these structures in different conditions.

For this study, we used a dual-beam FEI Strata DB235 with Ga$^{+}$ (atomic mass 31 amu) liquid-metal ion source focused ion beam (FIB) and a SFEG source for the scanning electron microscope imaging (SEM). We used ion currents within 15 and 1000~pA at 30~keV. The samples are kept orthogonal to the Ga ion beam. The FIB beam moves horizontally, from the top left part of the window and line-by-line, raster-scans the sample surface, down to the bottom right point.

For single scans (Figure~\ref{fig:figure1} a), the ion beam resolution is set to low (256 $\times$ 221 pixels), medium (512 $\times$ 442 pixels), and high (1024 $\times$ 884). At 50000$\times$ magnification (used for the milling process), the field of view has a side of 6080~nm. Thus, each pixel is respectively separated by 24, 12, and 6~nm for low, medium, and high resolution. The ion current was fixed at 30~pA, corresponding to a beam diameter of about 15~nm. Thus, only low-resolution modes do not result in a beam overlapping from one point to the other. As shown in Figure\ref{fig:figure1} a), different resolution impacts the final morphology with different features.

In analogy with what is found under strong laser illumination\cite{HER1998} and ion beam impact in metals\cite{repetto2012,repetto2015}, several factors are expected to contribute to the formation of ordered and disordered patterns, namely ion-induced melting of the a-Ge layer, curvature-dependent sputtering and/or swelling \cite{yanagisawa2007}. For Ge layers, a relevant role in the pattern formation was attributed to vacancy clustering \cite{rudawski2013}. Temperature and mass redistribution at surfaces are also expected to play a role as supported by observations in similar systems\cite{wei2009} and by the theory of pattern formation by ion bombardment \cite{bradley1988,bradley2010,cuerno1995,castro2005}. A comprehensive investigation and tuning of all these contributions are beyond the scope of this work, which instead focuses on selected conditions leading to interesting and potentially relevant properties, as discussed in the following.

\begin{figure}[t]

   \includegraphics[width=0.9\columnwidth]{fig3.png}
   \caption{
   \textit{Spectral density and assessment of the nearly-HU character}. a) Radial distribution of the normalized spectral density for a porous a-Ge sample averaged over ~50 images of different location in the same sample (see the Supplementary Information Figure 8). The shaded grey areas displays the variability of the spectral density measured over all the 50 repetitions of the pattern. b) Log-log plot for the curve in panel (a) focusing on the region $|\mathbf{k}|/K\leq1$ and showing the decay for $|\mathbf{k}|\rightarrow 0$.
    \label{fig:figure3}
\end{figure}

\begin{figure}[t]

   \includegraphics[width=1\columnwidth]{fig2.png}
   \caption{
   \textit{Disorder and correlations of spontaneous patterns}. a)  Minkowski functionals analysis $m_{1,2}(m_{0})$ for a 30~nm thick Ge layer exposed to multiple scans with a  2.7$\times$10$^{15}$ ions cm$^{-2}$ dose.  b) Same as a) for 8.02$\times$10$^{15}$ ions cm$^{-2}$ dose. Continuous lines represent the trend of $m_{1,2}(m_{0})$ for a Gaussian random field. The insets show SEM micrographs with details of the sample's surface. }
    \label{fig:figure2}
\end{figure}

\section{Pattern analysis}

Global features of the patterns as in Figure~\ref{fig:figure1} can be analyzed by looking at their spectral density $\widetilde{\chi}(\mathbf{k})$. It provides the generalization for the information conveyed by the structure factor $S(\mathbf{k})$ for point patterns \cite{TORQUATO20181}. The SEM images collected for the structures (c.f. Figure~\ref{fig:figure1}) well resolve the region occupied by the solid phase with respect to the surrounding vacuum. Therefore we consider the spectral density $\widetilde{\chi}(\mathbf{k})$ for a two-phase (solid-vacuum) system \cite{Torquato2016}. It is based on an indicator function:
\begin{equation}
\mathcal{I}^{(s)}(\mathbf{r})=
    \begin{cases}
    1, & \mathbf{r}\in \mathcal{V}_s, \\
    0, & \mathbf{r}\notin \mathcal{V}_s,
    \end{cases}
\end{equation}
with $\mathbf{r}\in \Omega$ the spatial coordinates and $\mathcal{V}_s$ the region with grey-scale values above a given threshold. $\widetilde{\chi}(\mathbf{k})$ then corresponds to the Fourier transform of the autocovariance function: \cite{torquato1999exact,Torquato2016}
\begin{equation}
    \chi(\mathbf{r})=\mathcal{S}_2^{(s)}(\mathbf{r})-(\phi^{(s)})^2,
\end{equation}
with
\begin{equation}\label{eq:Ix}
    \mathcal{S}_2^{(s)}(\mathbf{r})=\langle \mathcal{I}^{(s)}(\mathbf{r}')\mathcal{I}^{(s)}(\mathbf{r}'+\mathbf{r}) \rangle,
\end{equation}
the two-points autocorrelation function and $\phi^{(s)}=\langle \mathcal{I}^{(s)}(\mathbf{r}') \rangle$ the volume fraction of the solid phase. The spectral density can be computed directly as: \cite{torquato1999exact}
\begin{equation}
    \widetilde{\chi}(\mathbf{k})=\frac{1}{|\Omega|}|\mathcal{F}[\mathcal{I}^{(\rm{s})}(\mathbf{r})-\phi_i]|^2
\end{equation}
with $\mathcal{F}$ the Fourier transform and $|\Omega|$ the area of the images. The decay of $\widetilde{\chi}(\mathbf{k})$ for $\mathbf{k}\rightarrow 0$ conveys similar information to $S(|\mathbf{k}|)$ \cite{Torquato2016,TORQUATO20181}\footnote{Interestingly, results similar to the one reported in the following can be obtained by sampling the solid phase with randomly distributed points and computed $S(\mathbf{k})$.}.

The spectral density for the patterns obtained in Figure~\ref{fig:figure1} a) is shown in the first three panels of Figure~\ref{fig:figure1} c). The structures exhibit anisotropic arrangements with marked periodicities. The characteristic frequency $K$ (where the peak of $\widetilde{\chi}(\mathbf{k})$ is observed) is reported in the figures and points to characteristic length scales $2\pi/K$. For low resolution, the structures present a higher degree of disorder which is then reduced for larger resolution. The patterns are also found to exhibit slightly different predominant orientations $>$83$^\circ$ tilt with respect to the scan direction, medium resolution 75$^\circ$-83$^\circ$, and low resolution 61$^\circ$-70$^\circ$. As expected for (partially-) ordered structures, the corresponding spectral density show an exclusion zone at small wavevectors.
This analysis shows that, depending on how the ions are supplied, different anisotropic structures can be obtained with tunable disorder.

A behavior similar to the patterns illustrated above in terms of suppression of correlations at large wavelength (with an exclusion zone at small wavevectors), has been obtained through a multiple scans strategy, using an overall dose of $9.1 \times$10$^{15}$ ions cm$^{-2}$(Figure~\ref{fig:figure1} b). However, this time the final structure features a spectral density $\widetilde{\chi}(\mathbf{k})$ with isotropic angular distribution (last panel in Figure~\ref{fig:figure1} c) that corresponds to a disordered and uniform distribution of objects. Exploring the three different resolution configurations always provides the same morphology featuring connected and disordered structures, even when using an overall dose comparable to that of a single scan (see Supplementary Information Figure 7).

Beyond raster scans with a small spot (as shown in Figure\ref{fig:figure1} b)), disordered, porous a-Ge can be obtained also by ion impact on a single large spot (see the Supplementary Information Figure 1-4). This demonstrates that the fabrication process can be scaled to much larger surfaces as it does not require a high lateral resolution. Nonetheless, the possibility of creating these patterns over small and well-defined areas may allow to integrate them within other devices.

\begin{figure}[h!]

   \includegraphics[width=0.9\columnwidth]{fig4.png}
   \caption{
   \textit{Light reflection from porous a-Ge}. a) Reflection from 50~nm a-Ge on 30~nm thick SiO$_{2}$ normalized by the white tungsten lamp used for illumination. Open symbols correspond to flat a-Ge, whereas continuous lines are impacted areas with a multiple scan approach. Spectra are collected in bright field illumination with a 100$\times$ magnification objective lens, numerical aperture 0.9 (collection angle about $\pm$64 degrees, lateral, resolution $\sim$6~$\mu$m. Right inset: from top to bottom, bright-field images of the impacted areas with doses ranging from 0 to 2.3$\times$10$^{16}$ ions cm$^{-2}$ in steps of 5$\times$10$^{15}$ ions cm$^{-2}$. The highlighted areas are selected for spectroscopy. }

    \label{fig:figure4}
\end{figure}

We now further asses the disorder found in the case of uniform structures. We provide a quantitative assessment of the HU character of the uniform and disordered pattern in Figure~\ref{fig:figure1} b) and address their scalability by analyzing 50 separate images in different points in the sample (see a collage in the Supplementary Information Figure 8). In Figure~\ref{fig:figure3}, we report the normalized radial distribution of $\widetilde{\chi}(\mathbf{k})$ averaged over all the analyzed images. Grey areas show the range of variability of the results across the different images. Figure~\ref{fig:figure3}a shows a marked peak at about 0.16~nm$^{-1}$ corresponding to a characteristic length-scale of $\sim 39$~nm ($\pm 4$~nm). We observe a decay of $\widetilde{\chi}(\mathbf{k})$ toward zero, sampled up to the limits of the analyzed window, well approximated by $|\mathbf{k}|^{\alpha}$ with $\alpha\approx 3.5$ close to $|\mathbf{k}|/K=1$ (see log-log plot in Figure~\ref{fig:figure3}b).

An extrapolation of $\widetilde{\chi}(0)$ through fitting allows to estimate a H-metric value of $H\approx2\cdot 10^{-2}$ for the averaged $\widetilde{\chi}(\mathbf{k})$, and within a range $H\in [8 \cdot 10^{-3}-3 \cdot 10^{-2}]$ among the analysed samples. This value is compatible with a nearly-HU character \cite{Torquato2006,Kim2018} similarly to other systems where correlation and HU character emerge spontaneously \cite{wilken2020,xie2013hyperuniformity}.

The features of the disordered structures, beyond simple random fields, can be further accounted for by Minkowski functional analysis describing their topological features~\cite{MIN1989,SCH2014,salvalaglio2020}. We plot Minkowski functionals $m_i(\bar{\rho})=(1/|\Omega|)M_i(\mathcal{B}_{\bar{\rho}})$ %
with
$\bar{\rho}$ a given threshold defining a binary image of thresholded gray-scale images~\cite{Mantz_2008,salvalaglio2020} (Figure~\ref{fig:figure2}).
$m_0(\bar{\rho})$ is the fraction of $|\Omega|$ occupied by the non-zero region in $\mathcal{B}_{\bar{\rho}}$, as $M_0(\mathcal{B}_{\bar{\rho}})$ represents the area of  $\mathcal{B}_{\bar{\rho}}=1$. $m_1(\rho)$ is the average of the boundary length $U$ between the areas where $\mathcal{B}_{\bar{\rho}}=1$ and $\mathcal{B}_{\bar{\rho}}=0$, as $U=2 \pi M_1(\mathcal{B}_{\bar{\rho}})$. $m_2(\rho)$ is the averaged Euler characteristic $\chi$, as $\chi=\pi M_2(\mathcal{B}_{\bar{\rho}})$. With this method, the results are independent of image saturation and contrast~\cite{Mantz_2008,salvalaglio2020}.

$m_{1,2}(m_{0})$ assess the deviation from a random field (solid lines in Figure~\ref{fig:figure3} a) and b)) and eventual non-linearity\cite{bradley1988,alkemade2006} of the underlying pattern formation dynamics~\cite{Mantz_2008}. We observe significant deviations from a Gaussian random field of $m_{1}(m_{0})$ and $m_{2}(m_{0})$ for porous Ge irradiated with 8.0$\times$10$^{15}$ ions cm$^{-2}$ (Figure~\ref{fig:figure3} b)), which increase with increasing the dose.
Other FIB settings (e.g., lower beam energy), as in T. Bottger et al. \cite{bottger2013} result in disordered structures very close to a Gaussian random field (not shown).

\section{Optical properties}
\label{sec:application}

Finally, we address the changes in light reflection from the impacted areas featuring a nearly-HU character. A spectroscopic characterization of light reflection performed in bright-field illumination reveals strong changes in the affected areas' optical response with respect to the flat counterpart (Figure~\ref{fig:figure4}). The 50~nm thick a-Ge layer deposited atop 30~nm thick SiO$_{2}$ shows a drastic reduction in the reflected light intensity when impacted with a dose larger than $\sim$ 5$\times$10$^{15}$ cm$^{-2}$ in single scan mode. A red colorization appears for larger doses (from about 1.3$\times$10$^{16}$  cm$^{-2}$) owing to enhanced absorption at short wavelength, cutting the blue light reflection (up to about 25 times larger than the flat counterpart).
Qualitatively similar results have been obtained with 200~nm thick a-Ge, and with a-Ge deposited on thin and thick silicon on insulator (not shown).

Provided the presence of several layers with different optical constants, the overall response to illumination is determined by etalon effects and by the disordered network of filaments~\cite{galinski2017,galinski2017b}.

\section{Conclusions}
\label{sec:conclusions}
In conclusion, we showed that impacting a thin layer of Ge with an ion beam results in ordered and disordered structures with a nearly-HU character. Their typical size is the smallest reported so far for analogous HU systems~\cite{salvalaglio2020,CAS2017,tavakoli2020}. Provided the possibility to change stoichiometry by Si alloying~\cite{alkhaldi2016}, tune the size of the pores by changing ion dose~\cite{hooda2016,rudawski2013} or Ge film thickness~\cite{romano2012} and also obtain truly 3D structures~\cite{romano2012,rudawski2013}, this fabrication approach has considerable potential for overcoming the limitations of conventional methods for nanoarchitectures with HU character~\cite{soukoulis2011}.

\section*{Acknowledgments}
This research was funded by the EU H2020 FET-OPEN project
NARCISO (No. 828890). J. W. and J.-B. C. acknowledge
the European Research Council (ERC, Grant Agreement
No. 723241). We acknowledge the Nanotecmat platform of the IM2NP institute of Marseilles.

\providecommand{\noopsort}[1]{}\providecommand{\singleletter}[1]{#1}%