# Hierarchical Source-to-Post-Route QoR Prediction in High-Level Synthesis
with GNNs

Mingzhe Gao1, Jieru Zhao1∗, Zhe Lin2∗, Minyi Guo1 1Shanghai Jiao Tong
University, 2Sun Yat-sen University
<EMAIL_ADDRESS><EMAIL_ADDRESS>guo-
<EMAIL_ADDRESS>
###### Abstract

High-level synthesis (HLS) notably speeds up the hardware design process by
avoiding RTL programming. However, the turnaround time of HLS increases
significantly when post-route quality of results (QoR) are considered during
optimization. To tackle this issue, we propose a hierarchical post-route QoR
prediction approach for FPGA HLS, which features: (1) a modeling flow that
directly estimates latency and post-route resource usage from C/C++ programs;
(2) a graph construction method that effectively represents the control and
data flow graph of source code and effects of HLS pragmas; and (3) a
hierarchical GNN training and prediction method capable of capturing the
impact of loop hierarchies. Experimental results show that our method presents
a prediction error of less than 10% for different types of QoR metrics, which
gains tremendous improvement compared with the state-of-the-art GNN methods.
By adopting our proposed methodology, the runtime for design space exploration
in HLS is shortened to tens of minutes and the achieved ADRS is reduced to
6.91% on average. Code and models are available at https://github.com/sjtu-
zhao-lab/hierarchical-gnn-for-hls.

**footnotetext: Jieru Zhao and Zhe Lin are the corresponding authors.

##  I Introduction

High-level synthesis (HLS) enables users to create hardware designs at a high
level of abstraction using programming languages like C/C++, and provides a
series of pragmas, to tune hardware micro-architectures. Due to its superior
productivity, HLS is becoming a popular solution for next-generation hardware
development, especially for FPGAs. However, the turnaround time for HLS design
is significantly increased when post-route quality of results (QoR), such as
latency and resource utilization, are considered during optimization. To
obtain an accurate assessment of post-route QoR metrics, a complete C-to-
bitstream design flow should be invoked each time the HLS code is modified.
This runtime overhead undermines the benefits brought by HLS, especially when
it comes to choosing the solution with optimal QoRs from an extremely large
design space. Therefore, it is essential to provide a fast and accurate QoR
prediction framework at an early design stage like HLS to improve design
efficiency.

Existing prediction methods can be classified into analytical and learning-
based models. Analytical models are constructed based on domain expertise and
achieve high accuracy in latency prediction [1, 2]. However, they are tailored
for a specific subset of pragmas, resulting in limited scalability. Learning-
based models adopt machine learning techniques [3, 4, 5] and graph neural
networks (GNNs) [6, 7, 8] to predict latency and resource usage in HLS. We
compare representative methods in Table I. Zhong et al. [3] extract features
via source code profiling and use the gradient-boosted machine (GBM) to
predict post-HLS metrics that can be obtained from HLS reports after
synthesis. Then they conduct efficient DSE guided by their model. However, the
labels in their dataset, i.e., post-HLS metrics, may deviate significantly
from actual post-route QoR values. Hence, predicting post-HLS metrics may
mislead the DSE process and get sub-optimal designs. Dai et al. [4] and
Pyramid [5] also utilize ML models while predicting post-route QoR metrics
directly. However, their features come from HLS reports and necessitate the
execution of the HLS flow, incurring a higher time cost. Recent advances
introduce GNN models for QoR prediction due to their representation power of
C/C++ programs. GNN-DSE [6] and Ferretti et al. [7] represent source code as
graphs and predict post-HLS metrics, which induce the same issue as [3]. Wu et
al. [8] start from the HLS intermediate representation (IR) and predict post-
route QoRs. They first predict resource types of operations and then
incorporate the predicted resource types as node features to further estimate
post-route resource usage. However, their input graphs are constructed on top
of HLS IR, requiring running the HLS flow to generate related information.
More importantly, the samples in their dataset exhibit relatively simple code
structures, including randomly generated DFGs and loops without pragmas
applied. The missing pragma modeling makes it less suitable for DSE tasks.

Table I: Comparison between representative methods. |  Model |  Input |  Prediction Targets |  Pragmas Involved |  HLS Exec. Free
---|---|---|---|---|---
[3] | GBM | Source code |  Post-HLS | ✓ | ✓
[4] | XGB |  HLS reports |  Post-Route | ✗ | ✗
[6] | GNN | Source code |  Post-HLS | ✓ | ✓
[8] | GNN | HLS IR |  Post-Route | ✗ | ✗
Ours | GNN | Source code |  Post-Route | ✓ | ✓

To summarize, few of the previous methods achieve accurate source-to-post-
route QoR prediction. The problem is even more challenging if the code
structure is complicated and multiple HLS pragmas are considered. Different
from directly predicting QoR with the whole graph in previous methods, our
insight is that we can reserve the structural information of loops/functions
and predict post-route QoR metrics gradually from inner to outer hierarchies.
This will reduce the estimation difficulty and potentially improve the
prediction accuracy.

To this end, we propose a hierarchical prediction method using GNNs, which
distinguishes itself from prior arts with the following key features: 1)
efficiency: our method directly takes C/C++ source code as input and gives a
fast post-route QoR estimation without invoking any EDA design flow, inlcuding
HLS, during inference; 2) pragma embedding: we introduce an effective graph
construction approach that jointly represents the control and data flow graph
of source code and effects of HLS pragmas; 3) hierarchical training and
prediction: our method is capable of capturing the structural information of
code hierarchies and presents the ability to estimate QoR for different
loop/function hierarchies; 4) accuracy: our method presents an average
prediction error of $<$10% for post-route QoR metrics, outperforming the
state-of-the-art GNN methods. With efficient and accurate QoR estimation, the
DSE time is shortened to tens of minutes and the achieved ADRS is reduced to
6.91% on average.

##  II PRELIMINARIES

###  II-A Control and Data Flow Graph (CDFG)

The Control Flow Graph (CFG) is a graphical representation of a program that
describes the possible flows between basic blocks in the program. Likewise,
the Data Flow Graph (DFG) is another graphical representation that depicts the
dependencies between operations. By combining CFG and DFG, we obtain a Control
and Data Flow Graph (CDFG). Using CDFG to represent programs offers several
advantages: 1) CDFG can be generated with LLVM IR, allowing for the earliest
possible graph description of a program; 2) CDFG naturally reflects
hierarchies of loops, making it convenient to annotate loop-based information
at different levels; 3) CDFG contains control and data flow information that
can be used to distinguish whether loops can be parallelized and determine
whether resources can be shared.

###  II-B Graph Neural Network

GNNs process graph-structured data, where objects are nodes and relationships
are edges. GNNs learns node representation with a _message passing_ mechanism:
each node $v$ in the graph $G$ collects the embedding $h^{k}_{N_{v}}$ of its
$k$-hop neighbor ${N_{v}}$ to update its embedding $h^{k}_{v}$ after $k$ graph
convolution layers.

Since input C/C++ programs can be represented as graphs (i.e., CDFGs)
naturally, GNNs become a powerful technique to capture features of each
operation and recognize relationships between neighbor operations. Due to its
representation power, we make full use of GNNs in a hierarchical way for the
post-route QoR prediction.

Figure 1: Framework overview

##  III Methodology

Figure 1 illustrates the overview of our proposed approach. During the
training phase, we deploy the complete C-to-bitstream design flow to obtain
the ground-truth QoR considering various combinations of HLS pragmas. We
extract resource statistics from implementation reports after place and route
(PAR) and latency details from HLS reports. Our approach starts with the input
source code (C/C++) which is transformed into LLVM IR. Then we construct a
graph representation of the input code with different pragma configurations
applied (graph construction). After that, node and edge features are extracted
and annotated in the graph (feature annotation). By bundling the annotated
graphs and their corresponding ground-truth labels under the same pragma
configuration, our dataset can be generated (dataset generation). In the final
step, we train multiple GNN models for accurate post-route QoR prediction,
while incorporating the inherent hierarchical nature of loops (GNN model
training). During the inference phase, accurate predictions of post-route
performance are provided for a given input program, without invoking any HLS
or RTL implementation flow. Based on the fast feedback from the prediction
model, efficient DSE can then be conducted to search for the optimal
configuration of pragmas that achieves the best performance.

###  III-A Graph Construction

The input source code is first transformed to LLVM IR using the Clang front-
end. Then we construct the corresponding graph representation by extending a
CDFG generator, PrograML [9], which captures both control and data flow.
However, PrograML does not target FPGAs, and the generated program graphs need
to be extended to consider the effects of HLS pragmas. Figure 2 illustrates
our graph construction process considering different kinds of pragmas.

####  III-A1 Constructing graph with loop pipelining

In a pipelined loop, the execution of consecutive iterations can be initiated
after a specific interval, i.e., initiation interval (II). This is achieved by
inserting registers into the logic during the synthesis process of HLS tools.
The equivalent graph representation is not influenced or changed at the
compilation stage. Therefore, the graphs constructed for pipelined loops
remain the same as those of original loops without pipelining, as shown in
Fig. 2(a).

####  III-A2 Constructing graph with loop unrolling

The loop unrolling pragma facilitates the concurrent execution of iterations
and results in the creation of multiple replicas of hardware logic. In
accordance with this principle, when constructing graphs for unrolled loops,
the pertinent logic nodes in the unrolled region are replicated and
appropriate edges are established to connect to original predecessors and
successors. This procedure is illustrated in Fig. 2(b). By adding replicas of
logic nodes explicitly, the resulting graph achieves a natural representation
of loop unrolling and significantly alleviates the complexity of prediction.

####  III-A3 Constructing graph with array partitioning

HLS tools provide the array partitioning pragma to increase local memory
bandwidth. Arrays are divided into smaller ones based on user-specified
options. To reflect the impact of memory ports on performance estimation, we
incorporate I/O port nodes to the original CDFG, as highlighted in orange in
Fig. 2. When array partitioning is applied, the corresponding memory port
nodes are split into a certain number of nodes based on partitioning factors.
These newly introduced memory port nodes are connected to load and store nodes
that read from or write to corresponding memory locations. The connections are
determined by partitioning types as well as the memory access pattern of the
code. Specifically, for a $k$-dimensional array, assuming its partition
factors across the $k$ dimensions are $u_{1},u_{2},\cdots,u_{k}$, we add
$\prod_{i=1}^{k}u_{i}$ memory port nodes in the graph. Simultaneously, we
write LLVM passes to analyze the index values of each load and store operation
to determine which memory ports should be connected to them. If the array
index is dynamic or unpredictable, the load/store nodes will connect to all
memory ports.

Figure 2: Procedure of graph construction.

###  III-B Feature Extraction and Annotation

Effective features need to be extracted from the code and annotated to the
generated graph for model training. We categorize the features we have used
into two categories: (1) node features and (2) loop-level features.

####  III-B1 Node Features

Table II summarizes our node features. Optype denotes the operation type of
each node, such as `add`, `mul`, `div`, `load`, `store`. This is obtained by
analyzing LLVM IR instructions. #invocation denotes the number of times that
an operation or a node is executed in a loop, which is associated with the
loop tripcount and unroll factor. In and out degrees represent the number of
edges flowing into and out of the node, respectively, which can indirectly
reflect the potential resource usage of multiplexers. Delay and #cycle
represent the timing information of each operation/node, which helps latency
prediction. We build a latency and delay library for each operation type based
on profiling results from several micro-benchmarks. For LUT, DSP and FF, we
profile corresponding resource usage for arithmetic operations (e.g., `add`,
`mul`, `fadd`) using micro-benchmarks and build a library. For non-arithmetic
operations (e.g., `br`, `icmp`, `mux`), we simply set their resource-related
features to zero.

Table II: Node features Feature | Description | Values
---|---|---
Optype | operation type |  `load`, `etc.`
#invocation | the number of invocations | `int`
In degree | the number of incoming edge | `int`
Out degree | the number of outcoming edge | `int`
#cycle | the number of clock cycles | `int`
Delay | the delay (ns) of each operation | `float`
LUT | LUT usage of each operation | `int`
DSP | DSP usage of each operation | `int`
FF | FF usage of each operation | `int`
Figure 3: The illustration of our hierarchical source-to-post-route QoR
modeling approach.

####  III-B2 Loop-level features

Besides node features, we extract loop-level features so as to consider the
effect of HLS pragmas on loops. As discussed in Section III-A, the graph of a
pipelined loop remains the same as that of the loop without pipelining.
Therefore, additional features should be considered to differentiate
pipelining and non-pipelining. We use three features that determine the
latency of a pipelined loop, namely, iteration latency (IL), initiation
interval (II) and tripcount (TC). Specifically, IL can be approximated by GNNs
and TC can be obtained by analyzing the LLVM IR code, while II can be computed
by the equation[10]:

$\displaystyle II_{\textit{min}}$ $\displaystyle=$
$\displaystyle\max(II^{rec},II^{res})$ $\displaystyle=$
$\displaystyle\max(\max\\{\lceil\frac{\textit{Delay}_{p}}{\textit{Distance}_{p}}\rceil\\},\max\\{\lceil\frac{\textit{Access}_{m}}{\textit{Ports}_{m}}\rceil\\})$

where Delayp represents the latency between a pair of dependent instructions
from different iterations, Distancep refers to the difference between the
corresponding iteration numbers, Portsm is the number of ports and Accessm is
the number of accesses to array $m$.

###  III-C The Hierarchical Modeling Approach

After constructing and annotating graphs, the complete C-to-bitstream flow is
invoked to obtain the ground truth QoR considering various HLS pragmas. Then
the annotated graphs and their corresponding QoR labels form our dataset for
model training. The detailed setup is discussed in Section IV-A.

We propose a hierarchical source-to-post-route QoR modeling methodology, as
depicted in Fig. 3. To capture and incorporate the inherent hierarchical
nature of loops, we train multiple GNN models to estimate QoR results at
different loop hierarchies and gradually predict the overall latency and
resource usage of the application. Our method can be generalized to
applications with multiple nested loops. Take the code in Fig. 3 as an
example, the top function foo has two nested loops, each with two loop levels,
i.e., loop hierarchies. In the first nested loop, its sub-loop at the inner
hierarchy (j-level) is pipelined and partially unrolled, and its sub-loop at
the outer hierarchy (i-level) is partially unrolled with factor two. This
configuration creates a circuit with two replicas of the inner loop pipeline
within the outer loop. The QoR of inner replicas greatly impacts the overall
QoR of the nested loop. Thus, our strategy models the QoR of sub-loops at the
inner hierarchy using a local GNN, merges inner sub-loops to super nodes, and
integrates super nodes with operations at the outer hierarchy to model overall
QoR with a global GNN.

This hierarchical approach simplifies post-route QoR prediction by dividing
the problem into smaller ones. It also improves estimation accuracy by
capturing fine-grained features at different hierarchies, compared to previous
GNN methods that directly predict performance with the whole graph. To be more
specific, we divide loop hierarchies into two levels: inner and outer
hierarchies.

####  III-C1 Inner Hierarchy

In our method, a loop at the inner hierarchy denotes the loop that only
contains computing logic without sub-loops. There are four types of loops that
belong to the inner hierarchy: ① a single-level loop; ② a nested loop with
loop pipelining applied to the outermost level (in this case, all the inner
sub-loops will be fully unrolled); ③ a perfect nested loop with loop
flattening and loop pipelining applied to the innermost level (in this case,
all the loop levels will be flattened to form a deeper pipeline which is
equivalent to a single-level pipelined loop); and ④ a nested loop with all the
inner sub-loops fully unrolled. For example, in Fig. 3, the two j-level loops
in red and green boxes are loops at the inner hierarchy, corresponding to
category ② and ①, respectively. At the inner hierarchy, related loops are
first extracted and subgraphs are constructed for model training, as
illustrated in Fig. 3. Then we train two different GNN models, namely
$\textbf{GNN}_{p}$ and $\textbf{GNN}_{np}$, to predict QoR for pipelined and
non-pipeline loops, respectively. This is because execution models of
pipelined and non-pipelined loops are different and training GNN models
separately can improve accuracy based on our experiments. After predicting the
QoR of loops at the inner hierarchy, the nodes in a subgraph are merged into a
super node, which represents this loop at the inner hierarchy when modeling
the outer hierarchy.

####  III-C2 Outer Hierarchy

We define all the remaining logic and loop levels as components at the outer
hierarchy. Note that loops can contain more than two loop levels. After the
first stage, we condense each loop at the inner hierarchy into a super node.
These super nodes are connected to the nodes from the outer hierarchy to form
the complete graph representation of the top function foo. If the outer loop
is unrolled, corresponding logic and super nodes are replicated following
strategies discussed in Section III-A. In the entire program graph, we
annotate super nodes with predicted QoR as features. Each super node holds a
complete set of node features as shown in Table II, where features such as
`LUT`, `DSP`, `FF`, `latency` are derived from the previous prediction phase.
The calculation method used for other features remains consistent with those
of standard logic nodes. This program graph, containing super nodes and logic
nodes of the outer hierarchy, becomes the input graph of our global GNN model
named $\textbf{GNN}_{g}$. The global GNN model is trained to predict the final
QoR of the whole application. This approach not only assists us in predicting
the overall performance of real applications more accurately but also helps in
discerning serial and parallel relationships between super nodes, i.e., loops.

Figure 4: GNN architectures and the corresponding hierarchical training
process.

###  III-D GNN Architectures and Training Process

We develop three GNN models for the prediction of each QoR metric, namely
$\textbf{GNN}_{p}$, $\textbf{GNN}_{np}$, and $\textbf{GNN}_{g}$. Due to
different execution models on hardware, $\textbf{GNN}_{p}$ and
$\textbf{GNN}_{np}$ are deployed for pipelined and non-pipelined loops at the
inner hierarchy, respectively. $\textbf{GNN}_{g}$ is used for prediction at
the outer hierarchy. Figure 4 illustrates the model architectures and our
hierarchical training process. The training process of the three GNN models is
carried out in a hierarchical way. We initially collect the loops at the inner
hierarchy, categorize them into pipelined loops and non-pipelined loops, and
extract information from the two kinds of loops to generate corresponding
datasets. Then $\textbf{GNN}_{p}$ and $\textbf{GNN}_{np}$ are trained
accordingly. Once the training process is completed, the weights of both
$\textbf{GNN}_{p}$ and $\textbf{GNN}_{np}$ are frozen. The output of these two
GNNs is then used as features of super nodes, which are combined with features
of other operations at the outer hierarchy to generate the dataset for
training $\textbf{GNN}_{g}$.

Our GNN models have similar architectures, as depicted in Fig. 4. Each model
can be decomposed into 1) feature encoder, 2) propagation layer, 3) pooling
layer, and 4) MLP layer. The feature encoder converts the attributes of nodes
into initial feature vectors. We use one-hot encoding for the attribute
`optype`, while numerical attributes are directly incorporated into the
feature vector. Propagation layers implement the message-passing mechanism of
GNNs. We adopt the propagation layers from typical GNN models: GCN [11], GAT
[12], GraphSAGE [13], TransformerConv [14], and PNA [15]. After passing
through three propagation layers, node embeddings are obtained. Then we use a
pooling layer to generate graph-level embeddings from node embeddings. We use
two types of pooling: sum pooling $h_{sum}$ and max pooling $h_{max}$. The
embeddings from these two pooling mechanisms are concatenated to form the
graph-level representation. Finally, MLP layers are used to project the graph-
level embeddings into QoR prediction. For resource usage, we directly obtain
the estimates after feeding the embeddings of the GNN into an MLP layer. For
latency prediction, we handle the inner hierarchy ($\textbf{GNN}_{p}$ and
$\textbf{GNN}_{np}$) and the outer hierarchy ($\textbf{GNN}_{g}$) separately.
Regarding the inner hierarchy, first, $\textbf{GNN}_{p}$ and
$\textbf{GNN}_{np}$ utilize an MLP layer to predict iteration latency. Second,
another MLP layer takes the predicted iteration latency and loop-level
features (elaborated in Section III-B) as input and infers the latency of the
loops at the inner hierarchy. Regarding the outer hierarchy, we directly
predict the overall latency with an MLP.

##  IV Experiments

###  IV-A Experimental setup

Our feature construction and model development flow is fully automated and
implemented using Python and C++. We utilize 16 applications from the
Polybench, MachSuite and CHStone benchmark suite. Specifically, 12
applications are used for GNN model training and testing, and 4 applications
are used for the DSE experiment. The datasets to develop $\textbf{GNN}_{p}$
and $\textbf{GNN}_{np}$ are constructed based on sub-loops extracted from the
application source code, while the dataset for $\textbf{GNN}_{g}$ is
constructed based on the complete applications. Our experiments are conducted
on the AMD Ultrascale+ MPSoC ZCU102. We utilize Vitis-HLS 2022.1 and Vivado
2022.1 to collect ground-truth latency and resource usage as training labels.
By applying various pragma combinations, we produce 3102, 2300 and 6178 valid
designs to build $\textbf{GNN}_{p}$, $\textbf{GNN}_{np}$, and
$\textbf{GNN}_{g}$, respectively. 80% of the dataset is used for training, 10%
of the dataset is used for validation, and the remaining 10% of the dataset is
used for testing. We use the mean absolute percentage error (MAPE) to quantify
the prediction quality. Each model is trained with over 250 epochs.

Table III: MAPE of post-route QoR with different GNNs. GNN type |  | Latency |  | Iteration
---
Latency
DSP | LUT | FF
GCN | $\textbf{GNN}_{p}$ | 5.23% | 7.21% | 4.94% | 6.63% | 7.83%
$\textbf{GNN}_{np}$ | 8.41% | 8.41% | 9.35% | 5.88% | 8.39%
$\textbf{GNN}_{g}$ | 9.85% | N/A | 7.45% | 11.27% | 11.24%
GAT | $\textbf{GNN}_{p}$ | 4.91% | 7.89% | 4.67% | 9.12% | 9.20%
$\textbf{GNN}_{np}$ | 9.42% | 9.42% | 10.64% | 6.99% | 9.32%
$\textbf{GNN}_{g}$ | 10.88% | N/A | 8.18% | 12.34% | 11.57%
GraphSage | $\textbf{GNN}_{p}$ | 5.24% | 5.52% | 5.88% | 8.12% | 8.66%
$\textbf{GNN}_{np}$ | 7.47% | 7.47% | 9.83% | 6.25% | 7.22%
$\textbf{GNN}_{g}$ | 8.55% | N/A | 6.94% | 9.86% | 9.99%
Transformer | $\textbf{GNN}_{p}$ | 4.92% | 5.95% | 4.05% | 6.70% | 7.12%
$\textbf{GNN}_{np}$ | 8.77% | 8.77% | 12.02% | 5.83% | 8.65%
$\textbf{GNN}_{g}$ | 8.54% | N/A | 6.97% | 11.33% | 10.31%
PNA | $\textbf{GNN}_{p}$ | 6.74% | 7.43% | 4.11% | 8.92% | 9.93%
$\textbf{GNN}_{np}$ | 9.99% | 9.99% | 10.55% | 7.89% | 10.26%
$\textbf{GNN}_{g}$ | 9.29% | N/A | 7.84% | 9.65% | 10.55%

###  IV-B Evaluation of QoR prediction accuracy

As introduced in Section III-D, we adopt five different GNN models and compare
their QoR prediction accuracy, as shown in Table III. For loop-level models
($\textbf{GNN}_{p}$ and $\textbf{GNN}_{np}$), we present prediction errors for
both latency and iteration latency. For application-level models
($\textbf{GNN}_{g}$), the latency prediction error of complete designs is
reported. Among all the GNN models for application-level prediction,
Transformer, GraphSage, and PNA achieve the lowest estimation errors in
latency (8.54%), DSP usage (6.94%) and FF usage (9.99%), and LUT usage
(9.65%), respectively. Our hierarchical method leads to high prediction
accuracy of post-route QoR for pipelined loops ($\textbf{GNN}_{p}$), non-
pipelined loops ($\textbf{GNN}_{np}$) and the complete applications
($\textbf{GNN}_{g}$), showing prediction errors of less than 10% for latency
and all types of resources. Since our method directly takes the C/C++ source
code as input to predict the post-route QoR of FPGA designs without invoking
any EDA tool flow, the prediction accuracy is highly desirable.

###  IV-C Comparison with the state-of-the-art

We compare our method with the approach [8] which also predicts post-route QoR
using GNNs. Different from[8], we additionally consider various pragma
combinations, increasing the complexity of code structures. For a fair
comparison, we evaluate both methods on the dataset from [8] which does not
consider HLS pragmas and our dataset with pragmas, as shown in Table IV.
Without pragmas, our method and [8] achieve close prediction accuracy, while
on the dataset with pragmas applied, our method demonstrates superior
improvement in all QoR metrics. We attribute this notable improvement to the
following two reasons: 1) joint representation of code and pragmas during
graph construction: we deploy a graph construction method that represents the
CDFG of source code while taking into account the impact of pragmas. This
approach effectively embeds the pragma information into the graph
representation. 2) reservation of loop hierarchies: with our hierarchical
training approach, messages can propagate more effectively between different
levels of loop hierarchies, allowing GNNs to better understand the cascading
effects of changes in the design’s inner hierarchy on the outer hierarchy.

Table IV: Comparison of prediction error (MAPE). Method | Configuration | Latency | DSP | LUT | FF
---|---|---|---|---|---
[8] | w/o pragma | N/A | 8.26% | 5.10% | 7.58%
Ours | w/o pragma | 5.54% | 6.71% | 6.78% | 7.91%
[8] | w/ pragma | 35.81% | 57.31% | 27.14% | 29.03%
Ours | w/ pragma | 8.54% | 6.94% | 9.65% | 9.99%

###  IV-D Evaluation of Design Space Exploration

To demonstrate the effectiveness of our hierarchical QoR prediction method, we
perform DSE on four new applications (`bicg`, `symm`, `mvt` and `syrk`) that
are unseen by our model during training. The design space of each application
is constructed by varying pragma configurations. Specifically, we apply HLS
pragmas like loop pipelining, loop flattening and loop unrolling iteratively
from inner to outer loops, with unroll factors from {1, 2, 4, 8, 16}. As for
array partitioning, we set partitioning factors consistent with unroll
factors.

The objective of DSE is to find the optimal pragma configurations producing
the best QoR metrics, which form the Pareto frontier. To quantify the
optimality, the average distance from reference set (ADRS) can be computed:

$\text{ADRS}(\Gamma,\Omega)=\frac{1}{|\Gamma|}\sum_{\gamma\in\Gamma}\min_{\omega\in\Omega}f(\gamma,\omega),$

where $\Gamma$ is the exact Pareto-optimal set, $\Omega$ is the approximate
Pareto-optimal set, and $f(\cdot)$ computes the distance between two design
points. A lower ADRS means that the approximate Pareto set is closer to the
exact one.

To get the exact Pareto-optimal set, we exhaustively evaluate the actual post-
route QoRs for each design configuration by invoking the complete C-to-
bitstream flow. Vitis HLS and Vivado are utilized to generate the accurate
results. We then select the models with the highest prediction accuracy in
Table III to provide quick QoR estimation for DSE. For comparison, two other
GNN models from [8] and GNN-DSE [6] are tested similarly on the same design
space. GNN-DSE [6] is one of the state-of-the-art works in HLS DSE. It uses a
GNN-based post-HLS QoR prediction model to address the DSE problem. The
comparison of DSE results is shown in Table V. The design space contains up to
2796 design configurations. The DSE time shows the timing costs required to
finish the exploration with the QoR feedback from Vivado and from our model,
respectively. By utilizing our method, efficient QoR feedback is provided,
only necessitating constructing an IR graph from the C-level code followed by
several rounds of GNN inference. This effectively bypasses the tedious steps
of synthesis and implementation and shortens the DSE time from tens of days to
tens of minutes.

Compared to the work [8] and GNN-DSE [6], which achieve average ADRS values of
13.94% and 10.96% respectively, our method reaches an average ADRS of 6.91%,
demonstrating our superior ability to enable a high-quality DSE. The result
fully verifies that the high accuracy of our predictive model can lead to a
more precise approximation to the exact Pareto frontier. Regarding the
runtime, both GNN-DSE and our method complete the exploration within minutes,
while the method[8] requires one to two days to complete due to the invocation
of HLS tools. In summary, by incorporating our hierarchical prediction method,
fast and accurate performance feedback can be provided to boost the efficiency
of DSE.

Table V: DSE results on unseen applications. Kernels |  #Design Config | DSE time | ADRS(%)
---|---|---|---
Vivado | Ours | [8] | [6] | Ours
bicg | 2796 | 26 days | 12 min | 12.13 | 9.14 | 6.39
symm | 1972 | 44 days | 15 min | 13.15 | 11.50 | 5.45
mvt | 2116 | 25 days | 14 min | 16.01 | 12.17 | 9.31
syrk | 1972 | 40 days | 15 min | 14.45 | 11.03 | 6.49

##  V Conclusion

In this paper, we introduce a GNN-based hierarchical method for FPGA post-
route QoR modeling from C-level source code. This efficient source-to-post-
route QoR prediction strategy notably reduces the design turn-around time of
FPGA HLS. Specifically, we utilize loop hierarchies during model establishment
to achieve better prediction accuracy and propose an effective graph
construction method to jointly represent code and HLS pragmas. Experiments
show that our method outperforms existing works in the prediction of latency
and resource usage, which can better support the DSE task.

## Acknowledgement

This work is partially sponsored by the National Natural Science Foundation of
China (62102249, 62232015) and Shanghai Pujiang Program (21PJ1408200). Mingzhe
Gao and Jieru Zhao contributed equally to this paper.

## References

  * [1] J. Zhao _et al._ , “Comba: A comprehensive model-based analysis framework for high level synthesis of real applications,” in _ICCAD_ , 2017.
  * [2] Y.-k. Choi _et al._ , “Hlscope+,: Fast and accurate performance estimation for fpga hls,” in _ICCAD_ , 2017.
  * [3] G. Zhong _et al._ , “Design space exploration of fpga-based accelerators with multi-level parallelism,” in _DATE_ , 2017, pp. 1141–1146.
  * [4] S. Dai _et al._ , “Fast and accurate estimation of quality of results in high-level synthesis with machine learning,” in _FCCM_ , 2018.
  * [5] H. M. Makrani _et al._ , “Pyramid: Machine learning framework to estimate the optimal timing and resource usage of a high-level synthesis design,” in _FPL_. IEEE, 2019.
  * [6] A. Sohrabizadeh _et al._ , “Automated accelerator optimization aided by graph neural networks,” in _DAC_ , 2022.
  * [7] L. Ferretti _et al._ , “Graph neural networks for high-level synthesis design space exploration,” _ACM TODAES_ , vol. 28, no. 2, 2022.
  * [8] N. Wu _et al._ , “High-level synthesis performance prediction using gnns: Benchmarking, modeling, and advancing,” in _DAC_ , 2022.
  * [9] Cummins _et al._ , “ProGraML: A Graph-based Program Representation for Data Flow Analysis and Compiler Optimizations,” in _ICML_ , 2021.
  * [10] J. Zhao _et al._ , “Performance modeling and directives optimization for high-level synthesis on fpga,” _IEEE TCAD_ , vol. 39, no. 7, 2020.
  * [11] T. N. Kipf and M. Welling, “Semi-supervised classification with graph convolutional networks,” in _ICLR_ , 2017.
  * [12] P. Veličković _et al._ , “Graph attention networks,” in _ICLR_ , 2018.
  * [13] W. L. Hamilton, R. Ying, and J. Leskovec, “Inductive representation learning on large graphs,” in _NIPS_ , 2017.
  * [14] Y. Shi _et al._ , “Masked label prediction: Unified message passing model for semi-supervised classification,” _arXiv preprint_ , 2020.
  * [15] G. Corso _et al._ , “Principal neighbourhood aggregation for graph nets,” in _NeurIPS_ , 2020.


# Low-temperature high-frequency dynamic magnetic susceptibility of classical
spin-ice Dy2Ti2O7

S. Teknowijoyo1,2, K. Cho1, E. I. Timmons1,2, M. A. Tanatar1,2, J. W. Krizan3,
R. J. Cava3 and R. Prozorov1,2 1 Ames Laboratory, Ames, IA 50011, USA
2Department of Physics & Astronomy, Iowa State University, Ames, IA 50011, USA
3 Department of Chemistry, Princeton University, Princeton, New Jersey 08544,
USA<EMAIL_ADDRESS>
###### Abstract

Radio-frequency (14.6 MHz) AC magnetic susceptibility, $\chi^{\prime}_{AC}$,
of Dy2Ti2O7 was measured using self-oscillating tunnel-diode resonator.
Measurements were made with the excitation AC field parallel to the
superimposed DC magnetic field up 5 T in a wide temperature range from 50 mK
to 100 K. At 14.6 MHz a known broad peak of $\chi^{\prime}_{AC}(T)$ from kHz -
range audio-frequency measurements around 15 K for both [111] and [110]
directions shifts to 45 K, continuing the Arrhenius activated behavior with
the same activation energy barrier of $E_{a}\approx 230$ K. Magnetic field
dependence of $\chi^{\prime}_{AC}$ along [111] reproduces previously reported
low-temperature two-in-two-out to three-in-one-out spin configuration
transition at about 1 T, and an intermediate phase between 1 and 1.5 T. The
boundaries of the intermediate phase show reasonable overlap with the
literature data and connect at a critical endpoint of the first order
transition line, suggesting that these low-temperature features are frequency
independent. An unusual upturn of magnetic susceptibility at $T\to 0$ was
observed in magnetic fields between 1.5 T and 2 T for both magnetic field
directions, before fully polarized configuration sets in above 2 T.

††: J. Phys.: Condens. Matter

##  1 Introduction

Spin-ice compounds with a magnetic structure resembling the crystal structure
of water ice have attracted notable interest since their identification in the
1990s [1, 2, 3]. A rich variety of materials with the generic formula
A2B2(O,F)7 (here A is a trivalent cation and B is a tetravalent cation) with
pyrochlore lattice provides a fruitful playground for the search of exotic
magnetic phases, including gapless spin liquid, quantum spin-ice, and magnetic
monopoles, just to name a few [4]. In the classical spin-ice compounds,
Dy2Ti2O7 (DTO) and Ho2Ti2O7 (HTO), the magnetic moments of the (Dy, Ho) rare
earth elements are localized at the vertices of the corner-shared tetrahedra
of the pyrochlore lattice [2, 3, 5]. Because of a strong single-ion
anisotropy, the ground state of Dy3+ is well expressed by an Ising doublet
with local [111] quantization axes [6]. An effective ferromagnetic nearest-
neighbor interaction leads to a geometric frustration, which prevents the
system from reaching a long-range ordered state [7]. While the magnetic
moments are macroscopically disordered, their local configurations obey the
“ice rule”: two spins point outward and two spins point inward in each
tetrahedron, which is analogous to the Bernal-Fowler rule for water ice [8].
As a result, the lowest energy state is degenerate with a residual entropy
very close to that of water ice [3]. The resultant magnetically frustrated
configuration is extremely sensitive to the competition between short-range
and long-range magnetic interactions, leading to a complex magnetic behavior
[1, 2, 3, 9, 10, 11].

The extensive spin degeneracy is lifted when external magnetic field is
applied. For a field along [100] direction, all four spins in an elementary
tetrahedron align in a way such that they have a component of magnetic moment
along the field direction, without breaking the ice rule [12]. For a magnetic
field along [111] direction, however, the response is more complex. For low
fields, the moments first align following the ice rule so that one of the four
spins on each tetrahedron is parallel to the field and forms a layer of
triangular lattice which alternates with a layer of kagome spin lattice
configuration [12, 13]. When the applied magnetic field exceeds roughly 1 T,
the phase with two-in-two-out (2-2) configuration of the spins undergoes a
first order transition to a phase with three-in-one-out (3-1) configuration at
low enough temperatures ($T\leq$ 0.36 K) [12].

Magnetization measurements with precisely aligned magnetic field direction
found that a small deviation of a magnetic field from [111] direction towards
[112] direction splits the 1 T transition into two [14]. The resulting
intermediate state was identified in specific heat [13] and susceptibility
measurements [11, 14, 15]. To explain this new state, a theoretical model
taking into account the effects of dynamic lattice distortions has been
proposed [11].

In regards to the spin relaxation time, different spin configurations and
transitions have very different characteristic time scale(s), from seemingly
static to very slow activated behavior, to fast driven transitions. As a
result, early studies of Dy2Ti2O7 reported a variety of frequency dependencies
of magnetic susceptibility, which showed at least two distinct peaks. The
first peak at around 1 K is nearly frequency-independent at audio-frequency
[16] and its signature is also seen in zero-field heat capacity measurements
[3, 10]. At $T<1$ K, most of the relaxation channels “freeze” out and the
relaxation time scale grows exponentially on cooling. The possible dynamics of
magnetic monopoles in this part of the phase diagram is an interesting
suggestion, which has recently started to be explored experimentally [17, 18].
The second peak in AC susceptibility measurements exhibits strong temperature
dependence, from $\sim$0 K (DC) to 20 K (10 kHz AC) [16, 19, 7, 20]. In the
MHz range, a µSR experiment found a gap in muon polarization between $40-60$ K
[21] for which the authors attributed that it may or may not come from the
same relaxation channel as found in the susceptibility studies. Measurements
at even higher frequencies using nuclear forward scattering of synchrotron
radiation [22] and nuclear magnetic resonance [23] also identify response even
at higher temperatures.

Motivated by these studies and to bridge the gap between the audio frequency
AC susceptibility and the higher frequency experiments, here we report radio-
frequency (14.6 MHz) magnetic susceptibility study of the classic spin-ice
compound Dy2Ti2O7 using highly sensitive tunnel diode resonator (TDR)
technique. Our zero field data show that the higher-temperature peak shifts to
45 K which is in a good agreement with the µSR experiment [21]. Furthermore,
we show that our peak temperature follows the same Arrhenius-type activated
behavior as in low-frequency measurements. This observation supports the
notion that the same energy barrier is responsible for this activated behavior
across the wide frequency range. At low temperatures, we carefully explore the
$H-T$ phase diagram down to $T_{\mathrm{base}}\approx 50$ mK and up to 2 T in
both [110] and [111] directions. The frequency-insensitive 1 K peak grows
stronger and we observe an intermediate phase near the spin flip transition
for $H\parallel$[111]. We find that the magnetic action continues beyond the
boundary of the intermediate phase before reaching fully polarized state above
2 T.

##  2 Experimental

Single crystals of Dy2Ti2O7 were synthesized using optical floating zone
method, as described in detail elsewhere [24]. The crystals were oriented
using x-ray Laue diffractometry. The samples for the TDR measurements were cut
and polished into the shape of cubes with sides of 0.7 mm long oriented along
[110] or [111] directions. The accuracy of this cutting is estimated to be
within 2 degrees from the determined direction [25].

The tunnel diode resonator technique is widely used for precision radio
frequency (rf) magnetic susceptibility measurements of magnetic and
superconducting materials with $1:10^{7}$ frequency resolution measuring
better than 0.1 Hz frequency shift at 10 MHz [34, 35, 30, 31]. While in
metallic samples the signal is a combination of skin effect and spin
susceptibility, in insulating materials, such as the present case of spin-ice,
it gives directly the real part of rf AC susceptibility up to a constant
geometric factor. The TDR spectroscopy has been successfully used for studies
of field-dependent phenomena and quantum magnetism, for example in magnetic
molecules (e.g.Cr12Cu2) [26, 27], diluted magnetic systems [28, 29] and FM/AFM
materials [30, 31, 32, 33]. Please see reviews [34, 35, 36] for the detailed
description of the principles of the TDR technique used in the present work,
and Ref. [37] for the technical discussion of the implementation of the TDR in
a dilution refrigerator, which was used in our experiments in this work.

In brief, a tunnel-diode, biased to its regime of a negative differential
resistance, drives an $LC-$ tank circuit consisting of a capacitor, $C$, and
an inductor, $L$, with the fundamental frequency $2\pi f_{0}=(CL_{0})^{-1/2}$
around 10-20 MHz. In our particular dilution refrigerator setup,
$f_{0}\approx$ 14.6 MHz. Essentially, the tunnel diode compensates for the
losses and, when the impedances properly match, triggers spontaneous self-
oscillations at $f_{0}$. It is possible to maintain this state with minimal
external noise and excellent stability, provided that the circuit itself is
well stabilized and shielded [37, 38]. The outstanding sensitivity of the
technique comes from the fact that it operates in the frequency domain
(similar to microwave cavity perturbation technique and NMR), but without the
need to scan the frequency (or field) in search for the resonance, which is an
added benefit. A properly biased TDR system self-resonates and is always
locked on its own resonance. If a sample is inserted into the coil, the total
(coil + sample) inductance changes from $L_{0}$ to $L=L_{0}+\Delta L$. In all
practical setups, the inductance shift is small, $\Delta L\ll L_{0}$. The new
resonant frequency is $2\pi f=(C(L_{0}+\Delta L))^{-1/2}\approx f_{0}(1-\Delta
L/2L_{0}+3/2(\Delta L/2L_{0})^{2}+\mathcal{O}(L/2L_{0})^{3})$. Experimentally,
the difference $\Delta f=f-f_{0}$ is in the order of 1 kHz. Since $f_{0}\sim$
10 MHz, $\Delta f/f_{0}\approx 10^{-4}$ and therefore the quadratic and higher
order corrections can be safely omitted.

The exact sign and value of $\Delta f$ are dependent on the susceptibility of
the sample. A positive $\Delta f$ (i.e., negative $\Delta L$) implies a
diamagnetic response and the opposite occurs for a paramagnetic one. As a
function of temperature, the relation between the frequency change to the
susceptibility can be written as $\Delta f(T)=-G4\pi\chi^{\prime}_{AC}(T)$.
The constant of proportionality, $G=f_{0}V_{s}/2V_{c}(1-N)$, is determined by
the demagnetization factor $N$, the sample volume $V_{s}$ and the coil volume
$V_{c}$. This constant can be calibrated experimentally by measuring the
absolute value of the DC susceptibility of the same sample in a Quantum Design
Magnetic Properties Measurement System (MPMS). It can also be determined
directly by pulling the sample out of the coil. The important requirement for
the high resolution TDR measurements is the temperature stability of the coil
and the circuit within a mK range. This is achieved by thermally separating
the sample stage and the circuit parts as described in detail in [37]. The
noise level of our TDR setup is $\sim$0.05 Hz out of 14.6 MHz, which means
that the sensitivity of $\sim$3 ppb is achieved.

To cover an extended temperature range, we use two separate TDR setups in our
study, one mounted in a 3He refrigerator fitted with a 9 T superconducting
magnet (and equipped with a mechanical manipulator to perform sample pull-out)
and another in a dilution refrigerator working with a 14 T magnet [37]. In
both setups, the static field $H_{DC}$ is always parallel to the rf field
$H_{AC}$ ($\sim$20 mOe). Measurements in the dilution refrigerator are
performed at temperatures as low as 50 mK and up to 3 K. We therefore achieve
a good temperature overlap with the 3He system used from 400 mK to 120 K. For
Dy2Ti2O7, the low temperature region is of a particular interest because of
the slowing down and freezing of the spin dynamics [5, 19, 20] may allow the
system to reach a magnetically ordered ground state in an unusual way.

DC magnetic susceptibility measurements were performed in Quantum Design MPMS
SQUID magnetometer. These measurements were performed on the same samples as
used in TDR study, and were used for calibration purposes for radio frequency
measurements.

##  3 Results

The solid lines in Figure 1 show the temperature-dependent magnetic
susceptibility of Dy2Ti2O7 measured at 14.6 MHz in magnetic fields up to 9 T
applied along the [111] direction, using the 3He TDR setup. For reference, we
show the $\chi^{\prime}_{AC}(T)$ data by Snyder et al (dashed lines) measured
at 100 Hz [19]. The data for TDR spectroscopy are calibrated using the SQUID
DC magnetic susceptibility measurements (gray solid line). The previously
reported “15 K feature” in low-frequency ($<10$ kHz) susceptibility
measurements [16, 19] shifts to $\sim 45$ K at 14.6 MHz and becomes less
dependent on magnetic field. This temperature is also consistent with a
previously reported $40-60$ K gap in the muon polarization in a similar
frequency range [21]. As shown by the dashed line in the inset of Figure 1,
the position of the TDR peak (red star with yellow interior) extends the
frequency dependence found in the low-frequency measurements [16, 19] almost
in an exact linear fashion. The Arrhenius-type activation energy calculated
using the formula $f=f_{0}$ exp$(E_{a}/k_{B}T)$ yields $E_{a}/k_{B}=$
(slope)/log${}_{10}(e)\approx 230$ K and $f_{0}\sim 2$ GHz. The value of the
activation barrier is consistent with the theoretically calculated gap between
the ground state and the first excited doublets in the spin-ice system [39]
and therefore, the peak in $\chi^{\prime}_{AC}(T)$ has been attributed to the
Orbach relaxation process between these states [21]. These mechanisms are also
thought to be responsible for the results coming from experiments at even
higher frequencies such as nuclear forward scattering [22] and nuclear
quadruple resonance [23], which trace a higher energy barrier with their
combined average is shown as dotted line in the inset ($E_{a}/k_{B}\sim 365$
K). This difference could be related to the higher accessible excited crystal
energy levels (due to different experimental parameters such as higher
frequency and temperature), and the higher sensitivity to the local nuclear
spin and quadruple interactions around the Dy3+ ions of the of the NFS and NQR
techniques.

Figure 1: (Color online) Temperature-dependent radio frequency (14.6 MHz)
magnetic susceptibility $\chi^{\prime}_{AC}(T)$ of Dy2Ti2O7 (solid lines)
measured in magnetic fields $H\parallel$[111]. The gray solid line shows the
$\chi^{\prime}_{DC}(T)$ of the same sample measured in a DC SQUID magnetometer
for calibration. For reference, we show the 100 Hz $\chi^{\prime}_{AC}(T)$
data from Ref. [19] (dashed lines). The inset shows an Arrhenius plot of the
position of the high temperature peak on the measurement frequency. Peak
position in our RF measurements (red star with yellow interior) continues the
trend determined in audio frequency measurements (open squares [19]) with the
same activation energy ($E_{a}$ $\approx$ 230 K, dashed line). The data from
nuclear forward scattering (open magenta triangles) [22] and nuclear quadruple
resonance (open blue circles) [23] experiments deviate slightly to a higher
value of activation energy (dotted line).

In Figure 2, we compare TDR data for two experimental configurations,
$H\parallel$[111] (solid lines, left axis) and $H\parallel$[110] (dashes,
right axis) over a broad temperature range $0.4-120$ K. The curves for
different external DC magnetic fields are shifted vertically, from top to
bottom, $|H|$ = 0, 1, 2, 3, 4 and 5 T. The normalization factor between the
two orientations is determined by matching the zero field curve for [110] to
its [111] counterpart between 10 K and 50 K (around the Arrhenius-type spin
relaxation peak, because its position and shape are nearly independent of
field orientation and strength), and the same normalization factor is then
applied to the rest of the [110] curves. Above 50 K, the susceptibilities for
the two orientations differ slightly presumably because of a small difference
in the physical geometry of the two samples. A previous study on sample shape
dependence of the magnetization measurements reported that spheres are highly
desirable for the calculation of the exact (analytical) demagnetization
correction [40]. For our cubic samples, the simple demagnetization correction
$H_{i}\approx H-4\pi MN$ (where $H_{i}$ is the effective magnetic field sensed
by the sample and $N$ is the geometric demagnetization factor) is no longer
exact albeit still useful. This correction is applied for the region of our
interest near the $\chi^{\prime}_{AC}(T)$ peak at $T\leq 1$ K in
$H\parallel$[111] $\approx 1$ T (see Figure 6(a) and the Discussion section
below), which corresponds to a spin-flip transition from 2-in-2-out to
3-in-1-out configuration [12, 41].

Figure 2: (Color online) Temperature-dependent radio frequency (14.6 MHz)
magnetic susceptibility $\chi^{\prime}_{AC}(T)$ of Dy2Ti2O7 measured in
magnetic fields $H\parallel$[111] (solid lines, left axis) and
$H\parallel$[110] (dashed lines, right axis) directions. The curves for
different applied magnetic fields are shifted vertically for clarity. For the
sake of comparison of two independent measurements performed on crystals of
different orientation, the right axis is adjusted so that $H\parallel$[110]=0
curve matches its [111] counterpart at the 45 K peak feature. The other
dominant feature of $\chi^{\prime}_{AC}(T)$ is a sharp maximum observed at
$\sim$1 K only in $H\parallel$[111] configuration in fields close to 1 T. The
inset shows measurements on warming between 2.6 K to 50 K in $H\parallel$[111]
= 1.2 T magnetic field for both ZFC (black) and FC (red) conditions. No
obvious thermal hysteresis is found for both low-temperature and high
temperature features.

Our data show no obvious thermal history dependence, as shown in the inset of
Figure 2. Here we compare the data for measurements in $H\parallel$[111]
$=1.2$ T taken for both zero-field cooled (the magnetic field is switched on
after initial cool down at 2.6 K) and the consecutive in-field cooled run. The
shown ZFC and FC curves are taken on warming with the same heating protocol in
order to make their comparison straightforward.

Figure 3: (Color online) Radio frequency magnetic susceptibility of Dy2Ti2O7
measured in dilution refrigerator and 3He setups in 50 mK to 8 K temperature
range. Measurements in configuration $H\parallel$[111] are presented by solid
lines with respect to the left axis, measurements in $H\parallel$[110] are
shown with dashes with respect to the right axis. Panel (a) shows data for
fields $H\leq$ 0.8 T and (b) for fields $H\geq$ 0.9 T in the vicinity of the
spin-flip transition. Both left and right scales in (b) are $\sim$4 times
larger to accommodate the height increase of the peaks.

Figure 3 shows in more detail the combined $\chi^{\prime}_{AC}(T)$
measurements from the dilution refrigerator setup (50 mK $-$ 3 K) and the 3He
setup ($0.4-8$ K) in an external DC magnetic fields between $0-2$ T. It has a
similar construction to the Figure 2, where the solid lines, left axis (dashed
lines, right axis) correspond to the [111] ([110]) magnetic field direction
and the curves are offset vertically for clarity. The same normalization
factor used in Figure 2 is applied again to the [110] curves for consistency
in the comparison between [111] and [110] directions across different figures
and panels. Panel (a) shows data for fields below 0.8 T, where the difference
between the two experimental configurations is not very pronounced. The zero-
field curves for both configurations show a broad peak at around 2 K, followed
by an increase of $\chi^{\prime}$ response for $T\leq 0.5$ K. Both features
are of a broad cross-over type, consistent with no observed long-range
magnetic ordering down to 50 mK. Under increasing fields, the peak at 2 K
broadens and both the peak and the onset of the upturn shift slightly to
higher temperatures. The broad peak feature is similar to previous low-
frequency susceptibility measurements [7, 16, 19, 20] which is related to 2-2
to 3-1 spin flip transition. As $T\rightarrow 0$, the increase of
susceptibility in almost all the curves in the figure bears similarity to the
heat capacity data [10] whose nature has been suggested to be related to the
onset of phase transition to its (long-range ordered) ground state [10],
Schottky-type peak due to multitude of local orders [11], the dynamics of spin
monopoles [17, 18], or the presence of other spin relaxation channels [42,
43].

Figure 4: (Color online) (a) Temperature-dependent real part of the radio
frequency AC magnetic susceptibility of Dy2Ti2O7 taken in magnetic fields
along the [111] direction in the $1.2-2$ T range. Non-monotonic features in
the temperature dependencies are observed at $T\sim$ 0.4 K (circles),
suggestive of an intermediate phase between $1.2-1.5$ T range which is absent
in the 2 T curve. (b) An enlarged view of the susceptibility data in 1.3 T,
emphasizing the “kink” feature that is also clearly seen as a drop in the
temperature-dependent susceptibility derivative (right axis).

Figure 3(b) shows measurements for $0.9~{}\mathrm{T}\leq|H|\leq 2$ T. In the
vicinity of the spin-flip transition, $\chi^{\prime}_{AC}(T)$ in
$H\parallel$[111] configuration shows a significant increase (note the curve
at 1.1 T) and the build up of a strong peak at 0.7 K. The feature at 1.1 T is
consistent with a step in the DC magnetization $M(H)$ measurements, as
reported in the literature [12, 41]. For fields of 1.2 T and above, the
susceptibility gradually decreases, however, this decrease reveals some new
sharp “knee” features in the $\chi^{\prime}_{AC}(T)$ and its derivative
(Figure 4(a) and (b)), suggestive of a phase transition. None of these
features are observed for the $H_{\mathrm{DC}}\parallel$[110] $\geq 1$ T
configuration, for which the susceptibility shows a smooth evolution and
remains small. A specific heat Monte Carlo study has predicted a phase
transition at low temperature when an external field is applied in the [110]
direction [44] which is notably absent in our data and to the best of our
knowledge, still eludes an experimental verification [15].

Figure 5: (Color online) Field-dependent radio frequency magnetic
susceptibility as determined from $T$-sweep data of Figure 3 at characteristic
temperatures of 100 mK (black solid circles), 300 mK (yellow-red squares), 400
mK (blue solid down-triangles) and 500 mK (green open up-triangles). The curve
for 200 mK data is presented in Figure 6(a). The locations of the peaks in
susceptibility are marked as blue up-triangles in the phase diagram Figure
6(b). At 500 mK the curve forms a single peak with stronger relative intensity
as shown in the inset.

To get an additional insight into the nature of the kink feature observed
above 1.1 T, we plot in Figure 5 the temperature evolution of the field-
dependent $\chi^{\prime}_{AC}(H)$ as extracted from $T$-sweep data of Figure
4. The dependence shows a singular peak at 400 and 500 mK which loses
intensity and splits into two peaks at 1.2 T and 1.5 T below 0.3 K. This
evolution is similar to the two-peak features found in the low frequency AC
susceptibility data [11, 15], see below.

##  4 Discussion

Working with magnetic samples, even paramagnetic ones, there is always a
concern regarding the demagnetization correction. Most papers, however, do not
discuss this and it could be hard to follow their data to the actual magnetic
field sensed by the spins. Therefore any comparison between different sets of
data should be taken with this uncertainty in mind. In our case, we used cube-
shaped samples. Indeed, while true demagnetization correction is only valid in
ellipsoidal samples where magnetic induction is constant inside the sample, it
can still be used in the effective, integral sense to relate a total magnetic
moment of a given sample with the effective change of integral magnetic
susceptibility, $dM/dH$, where $H$ is the applied magnetic field. The shape-
dependent correction has been recently calculated for arbitrary-shaped
samples, including most frequent in the research cuboidal shapes [45]. Of
course, the next best after sphere is a cube, which at least has isotropic
corrections along three principal axes. This is the case in the present work.
The effective demagnetization factor for a cube along any of the primary
directions is $N=0.39$ [45]. Note that the sum of the demagnetization factors,
$\sum_{i}=0.39\times 3=1.17>1$. This sum is equal 1 only for ellipsoids and
this caused significant confusion in the literature trying to deal with non-
ellipsoidal shapes. Using cgs notation, the corrected magnetic field, “seen”
by the sample is $H_{i}=H-4\pi MN$, where $M=m/V$ is volume magnetization and
$m$ is the measured magnetic moment in emu=erg/G. In our case, the sample has
dimensions of $0.74\times 0.70\times 0.70$ mm3, so that its volume is
$V=3.63\times 10^{-4}$ cm3. Taking DC magnetization measured in Quantum Design
MPMS on this sample, the $M(H)$ curve saturates at $m_{\mathrm{sat}}=$ $0.22$
emu, roughly above $20$ kOe of the applied field at 5 K. In the region of
interest, at 12 kOe, it is already $m=$ $0.188$ emu which is not far below
$m_{\mathrm{sat}}$. This implies that we can use a simple approximation of
$\mu H/k_{B}T\gg 1$ for higher fields and lower temperatures. Importantly,
being close to saturation means that the demagnetization correction amounts to
a simple practically temperature and field - independent constant shift and,
therefore cannot contribute to any new features in our temperature - dependent
measurements. Using $m=0.188$ emu, the constant demagnetization correction is
roughly, $4\pi mN/V\approx 2.5$ kG. For smaller magnetic fields and higher
temperatures, this correction becomes smaller, because magnetization
decreases. Of course, for small fields, demagnetization correction will stay
temperature dependent down to the lowest temperatures, but it does not affect
higher-field constant correction.

Figure 6: (Color online) (a) Demagnetization correction of 0.25 T applied to
the 200 mK data. The corrected peak locations are indicated by arrows. (b)
Phase diagram of Dy2Ti2O7 as determined from heat capacity ([13], circles) and
low-frequency AC magnetic susceptibility ([15], squares) measurements in
$H\parallel$ [111] configuration. These points form the spin-flip transition
boundary (purple dashed line) and the kagomé-intermediate phase boundary
(green dashed line). Red-yellow down-triangles denote the kinks in the
$T$-sweeps from Figure 4a and and blue up-triangles denote the demagnetization
corrected positions of the peaks of $H$-sweeps from Figure 5.

Figure 6 shows the demagnetization correction applied to the 200 mK data and
the summary of the corrected special points in the $\chi^{\prime}_{AC}(T)$ and
$\chi^{\prime}_{AC}(H)$ data. The red down-triangles represent the points of
anomalies in the $T$-sweeps (Figure 4), while the blue up-triangles are from
the peaks in $H$-sweeps (Figure 5). Note that the splitting of the critical
point into two lines separated by approximately 0.25 T below 0.4 K is a
characteristic feature of the intermediate phase as found in heat capacity
(circles, from [13]) and low-frequency AC susceptibility data (squares, from
[15]). Our TDR data demonstrate reasonable overlap with these measurements
which suggests that, unlike the 15 K feature, the peak splitting at low
temperatures is largely frequency independent (up to the rf range). This
splitting feature is observed in magnetic fields slightly inclined from [111]
toward [112] direction and the magnitude of the splitting depends on the angle
of misalignment [11, 14]. The splitting in our data, $\sim$0.25 T, corresponds
to inclinations between 5 and 7 angular degrees [14] which is a reasonable
estimate of the combined accuracy of our TDR coil-sample alignment and the
sample crystallographic axes determination using Laue diffraction patterns.

Both our misalignment-induced splitting and the slight upturn at low
temperatures are consistent with the distortion-based dipolar spin-ice model
(d-DSM) [11], which builds upon previous theoretical models [46, 47, 48, 49].
This slight increase of magnetic susceptibility on cooling is observed in
magnetic fields of both [111] and [110] orientations, above 1.5 T boundary of
the intermediate phase (Fig. 3). It is generally believed that a spin
polarized state is achieved above the intermediate phase, in which case rf
susceptibility should be zero [31]. Non-zero value suggest that the magnetic
polarization is not complete, at least up to 2 T.

##  5 Summary

To summarize, we have extended the dynamic magnetic AC susceptibility
measurements into the radio-frequency range by using the tunnel diode
resonator technique to investigate a classical spin-ice compound Dy2Ti2O7. We
observe a peak in susceptibility at 45 K, adding a point that follows the
linear behavior in the Arrhenius plot ($\log_{10}f$ vs. $1/T_{\mathrm{peak}}$)
from the published low frequency data, thus indicating that a single
activation barrier around 230 K for the Orbach processes related to the
crystal field energy. This value is slightly lower than the barrier reported
by yet higher frequency nuclear forward scattering [22] and nuclear quadruple
resonance [23] experiments. The difference may be due to the difference in the
experimental parameters, contributions of higher order terms and the local
sensitivity of these probes. We also find that the AC susceptibility increases
as $T$ goes below 0.5 K for all fields including above the spin-flip
transition. This is inconsistent with the conventional spin-polarized state
which could be an onset of a phase transition or unusual spin relaxation
channels such as the proposed spin monopole dynamics. “Kink-like” and peak-
splitting features are visible in the susceptibility for a magnetic field
along the [111] crystallographic direction in a 1.2 to 1.5 T range. After
correcting for demagnetization effects, these points reveal a reasonable
overlap with the boundaries found in previous specific heat and magnetic
susceptibility studies on the $H-T$ phase diagram, suggesting that these
features arise from frequency-insensitive mechanisms. Our results are
consistent with the d-DSM model, which has been proposed as the underlying
microscopical theory of the classical spin ice Dy2Ti2O7.

## Acknowledgment

We thank R. Flint and M. J. P. Gingras for illuminating discussions. This
research was supported by the U.S. Department of Energy, Office of Basic
Energy Sciences, Materials Science and Engineering Division through the Ames
Laboratory [S.T., K.C., E.I.T., M.A.T., R.P.]. The Ames Laboratory is operated
for the U.S. Department of Energy by Iowa State University under Contract No.
DE-AC02-07CH11358. Work at Princeton University [J.W.K., R.J.C.] was supported
by DOE BES grant number DE-SC0019331.

## References

## References

  * [1] Gaulin B D, Reimers J N, Mason T E, Greedan J E and Tun Z 1992 Phys. Rev. Lett. 69(22) 3244–3247
  * [2] Harris M J, Bramwell S T, McMorrow D F, Zeiske T and Godfrey K W 1997 Phys. Rev. Lett. 79 2554–2557
  * [3] Ramirez A P, Hayashi A, Cava R J, Siddharthan R and Shastry B S 1999 Nature 399 333–335
  * [4] Gingras M J P and McClarty P A 2014 Rep. Prog. Phys. 77 056501
  * [5] Gingras M J P 2009 Science 326 375–376
  * [6] Rau J G, Petit S and Gingras M J P 2016 Phys. Rev. B 93(18) 184408
  * [7] Fukazawa H, Melko R G, Higashinaka R, Maeno Y and Gingras M J P 2002 Phys. Rev. B 65(5) 054410
  * [8] Bernal J D and Fowler R H 1933 J. Chem. Phys. 1 515–548
  * [9] Fennell T, Petrenko O A, Fåk B, Gardner J S, Bramwell S T and Ouladdiaf B 2005 Phys. Rev. B 72(22) 224411
  * [10] Pomaranski D, Yaraskavitch L R, Meng S, Ross K A, Noad H M L, Dabkowska H A, Gaulin B D and Kycia J B 2013 Nat. Phys. 9 353–356
  * [11] Borzi R A, Gómez Albarracín F A, Rosales H D, Rossini G L, Steppke A, Prabhakaran D, Mackenzie A P, Cabra D C and Grigera S A 2016 Nat. Commun. 7 12592
  * [12] Sakakibara T, Tayama T, Hiroi Z, Matsuhira K and Takagi S 2003 Phys. Rev. Lett. 90(20) 207205
  * [13] Higashinaka R, Fukazawa H, Deguchi K and Maeno Y 2004 J. Phys. Soc. Jpn. 73 2845–2850
  * [14] Sato H, Matsuhira K, Sakakibara T, Tayama T, Hiroi Z and Takagi S 2007 J. Phys.: Condens. Matter 19 145272
  * [15] Grigera S A, Borzi R, Slobinsky D G, Gibbs A S, Higashinaka R, Maeno Y and Grigera T S 2015 Papers in Physics 7 070009–070009
  * [16] Matsuhira K, Hinatsu Y and Sakakibara T 2001 J. Phys.: Condens. Matter 13 L737
  * [17] Toews W H, Reid J A, Nadas R B, Rahemtulla A, Kycia S, Munsie T J S, Dabkowska H A, Gaulin B D and Hill R W 2018 Phys. Rev. B 98(13) 134446
  * [18] Stöter T, Doerr M, Granovsky S, Rotter M, Goennenwein S T B, Zherlitsyn S, Petrenko O A, Balakrishnan G, Zhou H D and Wosnitza J 2020 Phys. Rev. B 101(22) 224416
  * [19] Snyder J, Slusky J S, Cava R J and Schiffer P 2001 Nature 413 48–51
  * [20] Snyder J 2004 Phys. Rev. B 69 064414
  * [21] Lago J, Blundell S J and Baines C 2007 J. Phys.: Condens. Matter 19 326210
  * [22] Sutter J P, Tsutsui S, Higashinaka R, Maeno Y, Leupold O and Baron A Q R 2007 Phys. Rev. B 75(14) 140402
  * [23] Kitagawa K, Higashinaka R, Ishida K, Maeno Y and Takigawa M 2008 Phys. Rev. B 77(21) 214403
  * [24] Krizan J W and Cava R J 2014 Phys. Rev. B 89 1–5
  * [25] Zeisner J, Bräuninger S A, Opherden L, Sarkar R, Gorbunov D I, Krizan J W, Herrmannsdörfer T, Cava R J, Wosnitza J, Büchner B, Klauss H H and Kataev V 2019 Phys. Rev. B 99(15) 155104
  * [26] Engelhardt L, Martin C, Prozorov R, Luban M, Timco G A and Winpenny R E P 2009 Phys. Rev. B 79(1) 014404
  * [27] Martin C, Engelhardt L, Baker M L, Timco G A, Tuna F, Winpenny R E P, Tregenna-Piggott P L W, Luban M and Prozorov R 2009 Phys. Rev. B 80(10) 100407
  * [28] Kim H, Bud’ko S L, Tanatar M A, Avdashchenko D V, Matovnikov A V, Mitroshenkov N V, Novikov V V and Prozorov R 2012 J. Supercond. Nov. Magn. 25 2371–2375
  * [29] Novikov V, Avdashchenko D, Bud’ko S, Mitroshenkov N, Matovnikov A, Kim H, Tanatar M and Prozorov R 2013 Philos. Mag. 93 1110–1123
  * [30] Vannette M D and Prozorov R 2008 J. Phys.: Condens. Matter 20 475208
  * [31] Vannette M D, Bud’ko S L, Canfield P C and Prozorov R 2008 J. Appl. Phys. 103 07D302
  * [32] Vannette M D, Sefat A S, Jia S, Law S A, Lapertot G, Bud’ko S L, Canfield P C, Schmalian J and Prozorov R 2008 Journal of Magnetism and Magnetic Materials 320 354–363
  * [33] Quirinale D G, Messina D, Rustan G E, Kreyssig A, Prozorov R and Goldman A I 2017 Phys. Rev. Appl. 8(5) 054046
  * [34] Prozorov R, Giannetta R W, Carrington A and Araujo-Moreira F M 2000 Phys. Rev. B 62(1) 115–118
  * [35] Prozorov R and Giannetta R W 2006 Supercond. Sci. Technol. 19 R41
  * [36] Prozorov R and Kogan V G 2011 Rep. Prog. Phys. 74 124505
  * [37] Kim H, Tanatar M A and Prozorov R 2018 Rev. Sci. Instr. 89 094704
  * [38] Van Degrift C T 1975 Rev. Sci. Instr. 46 599–607
  * [39] Rosenkranz S, Ramirez A P, Hayashi A, Cava R J, Siddharthan R and Shastry B S 2000 J. Appl. Phys. 87 5914–5916
  * [40] Bovo L, Jaubert L D, Holdsworth P C and Bramwell S T 2013 J. Phys.: Condens. Matter 25 3–9
  * [41] Aoki H, Sakakibara T, Matsuhira K and Hiroi Z 2004 J. Phys. Soc. Jpn. 73 2851–2856
  * [42] Shannon N, Sikora O, Pollmann F, Penc K and Fulde P 2012 Phys. Rev. Lett. 108(6) 067204
  * [43] Melko R G, den Hertog B C and Gingras M J P 2001 Phys. Rev. Lett. 87 067203–
  * [44] Ruff J P C, Melko R G and Gingras M J P 2005 Phys. Rev. Lett. 95(9) 097202
  * [45] Prozorov R and Kogan V G 2018 Phys. Rev. Appl. 10(1) 014030
  * [46] Siddharthan R 1999 Phys. Rev. Lett. 83 1854–1857
  * [47] den Hertog B C and Gingras M J P 2000 Phys. Rev. Lett. 84 3430–3433
  * [48] Bramwell S T and Gingras M J P 2001 Science 294 1495–1501
  * [49] Yavors’kii T, Fennell T, Gingras M J P and Bramwell S T 2008 Phys. Rev. Lett. 101 037204