\section{Credit card fraud detection entry point: dataset and evaluation}
\label{sec:intro}

\subsection{Credit card transactions dataset}

Payment data contains very sensitive private information about individuals and businesses and access to such data is highly restricted to only the data owners and the companies that manage the data. Which is the reason why there is no publicly available data set that would sufficiently reflect the magnitude and variety of real-world card payments and that could therefore be considered as a basis for studying the many interesting challenges present in this domain.

The only credit card transactions we were able to find\footnote{Besides the confidential one used during our joint research project with Worldline} was a small sample of 280.000 transactions with PCA-transformed features\footnote{\url{https://www.kaggle.com/mlg-ulb/creditcardfraud}, Last access: 10.05.2019.} published alongside a study which was conducted jointly together with the exact same payment service provider~\citep{dal2015calibrating}. This scarcity of publicly available real-world data is problematic for the research community, since it compromises transparency, comparability and reproducibility of published findings.

\begin{table} \centering
    \begin{tabular}{llp{0.35\linewidth}p{0.2\linewidth}}
    \toprule
    Category & Attribute & Description & Feature(s)\\
    \midrule
    Transaction & \underline{Time} & Date and time of the transaction  \vspace{0.3cm} & \texttt{weekday}, \texttt{hourOfDay}, \texttt{timestamp}\\
               		  & \underline{Type} & The kind of transaction issued. A real payment (e.g. e-commerce) or a virtual transaction (3D-secure ACS).  \vspace{0.3cm} & \texttt{ecom}, \texttt{3DSecureACS} \\
               		  & \underline{Authentication} & The kind of authentication experienced by the card holder: Chip-based PIN check, signature, 3D-secure. \vspace{0.3cm} & \texttt{emv}, \texttt{authentication}, \texttt{3DSecure} \\
               		  & Entry mode & The interaction modality between the card and the terminal: e.g. contactless, magentic stripe, manual entry.  \vspace{0.3cm} & \texttt{cardEntryMode} \\
               		  & Amount & The transacted monetary value.  \vspace{0.3cm} & \texttt{amount} \\
               		  & Fraud label & The manually assigned label indicating whether a transaction was a fraud. & \texttt{fraud}\\
    Card holder & \underline{Declarative profile} & Several attributes describing the card holder such as his country of residence, age or gender. \vspace{0.3cm} & \texttt{userCountry}, \texttt{age}, \texttt{gender}\\
    					  & \underline{Card} & Several attributes describing the credit card such as the credit limit, card type or expiration date. \vspace{0.3cm} & \texttt{creditLimit}, \texttt{cardType}, \texttt{expiryDate}\\
    Merchant & \underline{Declarative profile} & Two attributes describing the terminal: The merchant category code assigned to the terminal (e.g. shoe store, ATM, etc.) and the country the terminal is registered in (e.g. Belgium, Germany, etc.) \vspace{0.3cm} & \texttt{terminalCategory}, \texttt{terminalCountry}\\
     \bottomrule
    \end{tabular}
    \caption[Attributes of transactions]{Transactions are described along three groups of attributes: Transaction related, card holder related and merchant related. Underlined attributes are in fact collections of several attributes. Examples of attributes in the collection are listed in the description.}
    \label{tab:attributes}
\end{table}

As listed in~\cref{tab:attributes}, transactions are characterized by several attributes. We can broadly categorize them into attributes characterizing the transaction itself, attributes characterizing the card holder/card that issues a transaction and attributes characterizing the merchant who receives the transaction. In addition, each card holder and each merchant is assigned a unique identifier.

The fraud label is binary and it indicates the authenticity of a transaction. The label was assigned by human investigators after consultation with the card holder. The label takes the value \texttt{0} in most of the cases. Only when a transaction has been identified as having been issued by a third party and not the rightful card holder, is the label set to \texttt{1}, which flags the transaction as fraudulent. Throughout the paper we will refer to the label as the \textit{fraud label} or the \textit{class} of the transaction.

\subsection{Performance Measures} \label{sec:performance}

\subsubsection{Confusion matrix based evaluation}
\label{sec:confmat}
Traditional ways of evaluating machine learning classifiers use the confusion matrix describing the difference between the ground truth of the dataset and the model prediction.

\begin{figure}[h]
\centering
\includegraphics[width=1.0\linewidth]{confusion-matrix.png}
\caption[Confusion matrix]{A confusion matrix displaying the four possible outcomes of a binary decision (colored boxes) and an excerpt of measures that quantify different properties of the decision maker in terms of ratios of these outcomes. Type 1 errors (false positives) and type 2 errors (false negatives) are shaded in gray.}
\label{fig:contigency-table}
\end{figure}

Using the confusion matrix, one can compute a variety of metrics in order to compare classifiers performances:
\begin{description}
\item[Precision ] Proportion of well classified examples among examples of the positive class: $$prec=\frac{TP}{TP+FP}$$
\item[Recall ] Also called true positive rate, proportion of examples of the positive class well-classified among all positive examples: $$rec=\frac{TP}{TP+FN}$$.
\item[Accuracy] Proportion of well classified examples among all the testing examples (similar to precision but for every class): $$acc=\frac{TP+TN}{TP+TN+FP+FN}$$
\item[F1 score] Precision and recall have opposite behavior: forcing a good precision is detrimental for the recall and vice versa. F1 score is the harmonic mean between precision and recall and therefore gives equal importance to precision and recall. For this reason, this metric which takes values between 0 and 1 is often considered a good by default metric for unbalanced classification tasks: $$F1=2*\frac{Precision * Recall}{Precision + Recall} = 2* \frac{2TP}{2TP+FP+FN}$$
\item[MCC] Matthews correlation coefficient is another score that is robust to imbalance \cite{david2011}. It is related to the chi-square statistic and takes values between -1 and 1. For binary classification: $$MCC=\frac{TP \times TN - FP \times FN}{\sqrt{(TP+FP)(TP+FN)(TN+FP)(TN+FN)}}$$
\end{description}

Confusion matrix based metrics are a good starting point for evaluating classifiers' performances. However these metrics are point estimators and some users may want to have a fine grained knowledge of how the classifiers are performing over all the domains of the testing set. Parametric evaluation of classifier performance allows for that more fine grained look at the classifiers' behaviours.

\subsubsection{Parametric evaluation}
\label{sec:parametric}
Some machine learning classifiers like random forest, neural networks, SVM, logistic regression classifier, etc. output a set of $n$ predictions for each example corresponding to their confidence that each example belongs to each of the $n$ classes.
\newline

By changing the value of a threshold parameter through all the possible values of probabilities obtainable with the machine learning classifier, one can monitor the efficiency of the prediction for each examples: from the most obvious to classify to the least obvious. In fine, parametric evaluation can be seen as an evaluation of all the possible confusion matrices obtainable given a prediction vector.
\newline

When doing a parametric evaluation, two curves are usually plotted: the ROC (Receiving Operating Characteristic) curve and the precision-recall curve
\begin{itemize}
\item The ROC curve's axes are: true positive rate (i.e. recall) for the y-axis and false positive rate for the x-axis. Both axes are proportion of the class they are representing: positive for the true positive rate and negative for the false positive rate. Therefore, the ROC curve is insensitive to unbalanced datasets: differences in class effectives won't affect the overall shape of the curve.
\item The Precision-Recall curve axis are the precision for the y-axis and the recall for the x-axis. As for the axis of the ROC curve, the recall is a ratio of the positive class: it represents the proportion of example of the positive class identified. However, the precision is based on both positive class (TP) and negative class (FP) effectives. Therefore, a strong class imbalance (such as credit card fraud detection imbalance) would affect the precision values since the classification is harder. The precision recall curve is more punitive than ROC curve and consequently more challenging to improve: in a context of class imbalance, the precision value is closer to 0 and therefore the precision recall area under the curve decreases. Besides, we experienced that the precision recall curve is easier to read for industrial partners.
\end{itemize}

\begin{figure}[h]
\centering
\includegraphics[width=0.8\linewidth]{metrics.png}
\caption{Precision recall curve and area under the curve \textit{Top: Precision recall curve and precision at $K$, indicated by grey lines. Bottom: Area under the precision recall curve and area in the early retrieval range, calculated only over the interval $[0, 0.2]$.}}
\label{fig:aucpr}
\end{figure}

According to domain application needs, the area under a portion of the curve can also be calculated. For example, for credit card fraud detection evaluation, the need to calculate the AUC for values of recall below 0.2 was discussed since the experts can only investigate a limited number of alerts. Some others

Further details about the intricacies in interpreting ROC curves versus precision recall curves have been discussed in~ \citep{NoROCSaito2015} and~\citep{Davis2006}. Although the interpolation between points requires special care, the area under the precision recall curve (AUCPR) can serve as a single-valued aggregate measure of a classifier's performance~\citep{Hand2009}.
In~\cref{fig:aucpr} we illustrate several measures that can be derived from a precision recall curve:
\begin{itemize}
\item Precision at $K$ (P@K): The proportion of true-positives in the set of $K$ transactions with the highest scores. Each choice of $K$ corresponds to a specific level of recall on the x-axis.
\item Area under the precision recall curve (AUCPR): The integral of precision over all levels of recall over the interval $[0, 1]$.
\item Area in the early retrieval range (AUCPR@0.2): The integral of precision over recall levels from $0$ to $0.2$. Instead of reporting precision at several values of $K$, we use this early retrieval range to reflect the application-specific focus on high precision in the range of low recall.
\end{itemize}

\section{Learning to classify unbalanced datasets} \label{sec:imbalance}
Typical fraud detection tasks present a strong imbalance in the representations of the classes: Usually there is a majority genuine class that outweighs the minority fraudulent class \cite{Ali2019}.
\newline

Most learning algorithms are not adapted to a strong imbalance in the representations of classes in the dataset. Therefore, learning from unbalanced dataset is a challenging research question that needs to be considered (\cite{batista2000}). \cite{weiss2001} and \cite{estabrooks2004} observed that an unbalanced dataset leads to a decrease in the performance of standard learning algorithms such as SVM or random forest. In fact, class imbalance have been shown to produce a negative effect on the performance of a wide variety of classifiers (\cite{chawla2005}). \cite{japkowicz2002} have shown that skewed class distributions affects negatively the performance of decision trees and neural networks. \cite{visa2005} drew the same conclusion for neural networks. \cite{kubat1997}, \cite{mani2003} and \cite{batista2004} showed the weakness of k Nearest Neighbors algorithms to imbalanced datasets. \cite{yan2003} and \cite{wu2003} showed that SVM were also weak to skewed class distributions. Overall, most of the classifiers suffer from class imbalance, some more than others.
\newline

The strategies used to tackle this issue operate at two different levels: some aim to solve the imbalance of the dataset beforehand whereas others want to adapt the learning algorithm to make it robust to an unbalanced classification task (\cite{chawla2004}, \cite{pozzolo2015}).
\newline

The data level strategies are called like that since they happen during the preprocessing time, before any learning algorithm is applied. They consist in rebalancing the dataset in order to compensate the structural imbalance. The algorithmic level strategies involve methods where the algorithms are designed especially to deal with an unbalanced task and methods where the classification costs of the algorithm are adapted in order to reflect the priorities of the unbalanced classification task. The latter are called cost sensitive classifiers (\cite{elkan2001}).
\newline

The sampling based solutions to handle imbalance will be described in section \ref{sec:sampling_desc}. The model based methods will be briefly highlighted in section \ref{sec:algo_desc}

\subsection{Sampling methods}
\label{sec:sampling_desc}
Sampling methods consist in reducing the class imbalance by changing the class effectives in the dataset. In this section we will first present two naive methods that reduce class imbalance in the dataset by adjusting class ratios: undersampling and oversampling. Afterwards we will introduce the Synthetic Minority Oversampling Technique strategy (SMOTE) that aims to create new examples of the minority class (\cite{chawla2002}). The different strategies are illustrated in figure \ref{fig:SMOTE}.
\newline

Undersampling consists in rebalancing the dataset by taking less examples from the majority class in order to match the absolute number of examples of the minority class (\cite{drummond2003}). The assumption made with undersampling is that the majority class contains redundant information that can be removed. The main issue with undersampling is that, when the imbalance is too pronounced, too many examples from the majority class need to be taken away. This may cause a decrease in the performance of the algorithm due to a lack of data.
It is worth mentioning that undersampling speeds up the learning phase, which makes it an interesting choice when in presence of an unbalanced dataset.
\newline

Oversampling consists in rebalancing the dataset by repeating random examples of the minority class (\cite{drummond2003}). The main issue of oversampling is that repeating some examples doesn't add information and may lead to an overfitting of the learning algorithms. Besides, oversampling increases the learning time since it makes the train-set artificially bigger.
\newline

In 2002, \cite{chawla2002} introduce SMOTE, a strategy to reduce class imbalance by creating new examples in the neighbourhood of the minority class. The creation of new minority class example is done by averaging neighboring examples of the minority class. For categorical features, the averaging consists in a majority voting for the category: the category mostly represented in the neighbourhood becomes the category of the corresponding categorical feature of the newly created example. They show an increase in the performance of classifier when using this rebalancing strategy. However SMOTE does not consider the label of neighboring examples when creating examples from the minority class. This can lead to an overlap between minority and majority classes. ADASYN \cite{he2008} and Borderline-SMOTE \cite{han2005} tackle this issue by taking into account only neighboring examples from the minority class.

\begin{figure}[t]
\centering
\includegraphics[width=\textwidth]{SMOTE.png}

\caption[Sampling strategy for unbalanced classification.]{Sampling strategy for unbalanced classification. \textit{(The $+$ symbols refer to the majority class whereas the $-$ symbols refer to the minority class. The symbols colored in red are newly created symbols by the oversampling strategy or by the SMOTE strategy(\cite{chawla2002}))}}
\label{fig:SMOTE}
\end{figure}

\subsection{Model based methods}
\label{sec:algo_desc}

Some solutions were published in order to strengthen ensemble based methods when in presence of class imbalance. EasyEnsemble (\cite{liu2009}), UnderBagging (\cite{wang2009}) and SMOTEboost (\cite{chawla2003}) are ensemble based method that are adapted by design to skewed class distributions. SMOTEboost (\cite{chawla2003}) combines boosting with SMOTE in order to provide balanced datasets to each boosting weak learner. Similarly, EasyEnsemble (\cite{liu2009}) uses undersampling in order to create different balanced datasets for the boosting learners. In UnderBagging (\cite{wang2009}) undersampled datasets are used to train several learners whose predictions are aggregated with bagging.
\newline

An other leverage for dealing with class imbalance is to adapt the metrics used to evaluate the models and the costs used to train some of them in order to counterbalance class imbalance. Indeed, some real world imbalanced classification tasks such as credit card fraud detection give more importance to correctly identify the element of the minority class than those of the majority class. In the case of cost sensitive learners, the misclassification cost can be adapted in order to reflect this priority. \cite{shirazi2010} present cost sensitive SVMs. \cite{kukar1998} introduced earlier the possibility of cost adaptation for neural networks. Adapting the misclassification cost for imbalanced problem has also been studied for boosting methods (\cite{sun2007}, \cite{wei1999}).

\section{Dataset shift detection and adaptation} \label{sec:drift}

Many machine learning tasks make the assumption that the random variables in the data are i.i.d. (independent and identically distributed) among the train and the test sets. That property ensures that the decision function learned by the classifier on the train set is valid on the testing set. However in some case the i.i.d. hypothesis isn't valid: there are changes in the data distribution between the training and testing set. These changes of the data distribution can decrease drastically the performances of the classifiers (\cite{pozzolo2014}, \cite{abdallah2016}, \cite{rodriguez2008}) and need to be taken into account.
\newline

The term "dataset shift" is used in order to include the phenomenons where the i.i.d. hypothesis isn't valid. It is defined by \cite{torres2012} as "cases where the joint distribution of inputs and outputs differs between training and test stage"
\newline

Dataset shift of various types have been shown to affect classifiers performance. \cite{rodriguez2008} showed that complex classifiers are more robust to population drift than simple classifiers.

\subsection{Different types of dataset shift}
Dataset shift is a broad word that refers to a change in the joint distribution of the target variable and the input features $p(x,y)$. However, since $p(x,y)\quad = \quad p(y|x)*p(x)$, the causes of the change in the joint distribution can be different. Several shifting phenomenon can be identified (\cite{torres2012}, \cite{storkey2009}, \cite{gao2007}, \cite{shimodaira2000}):
\begin{itemize}
\item \textit{Covariate shift} was formerly called \textit{population drift}. It refers to the phenomenon when the conditional distribution of the classes with respect to the training data doesn't change from the train set to the test set, but the distribution of the data do. For example in the credit card fraud detection settings, the phenomenon where the buying behaviors change over seasons but the fraudulent strategies don't can be seen as \textit{covariate shift}. It is formally defined as:
\begin{equation}
p_{train}(y|x) = p_{test}(y|x) \quad \textrm{and} \quad p_{train}(x) \ne p_{test}(x)
\end{equation}
\item \textit{Prior probability shift} is an other subtype of dataset shift. It is also sometimes called \textit{shifting priors}. \textit{Prior probability shift} refers to a change in the class distributions: for example when the proportion of fraudulent transactions in the training set is way smaller than for the testing set. It may cause a over-representation or an under-representation of one class in the prediction of machine learning classifiers. It is formally defined as:
\begin{equation}
p_{train}(y|x) = p_{test}(y|x) \quad \textrm{and} \quad p_{train}(y) \ne p_{test}(y)
\end{equation}
\item  \textit{Concept drift} is the type of shift that appears the most in the credit card fraud detection literature.\textit{Concept shift} or more often \textit{concept drift} refers to the phenomenon where the conditional distribution of the target classes with respect to the input variables changes. For example when fraudster adapt their strategies and the classifier is not able to detect fraudulent transactions efficiently anymore. It is formally defined as:
\begin{equation}
p_{train}(y|x) \ne p_{test}(y|x) \quad \textrm{and} \quad p_{train}(x) = p_{test}(x)
\end{equation}
\end{itemize}

\subsection{Detecting dataset shift}

Most of the studies passively detect dataset shift when a decrease in classifiers performance on the test set is witnessed(\cite{pozzolo2015}, \cite{gomes2017},...). In order to prevent the decrease of the classifier performance, some teams aimed to estimate more realistically the expected performance of classifiers in presence of dataset shift. For example, \cite{webb2005} pointed out the weaknesses of the ROC evaluation in presence of \textit{prior probability shift}. Earlier, \cite{kelly1999} stated that the weaknesses of classifiers in presence of covariate shift and prior probability shift can be highlighted with the decrease of classifiers performance for their loan risk prediction task. Additionally, \cite{sugiyama2007} discussed the issues appearing when classifier performance expectation is set using a cross validation hyper-parameter optimization: cross validation expectations can't be trusted if the distribution changes between train and test set. They proposed a cost based approach by comparing the data distribution $p(x)$ in the train and test set and weighting more the data points of the train set that are similar to those of the test set for the cross validation optimization. If the data point is typical of the train and test distribution, then it is weighted more for the evaluation of the model.
\newline

\cite{barddal2019} provided interesting results where they showed how to describe in advance the concept shift in two dataset (Airlines and Sea \footnote{\url{https://moa.cms.waikato.ac.nz/datasets/}}) by looking directly at the changes in the decision function through the changes of feature importance in an adaptive random forest classifier (described in next paragraph) \cite{gomes2017}. They could therefore predict the direction of the dataset shift and the resulting change in the conditional distribution $p(y|x)$.  Earlier, \cite{lane1998} tried to quantify the direction of drift on a continuous one dimensional axis by looking at which examples get misclassified when concept drift happens. An optimization of the threshold value of their classifier was managed in the decision process of their user identification for computer security task. Similarly, \cite{wang2003} wanted to precisely characterize concept shift by comparing the class prediction between two classifiers trained before and after a concept shift (old training data vs recent training data). The examples that present a conflict for class prediction can then be studied in order to understand which changes happen between the conditional distributions $p_{old}(y|x)$ and $p_{new}(y|x)$
\newline

\citep{lucas2019} introduced a technique to detect and quantify the extent of covariate shift between credit card transactions from different days. For each pair of days $(a, b)$:
\begin{itemize}
\item they collect all transactions from these days $D = D_{a} \cup D_{b}$ and label the transactions according to the day they came from.
\item They partition the sets of transactions and train a classifier to distinguish the day from which each transaction comes.
\item They evaluate the classifier performance and use this value as a distance measure between days: the better the performance, the easier the differentiation of the days, the more different the transactions of these days, and vice versa.
\end{itemize}
Through manual inspection and after clustering, the authors can show that low ”day distinguishability performance” clearly corresponds to calendar clusters of weekdays, Sundays and work holidays or Fridays and Saturdays (see figure \ref{fig:lucas2019}).

\begin{figure}[H]
\centering
\includegraphics[width=0.95\textwidth]{zoom_DS.png}
\caption[Zoomed distance matrix (covariate shift) for the face-to-face transactions.]{Weekly periodicity observed in the dataset shift for face-to-face transactions. \textit{(centile based coloration: green $\Leftrightarrow$ similar days, red $\Leftrightarrow$ dataset shift} \cite{Lucas2020}}
\label{fig:lucas2019}
\end{figure}

\subsection{Strategies in presence of dataset shift}

\paragraph{Data based methods}
Some authors chose to handle dataset shift by adapting the training and testing set in order to erase the shift between these sets.

\cite{gretton2009} used reproducing kernel Hilbert spaces (RKHS) in order to weight instances in the training set so that the averages of the training set features correspond to those of the test set features. This kernel mean matching approach is mostly used in the context of covariate shift since it allows to make the distribution of the data in the train and test set more similar. \cite{bickel2007} leveraged this kernel mean matching approach for data streams. They derived the kernel mean matching approach in order to integrate it in the logistic regression classifier used and to update it over time. Similarly to the kernel mean matching approach, \cite{torres2013} proposed a method in order to achieve transfer learning between different cancer diagnosis datasets. They used genetic programming in order to create a regressor that transforms the features of one dataset to the features of an other dataset in order to be able to apply a pretrained classifier without adapting it.
\newline

In addition to kernel mean matching techniques, \cite{klinkenberg2004} propose to train SVMs with an adaptive time window in order to prevent covariate shift.

\paragraph{Classifier based methods}
An other way to prevent concept drift is to update regularly the models in order to adapt them to the dataset shift.
\newline

At first, sliding windows based approaches have been explored. The idea is that by retraining them regularly using fresh data, the classifiers will stay up to date with dataset shift events. \cite{widmer1996} used a sliding window strategy for their FLORA4 framework: they trained regularly new classifiers on recent and trusted data but stored the bad performing ones for hypothetical later usage. \cite{lichtenwalter2009} also leveraged the sliding window strategy for imbalanced data streams with concept drift: they re-trained regularly C4,5 decision trees classifiers in order to prevent covariate shift. Moreover, they adapted the information gain calculation used when building decision trees by incorporating the Hellinger distance  in order to increase the robustness of the C4,5 trees to class imbalance and to prior probability shift. The Hellinger distance is a symmetric and non-negative measure of similarity between two probability distributions. The equation \ref{fig:hellingerdist} corresponds to the Hellinger distance between two discrete probability distributions with equal effectives $P=(p_{1}, ..., p_{k})$ and $Q=(q_{1},..., q_{k})$.

\begin{figure}[H]
\begin{equation}
d_{H}(P, Q) = \sqrt{\sum \limits_{i \in V} (\sqrt{P})-\sqrt{Q})^2}
\end{equation}
\caption{Hellinger distance: Measure of distributional divergence between two probability distributions P and Q}
\label{fig:hellingerdist}
\end{figure}

\cite{pozzolo2015} proposed a sliding window method for data streams with delayed information on the same credit card transactions dataset used in this thesis (see figure \ref{fig:window}). Since the ground truth for the class of credit card transactions is delayed by a week in average, the labels of the most recent transaction aren't known. They proposed to aggregate the predictions of a classifier on the old ground truth example with the predictions of a classifier trained on the feedbacks of the investigators on alerts raised previously by the fraud detection system. They showed that aggregating predictions provided better fraud detection than pooling all the examples of the ground truth and investigators feedback together.

\begin{figure}[H]
\centering
\includegraphics[scale=0.35]{pozzolo-sliding_window.png}
\caption{Aggregating ground truth examples with investigators feedbacks (\cite{pozzolo2015})}
\label{fig:window}
\end{figure}

Instead of preventing dataset shift by constantly re-training the classifiers, one may want to retrain them only if they prove to be lackluster. \cite{kolter2003} proposed to train again the classifiers when their performance drop. In order to prevent missclassification for data streams, \cite{hoens2011} used Hellinger distance on all the features of the dataset in order to detect early a covariate shift and which features are affected by it. They would therefore be able to know which trees of a random forest classifier will be affected by the dataset shift and need to be re-trained.
\newline

More recently, authors in \cite{gomes2017}, \cite{barddal2019bis} proposed an autonomous architecture for random forest classifier introducing a variant to Hoeffding Adaptive Tree in order to include a drift detector inside decision nodes to monitor the internal error rates of the tree, thus, continuously selecting features and building its predictive model over time. The Hoeffding Adaptive Trees (also called Very Fast Decision Trees) leverage the observation that using only a subset of examples in order to build a decision tree is sufficient to obtain good prediction \cite{bifet2009}. This computationally efficient and memory efficient tree architecture is therefore an asset in dataset shift contexts where models or parts of the models have to be retrained regularly.

\section{Feature Engineering} \label{sec:features}

\subsection{Feature Encoding}
Most learning algorithms can not handle categorical variables with a large number of levels directly. Both supervised~\citep{kotsiantis2007} and unsupervised methods~\citep{lam2015clustering} require preprocessing steps to transform them into numerical features. In this section, we discuss several common solutions for encoding categorical attributes such that they can be used in learning algorithms.

\paragraph{Label Encoding and One-hot encoding} The most direct encoding of categorical attributes with $K$ different values is to create discrete variables $x \in \{1, 2, \dots, K\}$. Depending on the learning algorithm, label encoding might not be adequate to represent categorical attributes. The distance on integers typically does not reflect the distance on the attribute's values and would cause misleading results in distance-based algorithms. On the other hand, One-Hot encoding is particularly convenient and commonly used in machine learning literature~\citep{bishop2006pattern}. A categorical variable is represented by a $K$-dimensional vector $\mathbf{x}$, in which one element $x_k$ is set to $1$ and all remaining elements are set to $0$.
\paragraph{Likelihood Encoding} (also known as target encoding) The value of a nominal attribute is encoded by its relative frequency of showing up in the positive class in a training set. If the $k$-th state of a nominal attribute appears $n_p$ times in the positive class and $n -n_p$ times in the negative class, the state gets assigned the value $x_k = \frac{n_p}{n - n_p}$. This type of encoding requires class labels and it assumes that the nominal values' representation in both classes provides sufficient information to distinguish the values. Therefore, we consider this technique rather as a feature engineering method than an actual encoding technique.
\paragraph{Vector Embeddings}
The values of categorical attributes may not always be orthogonal as suggested by a one-hot coding scheme but they rather exhibit some notion of similarity relevant to the prediction task. For instance, an attribute like the terminal's country can be endowed with different notions of similarity based on geographical locations, economic alliances or cultural backgrounds. It is common practice in machine learning to account for such expert knowledge through purposefully encoded attributes. Even if the main motivation is not an integration of a-priori knowledge, categorical attributes with many different levels partition the set of examples into many small sets for which it might be challenging to obtain reliable parameter estimates. Vector embeddings are the outcome of unsupervised encoding techniques with the aim to embed the values of a categorical variable in a continuous valued vector space. The role of unsupervised methods is thereby to estimate a map $c: X \to V^d$ from categorical values $X$ to points in a $d$-dimensional vector space $V^d$ such that distances between points reflect some sought notion of similarity between categorical values. We keep this superficial description of vector embeddings because there is an abundance of techniques by which such maps can be created, each with its own definition and interpretation of $c$ and $V$. For instance, the rows and columns in a tf-idf - matrix~\citep{JONES1973619} in information retrieval, the projections on the first $d$ principal components in dimensionality reduction, the mixture coefficients in probabilistic models with mixture distributions~\citep{blei2003latent} or the weights in neural networks~\citep{mikolov2011}. Similar to the likelihood encoding, the continuity in the embedding space is assumed to permit reasonable interpolation between otherwise distinct categorical values. \cite{russac2018} for example, successfully used Word2Vec embeddings in order to replace one-hot encoding for a distance based classifier (logistic regression) in a credit card fraud detection task.

\subsection{Feature Engineering}

In credit-card fraud detection, a fraud is supposedly not exclusively a property of the transaction itself but rather a property of both the transaction and the context it appeared in. While a transaction might arouse suspicion if, for instance, the spent monetary amount is very high and it takes place at a particular type of merchant at a certain time of day, such fraud detection strategy considers transactions in isolation from each other. Neither the historic purchases of the customer nor other transactions from the merchant he interacted with are taken into consideration. Two transactions occurring in geographically distant locations at almost the same time are often considered as suspicious patterns ("collisions"). Also, it may be that the volume of transactions per time (in an account) is more indicative of fraud than the absolute attributes of any transaction~\citep{whitrow2008}. In order to go beyond transaction-level detection, it is crucial to choose a suitable transaction context and extract relevant features from it.

In general, most studies seem to agree that the Recency - Frequency - Monetary Value (RFM) framework provides the necessary information for discovering frauds. Recency refers to the time since the last transaction, frequency captures the volume of transactions per time interval and the monetary value corresponds to the total money spent within the interval. Some studies consider only frequency and monetary value \citep{krivko2010hybrid}, whereas others also integrate the time difference (\citep{vlasselaer2015} \citep{bahnsen2016}). From this perspective, aggregations are convenient because they allow to condition the RFM statistics on transaction types. It should be noted that there is no consensus on how to choose the most suitable time intervals. Usually, researchers apply reasonable rules of thumb and compute aggregates over days, weeks or months.

\paragraph{Rules}

\paragraph{Time differences} Since frauds tend to appear in bursts of several transactions within a short time window, the time difference between consecutive transactions seems to be a promising indicator of fraud. For any two consecutive transactions $\mathbf{x}_{t-1}$ and $\mathbf{x}_t$ with their time stamps denoted by $\mathbf{x}_{t-1}^{(\operatorname{Time})}$ and $\mathbf{x}_{t}^{(\operatorname{Time})}$ the \texttt{tdelta} feature of the $t$-th transaction is:
\begin{equation} \label{eq:tdelta}
\operatorname{tdelta}_t = \mathbf{x}_{t}^{(\operatorname{Time})} - \mathbf{x}_{t-1}^{(\operatorname{Time})}
\end{equation}
These time differences feature have been shown to provide relevant insight to classifiers in a credit card fraud detection problem (\citep{jurgovsky2018}).

\paragraph{Feature aggregations} Most feature engineering work in credit card fraud detection follows transaction aggregation strategies such as the ones described by \citep{whitrow2008}. The idea is to characterize the spending pattern in an account, the so called \emph{activity record}, by accumulating transaction information over time. Through exploratory data analysis, the authors identified several salient aspects of transactions as being relevant to fraud detection, such as the number of recently issued transactions and the total monetary value of these transactions. It is important to note that aggregates do not require any label information and, therefore, are entirely unsupervised. While aggregates function as account-level features, they can be turned into transaction-level features by updating the activity record with each new transaction. Most commonly, the statistics from the activity record are simply added as supplementary features to the basic feature set when training a fraud classifier. This aggregation strategy was found to work well in commercial applications and it is also an integral part of research in the academic field with broad usage in a number of studies (\citep{bhattacharyya2011},\citep{jha2012},\citep{bahnsen2013},\citep{sahin2013cost},\citep{pozzolo2014}).

Since then, the proposed aggregation strategy served as generic template for developing new features. For instance, \citep{jha2012} and \citep{bahnsen2016} do not aggregate over specific transaction types but instead on those transactions whose attribute values are identical to the current transaction. This difference is subtle but it has severe implications for the semantics of the aggregated value. For example, \citep{whitrow2008} calculate the number of transactions issued in the last 24 hours at a specific merchant (e.g. "MyFancyShoeStore") for all transactions in the dataset. The value of this aggregate is probably zero for most transactions of most users but its semantic is consistent across all transactions. The aggregates defined in \citep{jha2012}\citep{bahnsen2016} calculate the number of transactions issued in the last 24 hours at the \emph{same} merchant as the current transaction. Finally, feature aggregates can be easily integrated in any existing classification process and they seem to readily improve the prediction accuracy.

\paragraph{Trading entropy}  \citep{Fu2016} propose \emph{trading entropy}, a feature that is supposed to quantify the amount of information added to a user's purchase history when he issues a new transaction. As an example, they define the feature based on the amount spent and the merchant. In a user's history, they consider all but the most recent transaction and calculate how the total amount is distributed over the merchants. They quantify the uncertainty in this distribution by means of the Shannon entropy. The trading entropy is then the difference between the entropy of this distribution and the distribution they obtain after adding the most recent transaction to his history.

\paragraph{Anomaly scores} Another interesting approach for extracting expressive features is discussed in \citep{vlasselaer2015}. Apart from the detection system the authors maintain a bipartite graph in which nodes correspond to merchants or card holders and edges correspond to transactions between these entities. The edge weight is determined by the volume of transactions exchanged between a card holder and a merchant, and it decays exponentially with time. With a page-rank style algorithm the authors propagate the fraud label of a transaction to directly connected or transitively reachable card holder and merchant nodes. The authors extract network features from the graph which measure the current exposure of each node to fraud. These features include a score for the card holder, the merchant and the transaction, aggregated over short, medium and long time periods.

\paragraph{External data integration} \citep{wetice2017} explored external data sources for augmenting transactions with information beyond the raw attributes. They extracted indicators of public holidays from calendars and used them as additional features of transactions. They also extracted embeddings from the DBpedia knowledge graph for encoding the originally categorical country tokens as points in a continuous vector space. Their results suggest a better generalization of the trained classifier and they hypothesize that this improvement results from the continuum between semantically related countries. Later, \citep{jurgovskythesis} supported this hypothesis by encoding country tokens in terms of several demographic statistics, such as a country's population, its gross domestic product or the corruption index. This explicit encoding led to a comparable performance increase as the unsupervised country embedding approach. Overall, encoding categorical tokens with semantic embeddings that can be extracted in an unsupervised fashion from public knowledge bases appears promising because such an approach does not require any manual effort and it readily improves the detection performance.

\section{Sequence modeling for fraud detection} \label{sec:sequence}

Sequence modeling is a major machine learning research question. In this section we will split the field with respect to three main approaches: the similarity-based approach that aims to learn the joint distribution of the events in a sequence in order to compute a similarity score, the model-based approach that aims to find the labels of the events in a sequence using data-driven strategies and the counting-based approach that aims to assess a likelihood score by counting the number of occurrences in a reference dataset.

\subsection{Similarity-based approaches}
\subsubsection{Graphical models}
Hidden Markov models (HMMs) are the most used model in the field of anomaly detection for discrete sequences according to \cite{chandola2012}.

They comprise two matrices (see figure \ref{fig:HMM1}):
\begin{itemize}
\item The \textit{transition matrix} describes the succession of the hidden states. It reflects a multinomial distribution from all hidden states at time $t$ to all hidden states at time $t+1$. The hidden states obey to the Markov property. Indeed, the hidden state at the time $t+1$ only depends on the hidden state at the time $t$: There is no memory in the distribution of successive hidden states, only conditional probability given the immediate previous hidden state. Formally, it means that conditional distributions between hidden states $t$ in a sequence of size $n$ can be simplified as following:
$$p(t_{k}|t_{k-1},t_{k-2},..,t_{1},w_{k-1},w_{k-2},w_{1}) = p(t_{k}|t_{k-1}) \;  k\in\{1,..,n\}$$
\item The \textit{emission matrix} describes the conditional distribution of the observed variable $w_i$ given the hidden states $t_i$. The probability distribution changes with respect to the observed variable properties. For example it can be considered multinomial for categorical observed variables or gaussian for continuous observed variables.
\end{itemize}

\cite{srivastava2008} used multinomial HMMs in the context of credit card fraud detection. They modeled the sequence of amount of users as a categorical feature with 3 values: "low amount" for the transactions with amounts less than 50 \euro{}, "mid amount" for transactions with amounts between 50 and 200 \euro{} and "big amount" for transactions with amounts bigger than 200 \euro{}. The anomaly score raised is the inverse of the likelihood of the recent amounts to have been generated by an HMM trained on a set of genuine sequences. They generated a credit card transactions dataset in order to assess the qualities of their strategies. \cite{dhok2012} used a similar approach for a real world credit card transaction dataset, and \cite{robinson2018} extended that with success for a prepaid credit card fraud detection task.
\newline

\cite{dorj2013} also used HMMs in an anomaly detection task in electronic systems, but they modeled an abnormal sequence in addition to a genuine sequence. This is interesting since assessing the abnormality of outliers detected decreases the number of false positives. Not only do they flag outliers as such, but they also confirm that these outliers belong to the abnormal class.
\newline

\cite{lucas2019ter} proposed a protocol for creating sequential features describing temporal dependencies between the transactions in a supervised way. For this purpose they isolated 3 different perspectives. First, they compare the likelihood of sequence of transactions against historical honest sequences of transactions but also against historical sequences containing at least one fraudulent transaction. The reasoning is that to have a risk of fraud it is not enough for a sequence to be far from honest behavior but it is also expected to be relatively close to a risky behavior. Second, they create these features in order to describe the card-holders sequences but also the merchants sequences. Credit card transactions can be seen as edges in a bipartite graph of merchants and card-holders and they showed that the efficiency of detection was better when taking into account both card-holder and terminal sequences. Third, they considered the elapsed time between two transactions and the amount of a transaction as signal for building their features.
Combinations of the three binary perspectives give eight sets of sequences from the (training) set of transactions. An Hidden Markov Model is trained on each one of these sequences. The proposed features is the likelihood that each the sequence of recent transactions have been generated by each HMM. They have shown a 15\% increase in detection compared to the transactions aggregation strategy of \cite{whitrow2008}.

\subsubsection{Recurrent neural networks}

Recurrent neural networks (RNNs) are a class of artificial neural networks where connections between nodes form a directed graph along a temporal sequence. These connections allow to model temporal dynamic behaviors.
\newline

Recurrent neural networks  can be considered as an hidden state based model where the inner layers reproduce the hidden state based architectures. Furthermore, the connectivity between consecutive hidden states (or nodes) also follows the Markov property ($p(y_{t})$ only depends on $y_{t-1}$). RNNs are discriminant models that aims to predict the label of an event given a past sequence of events. They are detailed extensively in \cite{graves2012}.
\newline

HMMs for sequential anomaly detection have mostly been used in order to raise an anomaly score using the likelihood of the sequence to have been generated by the HMM. However, recurrent neural network (RNN) approaches and more precisely long short term memory networks were more varied. Similarly to the previously mentioned HMM approaches, LSTM were used in order to raise anomaly scores by \cite{bontemps2016} for an intrusion detection task. The anomaly scores of a window of events is the distance between this window of events and the window of events that would have been generated by the LSTM for the corresponding user. Since LSTM are discriminative models, a custom distance has been created by the authors for this generative task. \cite{malhotra2015} modeled the genuine distribution of sequence based dataset with LSTM in order to afterwards raise an anomaly score for electrocardiograms, space shuttle valve time series and power demand datasets. Alternatively, \cite{ergen2017} used LSTM to extend various length sequences in order to obtain sequences of equal length for raising an anomaly score afterwards with the help of One Class SVM.

\cite{wiese2009} and \cite{jurgovsky2018} modeled the sequence of transactions in a card holder's account with Long short-term memory networks. In contrast to generative sequential models, the LSTM-based approaches are trained discriminatively by jointly learning a map from the sequence of transactions to a sequence of latent states and a map from latent states to the binary output variables (i.e. the fraud label of a transaction). The applications of LSTMs to credit card fraud detection are based on the assumption that the succession patterns of transactions bear relevant information for detecting frauds among them. Each latent state is considered a "summary" of the sequence of previously issued transactions and thereby it might permit a better classification of a recent purchase in the light of its summarized purchase sequence. However, in their comparison with order-agnostic feature aggregations, Jurgovsky et al.\cite{jurgovsky2018} found that feature aggregations yielded the best results overall - even in the card-present scenario in which consistent purchase patterns are more prevalent.

LSTM were also used in order to raise anomaly scores by \cite{bontemps2016} for an intrusion detection task. The anomaly scores of a window of events is the distance between this window of events and the window of events that would have been generated by the LSTM for the corresponding user. Since LSTM are discriminative models, a custom distance has been created by the authors for this generative task. \cite{malhotra2015} modeled the genuine distribution of sequence based dataset with LSTM in order to afterwards raise an anomaly score for electrocardiograms, space shuttle valve time series and power demand datasets. Alternatively, \cite{ergen2017} used LSTM to extend various length sequences in order to obtain sequences of equal length for raising an anomaly score afterwards with the help of One Class SVM.

\subsection{Counting-based approaches}

Sequential anomaly detection aims to raise an anomaly score for events, subsequence of events or whole sequence with respect to a genuine distribution.
\newline

For sequential tasks, multiple definitions of anomalies exist that need to be handled differently (\cite{chandola2012}):
\begin{itemize}
\item An event of a sequence may be anomalous. For example, when a transaction is fraudulent within the sequence of transactions of a given card-holder.
\item A subsequence may be anomalous. For example, if there is a delay between the credit card theft and its reporting, the transactions issued with the credit card are not made by the card-holder but by the fraudster and must be identified as anomalous.
\item A sequence may be anomalous with respect to a set of non anomalous sequences. For example, a sentence in french will be anomalous in a set of english sentences.
\end{itemize}

One simple way to raise an anomaly score is counting the likelihood of the event or fixed length window of events in the genuine set of events.
\begin{equation}
L_{i} = \frac{\# subsequence(i)}{\# subsequences}
\label{eq:oddratio}
\end{equation}

\cite{forrest1999} proposed the t-STIDE technique (threshold based sequence time delay embedding) in order to raise anomaly scores for sequences. It consists in computing the ratio of sliding windows that have a likelihood (see equation \ref{eq:oddratio}) below a certain threshold $\lambda$. The anomaly score of the test sequence is proportional to the number of anomalous windows in the test sequence. The corresponding formula is:
\begin{equation}
A_{Sq} = \frac{|i: L(w_{i})< \lambda , 1 \le i \le t |}{t}
\end{equation}

This kind of window based counting can also be used as anomaly signal to calculate a 2nd order anomaly score. For example, \cite{hofmeyr1998} proposed to compute the average of the strength of the anomaly signal over a set of windows that cover the considered event or subsequence.
\begin{equation}
A(Sq) = \frac{1}{N} \sum \limits_{i=1}^N A(w_{i})
\end{equation}

\section{Conclusion} \label{sec:conclusion}

Credit card fraud detection presents several intrinsic challenges. (1) There is a strong difference in terms of absolute numbers between the positive and the negative class: fraudulent transactions are usually much more rare than genuine transactions. (2) Purchase behaviours and fraudster strategies may change over time, making a learnt fraud detection decision function irrelevant if not updated. This phenomenon named dataset shift (change in the distribution $p( x, y )$ ) may prevent fraud detection systems from obtaining good performances. (3) The feature set describing a credit card transaction usually ignores detailed sequential information which was proven to be very relevant for the detection of credit card fraudulent transactions.

Different approaches for tackling each of these challenges are highlighted in this survey and for each of these approaches at minimum one research work is described in full detail. The goal is to provide the reader with useful information on all the different research subjects introduced here.

\bibliographystyle{plainnat}


\section{Introduction}

At present, large-scale image datasets play a critical role in the deep learning area. ImageNet\cite{imagenet_cvpr09}, COCO\cite{COCO}, PASCAL VOC\cite{VOC} actually brought modern deep learning technology into new stage.

Bird species classification is a difficult problem not only because the similarity between different species of birds, but also because they have more than one type of plumage even for the same bird in different timings or areas. Besides, pictures of the birds usually show different poses and actions (e.g., birds in the water, in-flight or perching).

We present DongNiao International Birds 10000 (DIB-10K), a dataset contains not only different types but also a large number of different poses and gestures of birds. By aggregating tremendous bird images, this dataset would push the limits of the visual abilities for machine learning technologies. Currently, the whole dataset could be download from DongNiao website\footnote{\url{http://ca.dongniao.net/download.html}}.
\section{Dataset Specification and Collection}

DongNiao International Birds 10000 contains 4,876,536 thumbnails of 10,922 different bird species. It complies to IOC 10.1 taxonomy\cite{IOC10}, and covers all species of the birds in the world. The size of each thumbnail is 300 x 300 pixels, which contains only one object and it is in the centre of the thumbnail.

Images were crawled from the internet and then filtered by human reviewing. After that, all images were processed through a bird-detection tool written by using Tensorflow\cite{abadi2016tensorflow} with Faster R-CNN\cite{FasterRCNN} model. Each bird in images will be cropped out and resized to fit in a square of 300 x 300 pixels thumbnail.

The thumbnails were not stretched to the square directly but only one edge would be expanded or shrink to 300 pixels to keep the birds' aspect ratio. This step might cause black "frames" in the thumbnails, as Fig ~\ref{fig:samples} shows.

To get the correct labels or species of images from the internet, we used the image search engines to search for images with the scientific names and/or common names of each species as the keys. Images were filtered out by checking the title of the pictures. Thus all the images without scientific names or common names are dropped.

Three reasons make the DIB-10K dataset special from those popular image datasets such as ImageNet\cite{imagenet_cvpr09} and COCO\cite{COCO}.

{\bf Bird centred:} For every thumbnail, the bird is in the central position. This makes it much easier for a classification job but not suitable for object detection job.

{\bf Large number of categories:} DIB-10K has more than 10000 categories, which is far larger than CUB-200\cite{CUB200}'s 200 species, and ILSVRC\cite{ILSVRC15}'s 1000 categories. It is the biggest bird dataset in the status quo.

{\bf Severely unbalanced:} Different bird species have different population. Also, some birds are more difficult than others to take pictures of. Hence the different categories of DIB-10K have different numbers of pictures, as Fig ~\ref{fig:dist} shows. For example, the {\it Western Jackdaw} has more than four thousand images but the {\it Mauritius Blue Pigeon}, which has already extinguished, has only 20 images. This brings a big challenge to build a machine learning classification system.

\begin{table}
  \begin{center}
    \begin{tabular}{|l|c|c|}
      \hline
      Dataset & No. of images & No. of categories \\
      \hline\hline
      ImageNet & 14 million & 20,000 \\
      COCO & 328,000 & 91 \\
      PASCAL VOC & 20,000 & 20 \\
      CUB-200 & 6033 & 200 \\
      DIB-10K & 4,800,000 & 10,000 \\
      \hline
    \end{tabular}
  \end{center}
  \caption{Compare our dataset with other popular ones}
\end{table}

\begin{figure*}
  \caption{Randomly choose 65 categories to visual the distribution of numbers of each species}
  \begin{center}
    \includegraphics[width=\linewidth]{birds_dist}
  \end{center}
  \label{fig:dist}
\end{figure*}

\section{Conclusion}
DIB-10K has a total of 4,876,536 images over 10,922 bird species. It includes all the bird species around the world and hence become an interesting and challenging image dataset for fine-grain classification on ornithology and machine learning areas.

\begin{figure*}
  \caption{Example of five bird images from each of random categories}
  \begin{center}
    \includegraphics[width=\linewidth, height=\textheight]{18birds_0}
  \end{center}
  \label{fig:samples}
\end{figure*}

\begin{figure*}
  \begin{center}
    \includegraphics[width=\linewidth, height=\textheight]{18birds_1}
  \end{center}
\end{figure*}

{\small
\bibliographystyle{ieee_fullname}


\subsubsection{#1}}
\newcommand{\shead}[1]{\paragraph*{\textbf{#1}:}}

\graphicspath{{../Figs/}}

\begin{document}

\title{Evidence of the Berezinskii-Kosterlitz-Thouless Phase in a Frustrated Magnet}

\author{Ze~Hu}
\thanks{These authors contributed equally to this study.}
\affiliation{Department of Physics and Beijing Key Laboratory of
Opto-electronic Functional Materials $\&$ Micro-nano Devices, Renmin
University of China, Beijing, 100872, China}

\author{Zhen~Ma}
\thanks{These authors contributed equally to this study.}
\affiliation{National Laboratory of Solid State Microstructures and Department
of Physics, Nanjing University, Nanjing 210093, China}
\affiliation{Institute for Advanced Materials, Hubei Normal University, Huangshi 435002, China}

\author{Yuan-Da~Liao}
\thanks{These authors contributed equally to this study.}
\affiliation{Beijing National Laboratory for Condensed Matter Physics, and Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China}
\affiliation{School of Physical Sciences, University of Chinese Academy of Sciences, Beijing 100190, China}

\author{Han~Li}
\thanks{These authors contributed equally to this study.}
\affiliation{School of Physics, Key Laboratory of Micro-Nano Measurement-Manipulation and
Physics (Ministry of Education), Beihang University, Beijing 100191, China}

\author{Chunsheng~Ma}
\affiliation{Department of Physics and Beijing Key Laboratory of
Opto-electronic Functional Materials $\&$ Micro-nano Devices, Renmin
University of China, Beijing, 100872, China}

\author{Yi~Cui}
\affiliation{Department of Physics and Beijing Key Laboratory of
Opto-electronic Functional Materials $\&$ Micro-nano Devices, Renmin
University of China, Beijing, 100872, China}

\author{Yanyan~Shangguan}
\affiliation{National Laboratory of Solid State Microstructures and Department
of Physics, Nanjing University, Nanjing 210093, China}

\author{Zhentao~Huang}
\affiliation{National Laboratory of Solid State Microstructures and Department
of Physics, Nanjing University, Nanjing 210093, China}

\author{Yang~Qi}
\email{qiyang@fudan.edu.cn}
\affiliation{Center for Field Theory and Particle Physics, Department of Physics, Fudan University, Shanghai 200433, China}
\affiliation{State Key Laboratory of Surface Physics, Fudan University, Shanghai 200433, China}
\affiliation{Collaborative Innovation Center of Advanced Microstructures, Nanjing University, Nanjing 210093, China}

\author{Wei~Li}
\email{w.li@buaa.edu.cn}
\affiliation{School of Physics, Key Laboratory of Micro-Nano Measurement-Manipulation and Physics (Ministry of Education),
 Beihang University, Beijing 100191, China}
\affiliation{International Research Institute of Multidisciplinary Science, Beihang University, Beijing 100191, China}

\author{Zi~Yang~Meng}
\email{zymeng@hku.hk}
\affiliation{Department of Physics and HKU-UCAS Joint Institute of Theoretical and Computational Physics, The University of Hong Kong, Pokfulam Road, Hong Kong SAR, China}
\affiliation{Beijing National Laboratory for Condensed Matter Physics, and Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China}
\affiliation{Songshan Lake Materials Laboratory, Dongguan, Guangdong 523808, China}

\author{Jinsheng~Wen}
\email{jwen@nju.edu.cn}
\affiliation{National Laboratory of Solid State Microstructures and Department
of Physics, Nanjing University, Nanjing 210093, China}
\affiliation{Collaborative Innovation Center of Advanced Microstructures, Nanjing University, Nanjing 210093, China}

\author{Weiqiang~Yu}
\email{wqyu\_phy@ruc.edu.cn}
\affiliation{Department of Physics and Beijing Key Laboratory of
Opto-electronic Functional Materials $\&$ Micro-nano Devices, Renmin
University of China, Beijing, 100872, China}

\begin{abstract}
\end{abstract}

\date{\today}
\maketitle

\noindent {\bf{Abstract}}\\
\textbf{
The Berezinskii-Kosterlitz-Thouless (BKT) mechanism,
building upon proliferation of topological defects in 2D systems,
is the first example of phase transition beyond the Landau-Ginzburg paradigm
of symmetry breaking. Such a topological phase transition has long been
sought yet undiscovered directly in magnetic materials.
Here, we pin down two transitions that bound a BKT phase in an ideal 2D
frustrated magnet TmMgGaO$_4$,
via nuclear magnetic resonance under in-plane magnetic fields,
which do not disturb the low-energy electronic states and allow
BKT fluctuations to be detected sensitively.
Moreover, by applying out-of-plane fields, we find a critical scaling
behaviour of the magnetic susceptibility expected
for the BKT transition. The experimental findings can be explained
by quantum Monte Carlo simulations applied
on an accurate triangular-lattice Ising model of the compound
which hosts a BKT phase. These results provide
a concrete example for the BKT phase
and offer an ideal platform for future investigations
on the BKT physics in magnetic materials.
}

\noindent {\bf{Introduction}}\\
Topology
plays an increasingly important role in understanding different phases and phase transitions
in correlated quantum matters and materials.
One prominent example is the BKT mechanism in 2D systems
\cite{Berezinskii1971,Berezinskii1972,Kosterlitz1972,Kosterlitz1973,Kosterlitz1974},
which is associated with the binding and unbinding of topological defects.
The BKT transition cannot be characterized by conventional order parameters,
and constitutes the earliest example of phase transition beyond the Landau-Ginzburg
paradigm of spontaneous symmetry breaking.
Historically, the BKT mechanism was introduced in the $XY$ spin model and long predicted
to occur in magnetic thin films~\cite{Berezinskii1971,Kosterlitz1973}.
In experiments, signatures of the BKT transition have been observed in superfluid helium films~\cite{Bishop1978},
as well as 2D superconductor films~\cite{Hebard1980,Epstein1981} and arrays~\cite{Resnick1981}.
Regarding the original proposal in layered $XY$-type magnets,
despite intensive efforts~\cite{Heinrich2003,Cuccoli2003,Wawrzy2008,Wheeler2009,Tutsch2014,Kumar2019},
direct and unambiguous observation of the BKT transition is still lacking.
One major obstacle is the 3D couplings in the magnets,
although weak, will inevitably enhance the confining potential of vortices~\cite{Kumar2019}, leading to 3D ordering that masks the BKT transition.
Therefore, it is of fundamental interest to find and identify BKT materials that
could overcome the obstacle and study the topology-related low-energy dynamics.

\begin{figure}[t]\includegraphics[width=8.5cm]{fig1.pdf}
\caption{{\bf Phase diagram of TmMgGaO$_4$ under out-of-plane magnetic fields.} Under zero field, there are paramagnetic (PM), BKT, and antiferromagnetic (AFM) phases. The $T_{\rm U}$ ($T_{\rm L}$) is the upper (lower) BKT transition temperature, determined from the plateau structure in the NMR spin-lattice relaxation rate $1/^{69}T_1$ (see Fig.~\ref{fig:fig2}{\bf c} for details). The BKT phase between $T_{\rm U}$ and $T_{\rm L}$ is illustrated by the solid vertical line, while the AFM regime is indicated by the arrow. The contour background depicts the magnetic specific heat $C_{\rm m}/T$ at various fields and temperatures, with data adapted from Ref.~\cite{Li2020} and plotted in logarithmic scale. $T^*$ corresponds to the maximum of $C_{\rm m}/T$ at each field, signifying the position of strong magnetic fluctuations. $T'$($T''$) denotes the temperature at a specific field where a peak is found in the differential susceptibility $dM/dH$, shown in Fig.~\ref{fig:fig3} {\bf b}. The Curie-Weiss temperature $\theta_{\rm CW}$ is obtained from the $1/^{69}T_1T$ (see Supplementary Fig.~1). Remarkably, $T^*$, $T'$, $T''$ and $\theta_{\rm CW}$ all collapse to the same phase boundary between the BKT-like regime and AFM phase. A magnetically ordered phase is supposed to lie below the dome-like boundary. Errors represent one standard deviation throughout the paper.}\label{fig:fig1}
\end{figure}

\begin{figure*}[t]
	\includegraphics[width=12cm]{fig2.pdf}
	\caption{\label{fig:fig2}
	{\bf NMR spectra and spin-lattice relaxation rates of TmMgGaO$_4$.}
	{\bf a} $^{69}$Ga NMR spectra at different
	temperatures under a 3~T in-plane field,
	with zero frequency corresponding to
	${\gamma}H$=~30.692 MHz.
        Data are shifted vertically for clarity.
        At high temperatures, the spectra are roughly symmetric,
        while for $T \leq $ 2~K, a shoulder-like structure
        can be resolved on the left side of the main peak.
	{\bf b} NMR hyperfine shift
	$^{69}K_n = ({\bar{f}}/{\gamma H} - 1) \times 100\%$
	as a function of temperature,
	where $\bar{f}$ is the average frequency of each spectrum.
	{\bf c} NMR spin-lattice relaxation rate $1/^{69}T_1$
	versus temperature measured under in-plane fields of 3~T and 1~T.
	A plateau-like feature, characterizing strong magnetic fluctuations,
	is observed between $T_{\rm L} \simeq 0.9$~K (lower BKT transition)
	and $T_{\rm U} \simeq1.9$~K (upper BKT transition).
	{\bf d} $1/T_1$ data computed from the dynamical
	spin-spin correlation function with contributions from
	all momentum points
	[c.f., Eq.~(\ref{Eq:Ahf}), and see the hyperfine form factor
	in Supplementary Note~\B{2}]
	in the Brillouin zone (left scale)
	and from $\rm K'$ point in the vicinity of the K point (right scale),
	through large-scale
	QMC simulations (see Methods).
	}
\end{figure*}

Recently, a layered frustrated rare-earth antiferromagnet TmMgGaO$_4$~\cite{Cava2018,Shen2019,Li2020}
was reported to ideally realize the triangular-lattice quantum Ising (TLI) model~\cite{Li2020a}.
The relatively large interlayer distance of 8.3774~\r{A} along the $c$ axis
gives rise to excellent two dimensionality~\cite{Cava2018}, and no sign of conventional
3D phase transition was observed in either specific heat or magnetic susceptibility measurements.
Nevertheless, it was reported from neutron scattering that
TmMgGaO$_4$ ordered below $\sim$1~K into an antiferromagnetic phase
with a six-fold symmetry breaking~\cite{Shen2019,Li2020}, which closely resembles the ground state
of the TLI model originated from an order-by-disorder mechanism~\cite{R.Moessner2001,Isakov2003}.
At higher temperatures, the effective $XY$ degrees of freedom emerge,
and the BKT mechanism is expected to come into play~\cite{Isakov2003}.

In TmMgGaO$_4$, each Tm$^{3+}$, with a $4f^{12}$
electron configuration and a spin-orbit moment $J=6$,
forms a non-Kramers doublet due to the crystal-electric-field splitting.
The doublet is well separated from the rest levels by about 400~K~\cite{Li2020},
and can thus be regarded as an effective spin-1/2.
There is further a fine energy splitting within the doublet,
induced by the local trigonal crystal field~\cite{Cava2018},
acting as an intrinsic ``transverse field" applied on the effective spin.
From the magnetization measurements~\cite{Cava2018,Shen2019,Li2020},
one observes that Tm$^{3+}$ ions contribute highly anisotropic
Ising-type moments with $J_z=\pm 6$ along the $c$ axis,
resulting in an effective out-of-plane $g$-factor  $\sim$13.2~\cite{Li2020,Li2020a}.
On the contrary, the effective in-plane $g$-factor and
dipolar moment in the $ab$ plane are negligible.

Facilitated by this feature in TmMgGaO$_4$, in this work
we employed nuclear magnetic resonance (NMR),
a sensitive low-energy probe, to detect the BKT phase.
We applied a moderate in-plane field of 3~T,
which is adequate to collect the $^{69}$Ga NMR signals and,
at the same time, hardly disturbs the low-energy electronic states of the material.
This is important, since in the TLI model that is believed
to accurately model TmMgGaO$_4$~\cite{Li2020a},
the BKT phase can be fragile against out-of-plane fields
~\cite{Liu2019arXiv,Damle2015,Biswas2018},
thus posing a challenge to NMR measurements.
Taking advantage of the fact that in-plane moment in TmMgGaO$_4$
is mostly multipolar~\cite{Li2020,Li2020a},
our NMR experiments with in-plane fields, which merely couples to the nuclear spins,
can clearly identify the BKT phase in the material.

As shown in Fig.~\ref{fig:fig1},
from our NMR measurements of the spin-lattice relaxation rate $1/T_1$,
we identify $T_{\rm U} \simeq 1.9$~K and $T_{\rm L} \simeq 0.9$~K
which represent the upper- and lower-BKT transitions,
where a critical BKT phase
resides at zero magnetic field in between the high-$T$
paramagnetic and low-$T$ antiferromagnetic phases.
This finding is further substantiated by our scaling analysis
of the measured susceptibility data near $T_{\rm L}$,
as well as the simulated NMR and susceptibility data using large-scale
quantum Monte Carlo (QMC) calculations.

\bigskip \noindent {\bf Results}\\
\textbf{NMR probe of the BKT phase.}
The obtained NMR spectra with an in-plane magnetic field $\mu_0 H = 3$~T
are shown in Fig.~\ref{fig:fig2}{\bf a} at representative temperatures.
In order to better resolve the magnetic transition,
the hyperfine shifts $^{69}K_{\rm n}$ of the spectra were analyzed
and plotted in Fig.~\ref{fig:fig2}{\bf b} as a function of temperature.
Upon cooling, $^{69}K_{\rm n}$ peaks at about 0.8-0.95~K,
and then starts to drop at lower temperatures. Therefore, the
ordering temperature is determined to be $T_{\rm L}~{\simeq}~0.9$~K,
consistent with neutron scattering experiments~\cite{Shen2019,Li2020}.
Additionally, both the second moments (width of the NMR spectra)
and the third moments (asymmetry of the spectra)
of the spectra change dramatically below $\sim$~2~K,
suggesting the onset of local hyperfine fields
enhanced by the static or quasi-static magnetic ordering
(Supplementary Fig.~4). These two characteristic temperatures signal
the two-step melting of magnetic order
through two BKT transitions, suggesting an
intermediate floating BKT phase in the system. We suspect that there is some inhomogeneity of the local hyperfine fields, which
is very likely caused by the quenched disorder from
Mg/Ga site mixing~\cite{Li2020}, although
no significant influence on the electronic
and more importantly the magnetic
properties is seen (see more detailed discussions in
Supplementary Note~\B{2}).

The spin-lattice relaxation rate $1/T_1$ provides a highly sensitive detection
of low-energy spin fluctuations~\cite{Julien2000,Kitagawa2008,
Koutroulakis2015,Jeong2017,Moriya2003},
and thus the BKT transition. In Fig.~\ref{fig:fig2}{\bf c} we show the $1/^{69}T_1$
obtained under in-plane fields $\mu_0H = 3$~T  and $1$~T,
which reflects intrinsic spin fluctuations with zero out-of-plane field.
At 3~T, $1/^{69}T_1$ first decreases upon cooling from 10~K
then suddenly increases below $T_{\rm U}\approx$ 1.9~K, indicating
the onset of strong low-energy spin fluctuations. The data at 1~T show similar behaviours.
Below $T_{\rm L}\simeq0.9$~K, $1/^{69}T_1$ drops sharply, consistent with the onset
of the magnetic ordering as also inferred from the hyperfine shift.
Here, $1/^{69}T_1$ is dominated by the gapped
long wavelength excitations in the order state.
At the magnetic phase transition, a peaked feature
in $1/T_1$ develops, caused by the gapless low-energy
spin fluctuations with diverging correlation length.
Remarkably, at temperatures between
$T_{\rm U}\simeq1.9$~K and $T_{\rm L}\simeq0.9$~K, $1/^{69}T_1$ exhibits a
plateau-like structure, indicating the emergence of
a highly fluctuating phase with diverging spin correlations yet no true long-range order,
which is the hallmark of a BKT phase
~\cite{Berezinskii1971,Berezinskii1972,Kosterlitz1972,Kosterlitz1973,Kosterlitz1974}.
Therefore, it is for the first time that such
a phase is unambiguously observed in a magnetic crystalline material.

Our unbiased QMC simulations on the TLI model of the material (see Methods),
with accurate coupling parameters determined in Ref.~\cite{Li2020a},
quantitatively justifies the existence of the
BKT phase between $T_{\rm L}$ and $T_{\rm U}$.
We computed $1/T_1$ and compared with the experiment below.
Figure~\ref{fig:fig2}{\bf d} shows the calculated $1/T_1$ by including
fluctuations from all momentum points in the Brillouin zone
(cf. Supplementary Fig.~2), and compare to that from $\rm K'$
(around the K point at the corner of Brillouin zone).
The former shows a decreases upon cooling below 4~K,
and then an upturn above
$T_{\rm U}^* \simeq 2$~K, followed by a rapid decrease
below $T_{\rm L}^* \simeq 1$~K.
These behaviours are in excellent agreement with the
measured $1/^{69}T_1$.
The latter reflects gapless excitations of
the $XY$ degrees of freedom emergent in the BKT phase,
where the calculated $1/T_1$ from $\rm K'$ exhibits
an anomalous increase down to $T_{\rm U}^*$,
below which the increment slows down.
The contribution to $1/T_1$ near the K point reaches a
maximum at $T_{\rm L}^*$ and drops rapidly below it.
The absence of critical spin fluctuations at momentum
away from the K point suggests that the measured
$1/^{69}T_1$ below 2~K is mainly contributed by
excitations around the K point (see Supplementary Note~\B{3}).

Overall, the quasi-plateau behaviours in the QMC results,
as well as the two characteristic temperature scales,
are in full consistency with the NMR measurements.
This constitutes both strong support for the accurate
quantum many-body modeling of the material TmMgGaO$_4$
and also solid proof of the BKT phase therein detected by NMR.
Nevertheless, we note that there are still subtle
differences between the experimental and numerical data.
Needless to say, the real material is always more complicated
than our theoretical minimal model.
For example, influences from higher CEF levels
above the non-Kramers doublet,
the interlayer couplings not included in our model calculations,
and the lack of knowledge on the precise local
hyperfine coupling constant, etc.,
may explain the difference remained between the
panels (c) and (d) of Fig.~\ref{fig:fig2}.

\begin{figure*}[t]
	\includegraphics[width=18cm]{fig3.pdf}
	\caption{\label{fig:fig3}
	{\bf Uniform magnetic susceptibility of TmMgGaO$_4$ and scaling analysis.}
	{\bf a} dc susceptibility $\chi$ as functions of temperatures
		under small out-of-plane ($H~\ensuremath{/\kern -0.8em /}\xspace~c$) and in-plane ($H~\ensuremath{/\kern -0.8em /}\xspace~ab$) fields.
		The latter is multiplied by a factor of 20 for visualizing purpose.
		The deviation of data below 2~K indicates the entry to
		the BKT phase and the field-suppression of magnetic correlations.
		{\bf b} The differential susceptibility ${dM}/{dH}$
		under out-of-plane fields at different temperatures.
		The kinks at low fields, as denoted by the down arrows,
		suggest the transition from the BKT-like phase to the
		ordered phase (under the dome
		in Fig.~\ref{fig:fig1}) with increasing fields.
		The peaked features at high fields suggests a
		quantum phase transition to the polarized phase.
		{\bf c} Fits of ${dM}/{dH}$ to the power-law scaling
		function $dM/dH\sim H^{-\alpha}$ with
		$\alpha=2/3$ for the 0.4 and 0.8-K data,
		and $\alpha=0.123$ for 2.1~K data in the log-log scale.
		The 3~K data follow the $\alpha=0$ line in the
		paramagnetic phase.
		{\bf d}  ${dM}/{dH}$ by the QMC calculations in the same
		field and temperature range as in {\bf c}, and
                 fits to the power-law function with exponents $\alpha$,
		which give consistent results as experiments.
	}
\end{figure*}

\noindent\textbf{Universal magnetic susceptibility scaling.}
Magnetic susceptibility $\chi$ measurements are also performed to
strengthen the finding of the BKT phase.
In Fig.~\ref{fig:fig3}{\bf a}, we show the overall
temperature dependence of $\chi$
at different out-of-plane fields. For $T\gtrsim2$~K, $\chi$ increases monotonically
upon cooling, and barely changes with fields. However, for $T\lesssim 2$~K,
approximately the upper BKT transition $T_{\rm U}$ as obtained
from the $1/^{69}T_1$ measurements, $\chi$ increases as the field decreases,
suggesting the onset of peculiar magnetic correlations.
With further cooling, the susceptibility gets flattened with temperature.
The magnetization $M(H)$ is further measured at selected temperatures
(data shown in Supplementary Fig.~6),
and for the sake of clarity, the differential susceptibility ${dM}/{dH}$ is plotted
in Fig.~\ref{fig:fig3}{\bf b}. At around 2.5~T, a pronounced peak can be
observed at low temperature, indicating the existence
of a quantum phase transition, and the phase at lower fields
should be a magnetically order phase,
although its precise nature remains to be uncovered.
Besides the high-field feature, for temperatures at 0.8~K and above,
a kinked feature is clearly resolved on each ${dM}/{dH}$ curve at low fields;
whereas at 0.4~K, the low-field kink disappears,
which posts a question of whether there
is a quantum transition or merely a crossover from the zero-field
AFM phase to the finite-field ordered phase
under the dome in Fig.~\ref{fig:fig1}.
The temperature and field values indicated by the
down arrows in Fig.~\ref{fig:fig3} are denoted as
$T'$ and $T''$ in the phase diagram (Fig.~\ref{fig:fig1}).

Field-theoretical analysis of the TLI model~\cite{Damle2015,Biswas2018}
has predicted that upon applying a small out-of-plane field,
the differential susceptibility $dM/dH$ shall exhibit
a divergent power-law behaviour as ${dM}/{dH}\sim H^{-\alpha}$
in proximity to the BKT phase.
At $T_{\rm L}$, $\alpha$ has the value of ${2}/{3}$,
which corresponds to a critical exponent $\eta=1/9$ at the lower BKT transition
and is originated from the six-fold symmetry breaking~\cite{Damle2015}.
The exponent $\alpha$ gradually decreases as temperature increases,
and above an intermediate temperature between $T_{\rm L}$ and $T_{\rm U}$,
$\alpha = 0$ due to non-universal contributions.
This is exactly what we observe in Fig.~\ref{fig:fig3}{\bf c}.
We fit the $dM/dH$ with the power-law function at different temperatures,
with the fitting regime chosen between 0.6-0.9~T. At 0.8~K,
$\alpha$ is very close to the expected value of $2/3$ (and thus $\eta=1/9$) at $T_{\rm L}$,
which constitutes a remarkable fingerprint evidence for the BKT transition.
At lower temperatures such as 0.4~K, the exponent is also close to $2/3$ because the
susceptibility saturates with temperature, as shown in Fig.~\ref{fig:fig3}{\bf a}.
At high temperatures, $\alpha$ drops rapidly to a small value 0.12 at 2.1~K,
and becomes effectively zero at 3.0~K.

Therefore, the
susceptibility scaling gives the lower BKT transition at about 0.8~K and
upper transition probably between 2.1 and 3~K, in good agreement
with the $T_{\rm L}$ and $T_{\rm U}$ estimated from NMR.
These results are also fully consistent with
our QMC calculations on the susceptibility shown in
Fig.~\ref{fig:fig3}{\bf d}. At $T_{\rm L}$ or lower,
$\alpha$ is ${2}/{3}$, then decreases to a very small exponent $0.086$
at 2.67~K, and above 3~K, becomes zero within numerical uncertainty.
Such a power-law behaviour in ${dM}/{dH}$ again signifies the finite-temperature
window of the BKT phase with diverging magnetic correlations,
which gives rise to the universal power-law scaling of magnetic susceptibility.

\bigskip \noindent {\bf Discussion}\\
We believe the findings in this work are of various fundamental values.
Since the original proposal of a BKT phase in
magnetic films~\cite{Kosterlitz1972,Kosterlitz1973,Kosterlitz1974},
which also triggered the currently thriving
research field of topology in quantum materials,
tremendous efforts have been devoted to finding the
BKT phase in magnetic crystalline materials,
yet hindered by the obstacle outlined in the Introduction.
Here, benefiting from NMR as a sensitive low-energy probe,
and the nearly zero planar gyromagnetic factor
in a TLI magnet TmMgGaO$_4$, we are able to
reveal two BKT transitions and a critical BKT phase
with an emergent $XY$ symmetry.
Together with the power-law behaviour of the susceptibility,
and excellent agreement between
our QMC simulation and experiment data,
we unambiguously identify the long-sought BKT
phase in a magnetic crystalline material.

Many intriguing questions are stimulated,
based on the phase diagram in Fig.~\ref{fig:fig1}
obtained here. First, for example, what is the nature
of the ordered phase under finite fields, are there
further exotic phases and transitions in the phase diagram,
and will there be a field-induced quantum phase transition
at the high-field side of the dome -- these are all of great interests
to be addressed in future studies.
Second, it should be noted that the dynamical properties
obtained by QMC calculations in \Fig{fig:fig2} are computed on
a large, while finite-size, 36$\times$36 lattice, which already produces
$1/T_1$ data in very good agreement with the experimental measurements.
Such a great agreement is surprising,
given the possible existence of randomness
from Ma/Ga site mixing in the material TmMgGaO$_4$ \cite{Li2020},
and also the lattice disorder revealed by the large high-temperature
second moment of the NMR spectra (Supplementary note~\B{4}).
Although the random distribution in intrinsic transverse
fields and spin couplings does not seem to
alter the low-temperature spin ordered
phase as well as the sharp spin excitation line shapes
\cite{Shen2019,Li2020a}, its intriguing
effects on the finite-temperature phase diagram of TLI
and also the compound TmMgGaO$_4$ is calling for further studies.

Third, in the study of BKT transition in superfluid systems,
it has been observed experimentally and understood theoretically
that additional dissipations also appear above the
transition temperature due to fluctuations of vortices~\cite{Bishop1978}.
Hence, the plateau of $1/T_1$ we observe may also cover
regions slightly above the upper BKT transition temperature.
We leave this subtlety to future numerical and experimental efforts.
Lastly, in general terms, whether there are other rare-earth magnetic
materials in the same family of TmMgGaO$_4$ that,
acquire similar 2D competing magnetic interactions
from highly anisotropic gyromagnetic factor and unique triangular-lattice structures,
and also exhibit the BKT physics, is quite intriguing and calls for future investigations.
All these directions are ready to be explored from here.\\

\bigskip \noindent {\bf Methods}\\
\textbf{Crystal growth and susceptibility measurements.}
Single crystals were grown by the optical-floating-zone
method with an image furnace
(IR Image Furnace G2, Quantum Design).
The natural cleavage surface of the crystals is the $ab$ plane,
which allows us to align the field orientation
within an error of $2^{\circ}$. The dc susceptibility was measured
in a PPMS VSM (Quantum Design) for temperatures above
2~K and in a He-3 MPMS (Quantum Design) for temperatures ranging
from 0.4 to 2~K.  \\

\noindent
\textbf{NMR measurements.}
The $^{69}$Ga ($I$~=~${3}/{2}$, $\gamma$~=~10.219~MHz/T)
NMR spectra were collected with the standard spin-echo sequence,
with frequency sweep by a 50 kHz step using a topping tuning circuit.
The NMR hyperfine shift was obtained by calculating
the change of the first moment of the spectra.
The spin-lattice relaxation rate $1/^{69}T_1$
was measured by the inversion-recovery method, where a $\pi/2$ pulse
was used as the inversion pulse. The NMR data from 1.8~K and above were measured
in a variable temperature insert, and the data from 1.8~K and below were measured in
a dilution refrigerator. The weak NMR signal at low fields and the rapid decrease of
$^{69}T_2$ upon cooling (Supplementary Note~\B{5}) prevented us to measure the $1/^{69}T_1$
for in-plane fields less than 3~T, with temperature below 1.8~K.
Whereas for in-plane fields of 4~T and higher, the sample could not be held
in position because of the large anisotropy in the $g$ factor and unavoidable
sample misalignment ($\lesssim2^\circ$).   At $T=$1.2~K, we did not find any
change of $1/^{69}T_1$ with two different radio frequency excitation levels (14~mT and 24~mT),
and with different frequencies across the NMR line, within the error bar. \\

\noindent
\textbf{Triangular lattice Ising model.}
At zero field, the intralayer couplings
in TmMgGaO$_4$ can be described by the TLI Hamiltonian,
\begin{eqnarray}
\hspace{-9mm}&&H = J_1\sum_{
\langle i,j \rangle}S_i^zS_j^z + J_2\sum_{\langle\langle i,j\rangle\rangle}S_i^zS_j^z + \sum_i\Delta S_i^x,
 \label{h1}
\end{eqnarray}
where $J_1$ and $J_2$ are the superexchange interactions among Tm$^{3+}$,
$\langle i,j \rangle$ and $\langle\langle i,j \rangle\rangle$ refer to the
nearest neighbors and the next-nearest neighbors, respectively,
and $\Delta$ is the energy splitting within the non-Kramers doublet imposed by the crystal field.
We have shown that the parameter set $J_1$$=$ 0.99~meV,
$J_2/J_1$$\approx$~0.05 and $\Delta/J_1\approx$0.54
reproduces the experimental results of the
transition temperatures and the
inelastic neutron scattering spectra~\cite{Li2020a}.

In the TLI model [Eq.~(\ref{h1})],
we can define a complex field $\psi$ as a
combination of the Ising ($Z$) components
$m_{\rm A,B,C}^z$ on three sublattices, i.e.,
\begin{equation}
\psi = m^z_{\rm A} + {\rm e}^{i2\pi/3} \, m^z_{\rm B}  +  {\rm e}^{i 4\pi/3} \, m_{\rm C}^z.
\label{Eq:psi}
\end{equation}
Notably, $\psi = |\psi| {\rm e}^{i \theta}$ is a complex order parameter
that represents the emergent $XY$ degree of freedom
relevant to the BKT physics in the TLI model.\\

\noindent
\textbf{Quantum Monte Carlo calculations.}
QMC simulations were performed in the
path integral in the $S^{z}_{i,\tau}$
basis with discretization in space and time.
The lattice of $L\times L\times L_{\tau}$,
where $L=36$ and $L_{\tau}=\beta/\Delta\tau$
with $\Delta\tau=0.05$ and $\beta=1/T$,
were simulated with both local and Wolff-cluster
update schemes~\cite{YCWang2017,WLJiang2019}.
The $1/T_1$ results were obtained by first
computing the dynamical spin-spin correlation function
$\langle S^{z}_{i}(\tau)S^{z}_{j}(0)\rangle$,
and then acquiring its real-frequency dependence $S(\mathbf{q},\omega)$
from the stochastic analytic continuation~\cite{Sandvik2016}.
We then determined the $1/T_1$ either by
summing the contributions close to momentum
K or over the entire Brillouin zone, as discussed in the
Fig.~\ref{fig:fig2}d of the main text,
\begin{equation}
T^{-1}_1 (\bold{q}) = \frac{1}{L^2}\sum_{\bold q} |A_{\rm hf}({\bold q})|^2
S({\bold q},\omega \to 0),
\label{Eq:Ahf}
\end{equation}
where $A_{\rm hf}({\bold q})$ is the hyperfine coupling form factor
(see Supplementary Note~\B{2}). Similar analyses have been
successfully applied to the QMC computation
of NMR ${1}/{T_1}$ for the spin-1/2 and spin-1
chains~\cite{Sandvik1995,Capponi2019}.\\

\bigskip \noindent \textbf{Data availability}\\
The data that support the findings of this study are available from the corresponding author upon reasonable request.

\bigskip \noindent \textbf{Code availability}\\
All numerical codes in this paper are available upon request to the authors.\\