# Applied Machine Learning to Anomaly Detection in Enterprise Purchase
Processes

A. Herreros-Martínez
Veolia Group
Alicante, Spain
<EMAIL_ADDRESS>&R. Magdalena-Benedicto, J. Vila-Francés, A.J.
Serrano-López
Universidad de Valencia, Research Group IDAL
Valencia, Spain
<EMAIL_ADDRESS>&S. Pérez-Díaz
University of Alcalá de Henares, Research Group ASYNACS,
Madrid, Spain
<EMAIL_ADDRESS>
###### Abstract

In a context of a continuous digitalisation of processes, organisations must
deal with the challenge of detecting anomalies that can reveal suspicious
activities upon an increasing volume of data. To pursue this goal, audit
engagements are carried out regularly, and internal auditors and purchase
specialists are constantly looking for new methods to automate these
processes. This work proposes a methodology to prioritise the investigation of
the cases detected in two large purchase datasets from real data. The goal is
to contribute to the effectiveness of the companies’ control efforts and to
increase the performance of carrying out such tasks. A comprehensive
Exploratory Data Analysis is carried out before using unsupervised Machine
Learning techniques addressed to detect anomalies. A univariate approach has
been applied through the z-Score index and the DBSCAN algorithm, while a
multivariate analysis is implemented with the k-Means and Isolation Forest
algorithms, and the Silhouette index, resulting in each method having a
transaction candidates’ proposal to be reviewed. An ensemble prioritisation of
the candidates is provided jointly with a proposal of explicability methods
(LIME, Shapley, SHAP) to help the company specialists in their understanding.

##  1 Introduction

The Internal Audit department of a company (normally multinationals groups
and/or big-sized entities) is aimed to ensure the correctness and
effectiveness of the entities’ processes, its compliance to the approved
internal policies and to reduce risks in any form that could be presented [1].
In order to achieve this goal, the companies’ internal teams conduct audits
through on a regular basis defined audit engagements. During their missions,
the auditors identify, evaluate and document adequate information to achieve
the objectives of the engagement [2], carrying out interviews with the
auditees and performing a rigorous tracking of evidences supporting the audit
findings.

Currently, auditing still mainly relies on sampling the information
(registers, transactions, etc.) to assess the processes’ compliance during the
audit engagements [3]. Consequently, the so-called sampling-risk makes that
relevant information in the registers/transactions could remain out of the
sampling selection to be reviewed. Additionally, with the growing amount of
data, this traditional approach becomes obsolete, and the sampling risk is
aggravated [4].

Among the business processes, a special interest resides in searching for
anomalies or misbehaviours on purchases. Internal audit and purchase managers
need to prospect, evaluate, and select the methodologies and IT tools capable
of monitoring expenses and discovering relevant information that can highlight
an out-of-policy act or, even, fraud [5, 6]. The goal is to automate processes
within the company that help to prioritize the investigation activities
according to the level of suspicion of any fact.

In this business context, data analytics has been revealed as a key element in
supporting organisations in the challenge of processing and controlling huge
quantities of information. Among the various algorithms available, those
related to machine learning are known to offer good results in finding
relevant insights over structured and unstructured data. Machine learning and
data analytics are changing the auditing approaches [7] and are becoming an
important part of the control toolbox for companies.

A key drawback to full population tests, based on hand-crafted rules, is that
they will likely be able to find only errors, mistakes or deviations which
were already expected or anticipated. Machine learning approaches have the
ability to go beyond the results obtained through hand-crafted rules. They can
be used to find anomalies in large amounts of data [7] allowing a later
identification of errors, mistakes, circumvented processes or even fraud. Most
of the machine learning approaches that have been tested in auditing are based
on supervised learning [8, 9]. In this technique, an algorithm learns to map
input data to known targets (labelled data). In external auditing, due to
reporting regulations, data shares a similar structure across multiple
companies (e.g., financial statements). This fact makes easier to obtain data
labels for training supervised machine learning models. Regarding internal
auditing, on the contrary, labels for processes flagging fraudulent/non-
fraudulent cases are not available most of the time, given the huge amount of
data or the impossibility of previously labelling the suspicious or fraudulent
cases. In essence, traditional internal audit planning uses a risk-based
approach, selecting areas, departments, or processes to be audited according
to criteria such as the turnover in a business area or how long it has not
been audited. The unsupervised algorithms offer a promising choice in this
context. Within this algorithm’s category, and among others, clustering
techniques have been widely studied in anomaly detection [10, 11].

The goals of the present work are:

  1. i

investigate and propose a methodology able to detect anomalies in an actual
unlabelled purchase dataset,

  2. ii

prioritise the cases to be analysed by the organisation’s specialists, and

  3. iii

provide them an adequate explicability technique of such cases to support
their investigation.

##  2 Purchase Business Process

Most businesses rely on other businesses for products and services that help
them with their operations. The purchasing process enables an organization to
evaluate these business-to-business transactions for efficiency and track
their spending. It must be distinguished from the procurement process.
Procurement focuses on strategies, such as negotiations and researching while
purchasing focuses on the actual act of buying products and services. Despite
the control measures, an inadequate purchase process management can overlook
more deliberate harmful practices (i.e., purchase fraud) and causing
enterprise losses [12]). An effective purchasing process can help prevent
theft, fraud, or irregular spending since it requires documenting all business
transactions.

##  3 Anomalies detection

Anomaly detection is the process of finding outliers in a given dataset.
Outliers are the data objects that stand out amongst other data objects and do
not conform to the expected behaviour in a dataset. Anomaly detection
algorithms have broad applications in business, scientific, and security
domains where isolating and acting on the results of outlier detection is
critical [13]. For the identification of anomalies, there are many different
algorithms, such as classification, regression, and clustering. In case of
having elements with known anomalous behaviours, then supervised algorithms
can be used for anomaly detection. In addition to supervised algorithms, there
are specialized (unsupervised) algorithms whose whole purpose is to detect
outliers without any previous knowledge about the dataset [14].

Anomalies can be classified into three categories [15]:

  * •

Point Anomalies. Whether an individual data instance can be considered as
anomalous with respect to the rest of the data.

  * •

Contextual Anomalies. Whether a data instance is anomalous in a specific
context, but not otherwise. For instance, within a cluster, the points more
distant from the centroid can be considered anomalies within the cluster but
not in the entire dataset.

  * •

Collective Anomalies. Whether a collection of related data instances is
anomalous with respect to the entire data set or not. For instance, all the
elements belonging to a group in a classification problem.

###  3.1 Anomaly Detection Techniques (ADTs)

####  3.1.1 Outlier Detection Using Statistical Methods

Outliers in the data can be identified by creating a statistical distribution
model of the data and identifying the data points that do not fit into the
model or data points that occupy the ends of the distribution tails [16].
Outliers can be detected based on where the data points fall in the standard
normal distribution curve. A threshold for classifying an outlier is
specified. Other statistical techniques take multiple dimensions into account,
computing other distance metrics (multivariate anomaly).

####  3.1.2 Outlier Detection Using Data Science

Normally outliers show a specific set of characteristics that can be exploited
to locate them. According to the unique set of characteristics that identify
an outlier, the detection techniques are:

  * •

Distance-based: By nature, outliers are distant from other data points. If the
average distance of the nearest N neighbours is measured, the outliers will
have a higher value than other normal data points. Distance-based algorithms
utilize this property to identify outliers in the data.

  * •

Density-based: The density of a data point in a neighbourhood is inversely
related to the distance to its neighbours. Consequently, outliers are located
in low-density areas, while the regular data points are concentrated in high-
density areas.

  * •

Distribution-based: Outliers are the data points that have a low probability
of occurrence, and they occupy the tail ends of the distribution curve.

  * •

Clustering: Outliers, by definition, are not similar to normal data points in
a dataset. The outliers can be the lone data points that are not clustered or
the outliers forming a far and low-populated cluster.

  * •

Classification techniques: Whether labelled data is previously available, a
classification model can be built based on a training (labelled) dataset and
tested through the prediction of an unknown dataset (not labelled). The
challenge in using a classification model is the availability of previously
labelled data.

Each of these techniques has multiple parameters and, hence, a data point
labelled as an outlier in one algorithm may not be an outlier in another one.
Thus, it is a good practice to rely on multiple algorithms before labelling
the outliers. This approach is applied in the present work which, considering
the problem statement and the unlabelled input data, will be focused on
clustering (k-Means, DBSCAN) and graph techniques (Isolation Forest).

###  3.2 Algorithms

####  3.2.1 k-Means

k-Means clustering is a distance-based clustering method where the dataset is
divided into k number of clusters. It is one of the simplest and most commonly
used clustering algorithms. In this technique, the user specifies the number
of clusters (k) that the dataset need to be grouped in. The objective of
k-Means clustering is to find a prototype data point (centroid) for each
cluster; all the data points are then assigned to the nearest prototype, which
then forms a cluster. The centre of the cluster can be the mean of all data
objects in the cluster, as in k-Means, or the most represented data object, as
in k-medoid clustering [17]. The evaluation of k-Means clustering has been
done using two methods: i) computing the total of the Sum Squared Errors
(SSE), where good models will have low SSE within the cluster and low overall
SSE among all clusters. SSE can also be referred to as the average within-
cluster distance and can be calculated for each cluster and then averaged for
all the clusters. ii) the Silhouette coefficient, proposed by Peter Rousseeuw
in 1987 [18]. The Silhouette value is a measure of how similar an object is to
its own cluster (cohesion) compared to other clusters (separation). The
Silhouette ranges from -1 to +1, where a high value indicates that the object
is well matched to its own cluster and poorly matched to neighbouring
clusters. If most objects have a high value, then the clustering configuration
is appropriate. If many points have a low or negative value, then the
clustering configuration may have too many or too few clusters.

####  3.2.2 Isolation Forest

The Isolation Forest [19] algorithm takes advantage of two quantitative
properties of anomalies: (i) they are the minority consisting of few
instances, and (ii) they have attribute-values that are very different from
those of normal instances. These characteristics make them more susceptible to
being separated (isolated) from normal instances. The isolation process is
based on a binary tree structure called the isolation tree (iTree). Because of
the susceptibility to isolation, anomalies are more likely to be isolated
closer to the root of an iTree, whereas normal points are more likely to be
isolated at the deeper end of an iTree. The Isolation Forest method (iForest)
builds an ensemble of iTrees for a given dataset; anomalies are those
instances which have short average path lengths on the iTrees [19].
Furthermore, the algorithm has a linear time complexity with a low constant
and a low memory requirement, which works well with high-volume data such as
in the case of the present study [15].

##  4 Materials and Methods

The present work has been carried out using the KNIME Analytics Platform. The
tool KNIME v4.5.2- build March 2022 (Konstanz Information Miner), is a free
and open-source data analytics platform. It allows an easy implementation of
data process and repeatability in case of feedback improvements. KNIME
integrates a wide range of components for machine learning and data mining
through its modular data pipelining "Building Blocks of Analytics" concept
[20]. A graphical user interface and use of JDBC allows assembly of nodes
blending different data sources, including pre-processing (ETL: Extraction,
Transformation, Loading), for modelling, data analysis and visualization
without, or with only minimal, programming.

###  4.1 Data Acquisition & Pre-processing

The data used in the experiment was provided by two companies of a
multinational Group and came from real procurement transactions occurred
during the year 2021. Since this data was related to real transactions, to
avoid any secrecy problems, categorical data that had an associated numerical
identification code were excluded from the dataset. As the file had the
identification codes for all of them, this deletion did not mean any
information loss. Finally, personal information concerning the requester, the
buyer and the approvers of each transaction was received anonymized in order
to respect the personal data protection directives of the European GDPR
regulation. The initial source of data was composed globally of 65.712 records
with 17 columns, split into two datasets, one per company. Table 1 contains
the extended descriptive figures, which further reveals the different
statistical nature of the transactions, inherent in the internal process
organisations within each entity. For example, the number of requesters
placing orders is significantly higher in Company1. On the contrary, the
number of buyers is higher in Company2. This implies a quite different average
of purchase orders and items per requester and buyer.

Table 1: Source Dataset Description Figure/ Company |  Company1 |  Company2
---|---|---
#Records |  27779 |  38162
#Purchase Orders (OrderID) |  6898 |  15455
#Items (ItemID) |  25961 |  36300
#Group Categories |  3 |  12
#Material Categories |  23 |  26
#Unique Item descriptions |  17198 |  12018
#Vendors |  988 |  1126
#Requesters |  122 |  63
#Buyers |  11 |  51
#Approvers |  76 |  101
Min/Max/Mean #OrderID / Requester |  1/56/870 |  1/245/1492
Min/Max/Mean #ItemID / Requester |  1/212/11884 |  1/576/4882
Min/Max/Mean #VendorCode / Requester |  1/14/102 |  1/31/137
Min/Max/Mean #OrderID / Buyer |  1/627/2141 |  1/303/1492
Min/Max/Mean #ItemID / Buyer |  1/2360/11907 |  1/712/4883
Min/Max/Mean #VendorCode / Buyer |  1/119/382 |  1/38/158
Min/Max/Mean #OrderID / Approver |  1/95/819 |  1/153/835
Min/Max/Mean #ItemID / Approver |  1/357/11574 |  1/359/3531
Min/Max/Mean #VendorCode per Approver |  1/25/107 |  1/28/174
Total Purchase Amount |  98,9MEUR |  243,9MEUR

The preparation of the data (cleaning and transformation) included several
tasks aimed to make it ready to the modelling phase.

The statistical profiling performed to understand the nature of the data
disclosed the high number of categorical features, which represents a
limitation to applying well-known and robust clustering algorithms. To solve
this, different Data Coding techniques were used to transform categorical
information into numeric.The categorical data encoding used was the Target
Encoding technique [21], consisting in substituting each group in a
categorical feature with the average response in another target and meaningful
variable of the dataset (the order amount in our case). This represents a
powerful solution because it avoids generating a high number of features, as
is the case of One-Shot Encoding, keeping the dimensionality of the dataset as
the original one.The idea behind was to analyse the behaviour of these four
strategies of target encoding with the algorithms and compare their results.

Other aspects revealed in these pre-analyses were: i) numerical attributes do
not follow a normal distribution, ii) the amount of a purchase order is the
most relevant numeric information in the dataset, so it would be used to
characterise categorical attributes, iii) the few rows removed due to the
presence of missing values, iv) no redundant information was found through
correlation measures so feature selection to reduce the problem dimensionality
was not possible.

Then, for every dataset it has been applied the different anomaly detection
methods, testing the several options for data grouping (median, avergage,
count…) in order to perform feature engineeerin with data, both in univariate
and multivariate approach.. This process will point out the more disciminant
way and method to detect the anomalies

##  5 Results

###  5.1 Univariate Anomaly Detection

A first approximation of this work has been using a univariate approach where
four techniques were selected and applied individually upon the 17 columns:

  * •

Numeric outliers: those transactions whose column value is outside of the IQR
rank (Lower bond: (Q1-1.5IQR) – Upper bond: (Q3+1.5IQR)) has been flagged as
an outlier.

  * •

z-Score: those transactions whose z-Score is greater than 2,5 has been flagged
as an outlier.

  * •

DBSCAN: those transactions that do not have at least three neighbours within a
distance of the unity (Core Points) nor are within a distance of the unity
from a Core Point (but have less than 3 neighbours), have been flagged as an
outlier.

  * •

Isolation Forest: (KNIME node executing the implementation in the
sklearn.ensemble Python package) [22].

IQR and IForest techniques performed excellently, but they were discarded as
they disclosed a high number of outliers (the adopted criteria were to provide
a manageable set of transactions to be reviewed by the company specialists).
Furthermore, z-Score and DBSCAN are based on different approaches (distance to
mean value and density, respectively) and, according to the results, can
provide a complementary point of view. As a result, the transactions (rows)
detected as outlier by z-Score or DBSCAN and for, at least, one attribute,
have been flagged as univariate outlier. Applying this rule, the resulting
dataset for each encoding has a different size (number of transactions
-rows-), as shown in Table 2.

Table 2: Number of transactions per dataset flagged by z-score or DBSCAN as an univariate outlier vs normal transactions. Company | Encoding | #UnivariateOutliers | #NormalElements | #TotalDataset
---|---|---|---|---
Company1 | Count | 2286 | 25435 | 27721
Company1 | Mean | 4546 | 23175 | 27721
Company1 | Median | 2870 | 24851 | 27721
Company1 | Mode | 4721 | 23000 | 27721
Company2 | Count | 2889 | 35049 | 37938
Company2 | Mean | 3656 | 34282 | 37938
Company2 | Median | 3431 | 34507 | 37938
Company2 | Mode | 3741 | 34197 | 37938

###  5.2 Multivariate Anomaly Detection

An iterative methodology has been followed, starting with the application to
the datasets a well-known and proven algorithm and, in successive steps, new
methods and algorithms have been added, enriching the final result.
Furthermore, a clustering quality optimization and a preliminary model
interpretability testing have been carried out.

####  5.2.1 k-Means – Optimised - With univariate outliers

After performing the pre-processing data tasks, the k-Means model has been
applied to the dataset. The univariate outliers have been kept in the dataset
as our goal was to detect anomalies, and the univariate flagged transactions
could represent a specific cluster. The clustering algorithm uses the
Euclidean distance on the selected attributes. As the data is not normalized
by the node, the "Normalizer" node is used as a pre-processing step to compute
a Gaussian normalization (Z-score).

The Elbow Curve and the Silhouette analysis have been the methods employed to
choose an appropriate value for the number of clusters. The elbow method runs
k-means clustering on the dataset for a range of values of k (2-25) and, for
each one, the average distances to the centroid across all data points is
computed. The cut-off point where the diminishing SSE is no longer worth the
additional number of clusters, indicates the optimal number of clusters. A
comprehensive analysis has been performed to find the optimal number of
clusters and the encoding strategy (see Figure 1) for results in Company1.

Figure 1: Comparison of the Elbow Curve by encoding strategy (frequency, mean,
median and mode). Company 1

Because of the visual criteria’s subjectivity inherent to the Elbow method,
Silhouette coefficient has been used the analytical method to contrast the
result. The higher the coefficient, the more compact are the clusters and
easily differentiable from the rest. The interpretation of the silhouette
coefficient adopted on this work has been the following [20]:

  * •

Sc > 0,7: strong cluster structure

  * •

Sc > 0,5: reasonable cluster structure

Figure 2 shows the silhouette coefficient computed for each number of clusters
tested (parameter k) and the encoding strategy for Company1. In general, the
Median, Mean and Mode strategies provides better results as their overall
coefficient is higher than the Count strategy with any k value. This means
that the clusters are more compact and easily differentiable from the rest.

Figure 2: Comparison of Silhouette coefficient by encoding strategy
(frequency, mean, median and mode). Company 1

According to the adopted criteria, at this stage of the work, the k-Means
model chosen for Company1 has been the dataset encoded by Mode strategy and 8
clusters (SSE=208.459 Silhouette=0,41) and, for Company2, the dataset encoded
by Median strategy and 5 clusters (SSE=267.784 Silhouette=0,34). In both
cases, the coefficients obtained are below the threshold of a reasonable
cluster structure.

Due to its high-demanding computational cost, as approximation to efficiently
compute the coefficient in large datasets, a stratified random sampling (10%
of the total population) has been performed before its calculation [23]. This
sampling method provides an unbiased and accurate representation of the a
priori distribution of a feature in the underlying population. Given the poor
values obtained (Silhouette was below 0,5 for every number of clusters between
2 and 25 and for any of the encoding datasets), a coefficient’s full
computation for the Company1 dataset (including univariate outliers) has shown
quite similar results so the sampling method used for the Silhouette
computation shows as a good approximation.

####  5.2.2 k-Means – Optimised - Without univariate outliers

As seen in the previous step, the indicators of the cluster validity were not
suitable according to commonly adopted criterion (reasonable cluster
structure: Sc > 0.5). For this reason, the univariate outliers have been
segregated from the datasets and the process has been repeated (previously, a
re-encoding of the categorical features has been computed). Figures 3 and 4
show the Elbow curves and the Overall Silhouette coefficient in relation to
the number of clusters and for each encoding strategy for the Company1 dataset
with the univariate outliers segregated. The Elbow curve for Count remains
approximately at the same SSE levels while Mean, Median and Mode curves have
decreased their sum of squared distances to centroids closer to Count values.
As expected, filtering out the univariate outliers has a direct impact on the
encoding strategies depending on the central trending of the features/columns.

Figure 3: Company1 - Comparison of the Elbow Curves by encoding strategy
(frequency, mean, median and mode) with univariate outliers segregated

Figure 4: Company1 - Comparison of the Overall Silhouette coefficient – 10%
sampled by encoding strategy (frequency, mean, median and mode) with
univariate outliers segregated

The Silhouette coefficient almost reaches an acceptable level in k=2-3 for the
Mean and Median datasets although, in comparison with the results with no
segregation of outliers (Figure 4), all the datasets decrease dramatically
from this value of k. In this case, a low number of clusters is not adequate
for our problem, given that they will contain a high number of transactions
(e.g. Median dataset with k=2 clusters, the smaller cluster with 2217
transactions) whose review by an expert is unfeasible. Similar results have
been obtained for Company2.

Filtering out the univariate outliers from the datasets should improve the
quality indicators as long-distant transactions are no longer considered for
the computation. However, in terms of clustering quality indicators, a
significant improvement has not been achieved and the hyperparameters selected
to model each company are those adopted without segregation of the univariate
outliers).

####  5.2.3 Isolation Forest

The KNIME nodes used are implemented using H2O [24]:

  * •

Node H2O Isolation Forest Learner. The model is learned unsupervised from the
input. The main training hyper-parameters used are the number of trees to
build (nTree=100) and the subsampling size to train each tree
(sample_size=256); both are the defaults of the original paper [19]. A static
random seed has been used.

  * •

Node H2O Isolation Forest Predictor. The node applies the model resulting from
the learner to an input dataset to predict anomalies or outliers. The output
of the node will consist of:

    * –

Prediction expressed as a normalized anomaly score. The higher the score, the
more likely it is an anomaly.
$Prediction=\frac{MaxPathLenth-MeanPathLength}{MaxPathLength-MinPathLength}$

    * –

Mean length of the predicted decision tree paths of each observation. The
shorter, the more likely it is an anomaly.
$MeanLength=\frac{PathLength}{nTree}$

Given that the dataset is not labelled and, thus, unsupervised algorithms are
required, the learning and predictor stage of the algorithm has been applied
upon each one of the whole companies’ datasets.

Due to the unknown proportion of outliers in the unlabelled datasets used,
iForest cannot accurately set a threshold to determine whether a certain
transaction is put into the anomaly candidates set. In our case, the anomaly
threshold for the IForest algorithm is established as the minimum score of the
1% of the transactions ordered descending by the prediction score (with a
maximum of 500 transactions). This approach is based on the limited
availability of the business specialists, and it can be revisited according to
the actual effort they can dedicate. This threshold is computed as the 99th
quantile of the Prediction provided by IForest and flagging true/false
according to each transaction .

###  5.3 Anomalies Prioritisation - Ensemble Algorithm

As a preliminary step to establish a global anomalies prioritisation, the
dataset with an optimal k-Means model has been chosen and, upon these, the
adopted criteria to prioritise the anomalies have been the following:

  1. i

Low populated clusters. Initially, the clusters populated less than 1% of the
total dataset have been defined as the analytical threshold. The transactions
belonging to those clusters have been flagged as “k-Means anomaly”.

  2. ii

Single transactions which differ from others in the same cluster. In other
words, outliers within the clusters. The transactions with an individual
Silhouette coefficient negative have been flagged as “Silhouette anomaly”.

  3. iii

Isolation Forest prediction. According to the effort-based approach, already
explained in section 5.2.3, those transactions beyond the 99th quantile of the
algorithm prediction have been flagged as “iForest anomaly”. An absolute
maximum of 500 transactions is established.

  4. iv

Univariate anomalies. In the case of not having segregated the univariate
outliers from the dataset, those transactions with at least one univariate
anomaly -according to z-Score or DBSCAN - have been flagged as “univariate
outlier anomaly”.

Criteria i) and ii) correspond to the 2nd and 3rd categories of clustering-
based techniques for anomaly detection proposed by [15]. Some studies have
adopted this categorisation to establish the key assumptions to identify
anomalies [25, 26].

In a second step, a re-prioritisation can be carried out in case of an
excessive number of anomalies -with respect to the specialist’s availability
to review the results-. For example, if the number of “k-Means anomalies” must
be reduced, they can be ordered by the individual Silhouette coefficient (the
lower, the more critical) or by the iForest prediction (the higher, the more
critical).

A simple global prioritisation, a kind of independent ensemble, has been
implemented flagging the transactions with the number of anomaly groups they
belong to and ordering them (descending). The resulting prioritised groups are
shown in Table 3 and Table 4. For a better understanding of the anomaly
location, the distribution of the prioritised groups within the k-Means
clusters is provided.

Table 3: Company1 - Mode dataset. k=8. Global prioritisation groups by type of anomaly Anomaly Priority |  k-Means Anomaly |  Silhouette Anomaly |  iForest Anomaly |  Univariate Anomaly |  cluster 0 |  cluster 1 |  cluster 2 |  cluster 3 |  cluster 4 |  cluster 5 |  cluster 6 |  cluster 7 |  Total
---|---|---|---|---|---|---|---|---|---|---|---|---|---
2 | No | No | Yes | Yes | 1 | 19 | 3 | 4 | 14 | 0 | 0 | 3 | 44
2 | No | Yes | No | Yes | 34 | 0 | 3 | 11 | 0 | 0 | 3 | 100 | 151
2 | No | Yes | Yes | No | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 7 | 7
1 | No | No | No | Yes | 1411 | 1484 | 305 | 335 | 290 | 226 | 202 | 273 | 4526
1 | No | No | Yes | No | 0 | 0 | 0 | 14 | 158 | 6 | 9 | 45 | 232
1 | No | Yes | No | No | 0 | 0 | 0 | 89 | 0 | 0 | 9 | 768 | 866
0 | No | No | No | No | 0 | 0 | 0 | 5150 | 5263 | 4352 | 4183 | 2947 | 21895
Table 4: Company2 - Median dataset. k=5. Global prioritisation groups by type of anomaly Anomaly Priority |  k-Means Anomaly |  Silhouette Anomaly |  iForest Anomaly |  Univariate Anomaly |  cluster 0 |  cluster 1 |  cluster 2 |  cluster 3 |  cluster 4 |  Total
---|---|---|---|---|---|---|---|---|---|---
3 | Yes | Yes | No | Yes | 60 | 0 | 0 | 0 | 0 | 60
2 | No | No | Yes | Yes | 0 | 10 | 0 | 21 | 7 | 38
2 | No | Yes | No | Yes | 0 | 4 | 0 | 0 | 71 | 75
2 | Yes | No | No | Yes | 28 | 0 | 0 | 0 | 0 | 28
1 | No | No | No | Yes | 0 | 919 | 633 | 796 | 882 | 3230
1 | No | No | Yes | No | 0 | 147 | 10 | 173 | 11 | 341
1 | No | Yes | No | No | 0 | 0 | 0 | 0 | 355 | 355
0 | No | No | No | No | 0 | 8990 | 8437 | 8456 | 7928 | 33811

Other approaches are possible to set analytical rules and merge the individual
algorithms’ results. [27] uses these kind of approaches to combine anomalies
detected in audit logs for application systems by LOF and DBSCAN algorithms.

###  5.4 Model Interpretability

Interpretable Machine Learning (IML) methods can be used to discover
knowledge, debug or justify the model and its predictions, and control and
improve the model [28]. Methods that study the sensitivity of an ML model are
mostly model-agnostic and work by manipulating input data and analysing the
respective model predictions. These IML methods often treat the ML model as a
closed system that receives feature values as an input and produces a
prediction as output.

Some of the available explicability tools (Shapley values, SHAP and Lime,
[28]) has been tested to provide an overview of how they can help the business
specialists in the task of understanding why an outlier has been flagged [28].
A further investigation needs to be accomplished to define the better and most
reliable explicability method to support the specialists in their analysis.
Figure 5 shows a test using the SHAP package in Python. In the sklearn
implementation of iForest, the lower the score is, the more abnormal [22]. For
the example shown in the figure, the algorithm has been applied to the
Company2 dataset, disclosing that the requester, the buyer and the department
(organization code) have the major contributions to the anomaly deviation with
respect to the dataset average prediction.

Figure 5: Company2. Example of SHAP values. (left) Force plot for individual
explicability. (right) Bar plot for feature explicability

##  6 Conclussion

In this paper, Machine Learning based anomaly detection has been performed on
real purchase datasets, in order to provide a disciminant tool for fraud
detection in auditory processes. A combination of several anomaly detection
techniques (ADT) is proposed: univariate detection (through z-score and
DBSCAN), low populated clusters and the negative Silhouette coefficients
defined by k-Means, and the Isolation Forest.

The most suitable models generated by the k-Means algorithm presented a low
clustering quality according to the overall Silhouette index (below 0,5).
Complementary analysis, as using datasets without univariate outliers or
performing a full computation of the Silhouette coefficient, has not improved
the quality measures.

Concerning the adopted assumptions to identify anomalies, low populated
clusters were assumed as collective outliers, flagging transactions belonging
to clusters representing less than 1% of the total dataset. To provide further
information and to better assess the results, the rows (transactions) with
negative Silhouette index were flagged as potential contextual anomalies. To
complement the cases with a different approach and enrich the model with a
wide variety of the detectable anomalies’ nature, the transactions with the
highest prediction computed through Isolation Forest algorithm were selected
(those higher than 99th quantile). Finally, the preliminary univariate outlier
identification performed during the EDA phase was retained in the output,
flagging transaction with, at least, one univariate outlier on its features.

As showed in this work, the use of the KNIME has been revealed as an easy-to-
use and well-documented tool. It enables the use of a wide range of Machine
Learning techniques without much effort, and makes possible to automatize
processes, being a key aspect in iterative investigations. Python integration
has allowed to verify some the results and improve the output’s look and feel.

The work has revealed the extensive possibilities that Machine Learning, and
particularly clustering, offers to unsupervised learning problems as the
procurement process. Among others, some axes of future investigation are
testing of other encoding options for categorical features, performing a time
analysis for the purpose of detecting stational behaviours (splitting purchase
transactions in a short period of time is among the major internal control
concerns in a company), extending the use of clustering algorithms to others
(such as DBSCAN, Fuzzy c-Means - overlapping clustering -, hierarchical
clustering or Self-Organised Maps -SOM-). This certainly will extend the
nature of the anomalies addressed by the proposed methodology.

## Acknowledgment

This work is partially supported by MCIN/AEI/ 10.13039/501100011033 by “ERDF A
way of making Europe”, grant PID2021-127946OB-I00. The author, S. Pérez-Díaz,
is partially supported by Ministerio de Ciencia, Innovación y Universidades -
Agencia Estatal de Investigación/PID2020-113192GB-I00 (Mathematical
Visualization: Foundations, Algorithms and Applications).

## References

  * [1] IIA,“Definition of Internal Auditing,” [Online]. Available: https://www.theiia.org/en/standards/what-are-the-standards/definition-of-internal-audit/. Accessed on: Dec. 27 2023.
  * [2] IIA, “International Standards for fhe Professional Practice of Internal Auditing (Standards),” [Online]. Available: https://www.theiia.org/globalassets/site/standards/mandatory-guidance/ippf/2017/ippf-standards-2017-english.pdf Accessed on: Dec. 27 2023.
  * [3] P. Byrnes, A. Al-Awadhi, B. Gullkvist, H. Brown-Liburd, R. Teeter, J. Warren and M. Vasarhelyi, “Evolution of Auditing: From the Traditional Approach to the Future Audit: Theory and Application,” 2018.
  * [4] V. Chiu, Q. Liu and M. A. Vasarhelyi, “The Development and Intellectual Structure of Continuous Auditing Research,” in _Continuous Auditing_ , Emerald Publishing Limited, 2018, p. 53–85.
  * [5] P. Hajek and R. Henriques, “Mining corporate annual reports for intelligent detection of financial statement fraud – A comparative study of machine learning methods,” _Knowledge-Based Systems_ , vol. 128, p. 139–152, July 2017.
  * [6] Q. Liu, “The application of exploratory data analysis in auditing,” M.S. thesis, Graduate School-NewarkRutgers, New Jersey, 2014.
  * [7] J. Kokina and T. Davenport, “The Emergence of Artificial Intelligence: How Automation is Changing Auditing,” _Journal of Emerging Technologies in Accounting,_ vol. 14, April 2017.
  * [8] J. Nonnenmacher and J. Gómez, “Unsupervised anomaly detection for internal auditing: Literature review and research agenda,” _The International Journal of Digital Accounting Research,_ pp. 1-22, January 2021.
  * [9] A. Sharma and P. K. Panigrahi, “A Review of Financial Accounting Fraud Detection based on Data Mining Techniques,”, /emphInternational Journal of Computer Applications, September 2013.
  * [10] M. Jans, N. Lybaert and K. Vanhoof, “Data Mining for Fraud Detection: Toward an Improvement on Internal Control Systems?,” February 2008.
  * [11] M. Jans, N. Lybaert and K. Vanhoof, “Internal fraud risk reduction: Results of a data mining case study,” _International Journal of Accounting Information Systems,_ vol. 11, p. 17–41, March 2010.
  * [12] PwC, “Global Economic Crime Survey 2022. Technical report,” 2022.
  * [13] V. Kotu and B. Deshpande, _Data Science: Concepts and Practice,_ Elsevier Science & Technology Books, 2019, p. 568.
  * [14] E. Soria, J. D. Martín, A. J. Serrano and D. Aguado, _Análisis de datos experimentales,_ Ed. UPV, Valencia, Spain, 2007.
  * [15] V. Chandola, A. Banerjee and V. Kumar, “Anomaly Detection: A Survey,” _ACM Comput. Surv.,_ vol. 41, July 2009.
  * [16] V. Kotu and B. Deshpande, “Anomaly Detection,” in _Data Science,_ Elsevier, 2019, p. 447–465.
  * [17] M. D. J. Bora and D. A. K. Gupta, “Effect of Different Distance Measures on the Performance of K-Means Algorithm: An Experimental Study in Matlab,” _International Journal of Computer Science and Information Technologies,(IJCSIT),_ Vol. 5 (2), 2501-2506, May 2014.
  * [18] P. J. Rousseeuw, “Silhouettes: A graphical aid to the interpretation and validation of cluster analysis,” _of Computational and Applied Mathematics,_ vol. 20, p. 53–65, November 1987.
  * [19] F. T. Liu, K. M. Ting and Z.-H. Zhou, “Isolation-Based Anomaly Detection,” /emphACM Trans. Knowl. Discov. Data, vol. 6, March 2012.
  * [20] KNIME, “Introduction to Machine Learning Algorithms,” 1 January 2021. [Online]. Available: https://www.knime.com/sites/default/files/2021-01/2021-01-slides-l4- ml.pdf. Accessed on: Dec. 27 2023.
  * [21] E. Zuccarelli, “Handling Categorical Data, The Right Way,” 1 August 2020. [Online]. Available: https://towardsdatascience.com/handling-categorical-data-the-right-way-9d1279956fc6. Accessed on: Dec. 27 2023.
  * [22] Scikit-Learn, “Isolation Forest Algorithm,”” [Online]. Available: https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.IsolationForest.html. Accessed on: Dec. 27 2023.
  * [23] V. Palacios, M. Widmann and R. Cadili, “The pros and cons of statistical sampling methods plus how to find the right strategy,” 8 November 2021. [Online]. Available: https://www.knime.com/blog/statistical-sampling. Accessed on: Dec. 27 2023.
  * [24] H2O, “Isolation Forest - H2O implementation,” [Online]. Available: https://docs.h2o.ai/h2o/latest-stable/h2o-docs/data-science/if.html. Accessed on: Dec. 27 2023.
  * [25] S. Thiprungsri and M. Vasarhelyi, “Cluster Analysis for Anomaly Detection in Accounting Data: An Audit Approach,” _The International Journal of Digital Accounting Research,_ vol. 11, 2011.
  * [26] M. Ahmed, A. N. Mahmood and M. R. Islam, “A survey of anomaly detection techniques in financial domain,” in _Future Generation Computer Systems,_ vol. 55, p. 278–288, February 2016.
  * [27] H. D. Kuna, R. Garcia-Martinez and F. R. Villatoro, “Outlier detection in audit logs for application systems,” in _Information Systems,_ vol. 44, p. 22–33, August 2014.
  * [28] C. Molnar, “Interpretable Machine Learning,” 2022. [Online]. Available: https://christophm.github.io/interpretable-ml-book. Accessed on: Jan. 1 2024.


# Superfluid Fraction and Leggett Bound in a Density Modulated Strongly
Interacting Fermi Gas at Zero Temperature

G. Orso<EMAIL_ADDRESS>Université Paris Cité, Laboratoire
Matéeriaux et Phéenomènes Quantiques (MPQ), CNRS, F-75013, Paris, France  S.
Stringari<EMAIL_ADDRESS>Pitaevskii BEC Center, CNR-INO and
Dipartimento di Fisica, Università di Trento, Via Sommarive 14, 38123 Povo,
Trento, Italy

###### Abstract

We calculate the superfluid fraction of an interacting Fermi gas, in the
presence of a one-dimensional periodic potential of strength $V_{0}$ and wave-
vector $q$. Special focus is given to the unitary Fermi gas, characterized by
the divergent behavior of the s-wave scattering length. Comparison with the
Leggett’s upper bound $(\langle n_{1D}\rangle<1/n_{1D}>)^{-1}$, with $n_{1D}$
the 1D column density, explicitly shows that, differently from the case of a
dilute interacting Bose gas, the bound significantly overestimates the value
of the superfluid fraction, except in the phonon regime of small $q$. Sum rule
arguments show that the combined knowledge of the Leggett bound and of the
actual value of the superfluid fraction allows for the determination of
curvature effects providing the deviation of the dispersion of the Anderson-
Bogoliubov mode from the linear phonon dependence. The comparison with the
predictions of the weakly interacting BCS Fermi gas points out the crucial
role of two-body interactions. The implications of our predictions on the
anisotropic behavior of the sound velocity are also discussed.

##  I I. Introduction

In 1970 Leggett [1] derived an upper bound to the superfluid fraction
$f_{S}=n_{S}/\bar{n}$, given by the ratio between the superfluid density
relative to the motion of the fluid along the $x$-direction and the mean
density, in terms of the spatial average of the inverse of the 1D column
density $n_{1D}(x)=\int n(\mathbf{r})dydz$ [2]:

$f_{S}\leq(\langle{n}_{1D}(x)\rangle\langle\frac{1}{n_{1D}(x)}\rangle)^{-1}\equiv
f^{L}_{S}.\;$ (1)

Result (1) is particularly relevant at zero temperature, where the superfluid
fraction is affected by the density modulations induced by the breaking of
translational invariance, and holds independent of whether the breaking occurs
spontaneously or it is induced by an external force. At finite temperature the
superfluid fraction is instead affected by thermal effects even for uniform
systems, with important consequences on the propagation of second sound [3, 4,
5]. Under the assumption that the many-body wave function is described by a
single phase equal for all the particles, an assumption applicable to a dilute
Bose gas at zero temperature described by the Gross-Pitaevskii equation,
Leggett derived a further lower bound [6] to the superfluid fraction. The
Leggett upper and lower bounds coincide if the 3D density depends only on the
$x$ variable, so that in this case the inequality (1) reduces to an equality.
In a recent experiment [7] the Leggett’s integral entering Eq.(1) was actually
measured in an ultracold Bose gas confined in a 2D box, in the presence of a
periodically modulated external 1D potential, causing important
disomogeneities in the density profile, and the bound (1) was shown to agree
with good accuracy with the result for the superfluid fraction measured
through the velocity of sound, in agreement with the predictions of Gross-
Pitaevskii theory.

The question of the accuracy of Leggett bound in other interacting many-body
systems is basically unexplored and the aim of our work is to address the
question in the case of an interacting Fermi gas at zero temperature, with
special focus on the so-called unitary limit, where the $s$-wave scattering
length is infinite and the system behaves as a strongly interacting
superfluid. We will consider a gas confined in a box in the presence of a
periodic potential of the form

$V_{ext}=V_{0}\cos(qx)$ (2)

which induces density modulation of period $d=2\pi/q$ in the ground state
density profile. We will attack the problem following two different lines. In
Sect. II we develop a sum rule approach to explore the behavior of the
superfluid fraction in the presence of a weak perturbation (small $V_{0}$) and
calculate the Leggett’s integral as well as the exact value of the superfluid
fraction. The sum rule results provide the proof that only if the excitation
spectrum of the unperturbed uniform superfluid consists of a single excitation
exhausting the relevant sum rules of the dynamic structure factor, will the
inequality (1) reduces to an identity. This is the case of the dilute BEC gas
for all values of the wave vector $q$, but remains an inequality in the most
general case. The results of the sum rule approach are also used to infer on
the curvature of the dispersion of the Anderson-Bogoliubov collective mode of
a Fermi superfluid at long wavelengths. In Sect. III we first provide a
quantitative implementation of the sum rule approach developed in Sect II by
applying the dynamical Bardeen–Cooper–Schrieffer (BCS) theory [8]. We then
present our numerical results for the superfluid fraction of the unitary Fermi
gas by solving the Bogoliubov-de Gennes (BdG) equations in the presence of the
external periodic potential (2) in order to provide a quantitative comparison
between the Leggett bound and the actual value of the superfluid fraction,
beyond the weak perturbation limit of small $V_{0}$, in a wider range of
parameters of possible relevance for future experiments. In Sect IV we discuss
how our predictions for the superfluid density could be tested by measuring
the effects of anisotropy in the velocity of sound.

##  II II. Sum rule approach

At zero temperature, the dynamic structure factor $S(\omega,q)$ of a many-body
system is related to the imaginary part of the density response function
$\chi(\omega,q)$ by

$S(\omega,\mathbf{q})=-\frac{1}{\pi}\textrm{Im}\chi(\omega+i0^{+},\mathbf{q}),$
(3)

so that its spectral representation is

$S(\omega,\mathbf{q})=\frac{1}{V}\sum_{n}|\langle
0|\rho_{\mathbf{q}}|n\rangle|^{2}\delta(\omega+E_{0}-E_{n}),$ (4)

where $V$ is the volume of the box,
$\rho_{\mathbf{q}}=\sum_{k}\exp(i\mathbf{q}\cdot\mathbf{x}_{k})$ is the
Fourier transform of the density operator, $|n\rangle$ is a complete set of
eigenstates of the many-body Hamiltonian and $(E_{n}-E_{0})$ are the
corresponding excitation energies. Here and in the following we use the
convention $\hbar=1$ and consider configurations where the dynamic structure
factor $S(\omega,\mathbf{q})$ depends on the wave-vector only through its
modulus $q=|\mathbf{q}|$. From the dynamic structure factor one can extract
the energy weighted moments

$m_{j}(q)=\frac{1}{\bar{n}}\int_{0}^{\infty}d\omega\omega^{j}S(\omega,q)$ (5)

where $j$ is an integer and $\bar{n}=N/V$ is the average fermion density. The
energy weighted moment $m_{1}(q)$ is model independent and equal to
$q^{2}/(2m)$, as follows from the number conservation law (f-sum rule), with
$m$ the particle mass. All the other moments are generally model dependent and
do carry crucial information on the many-body properties of the system, as we
shall see below.

###  II.1 A. Superfluid fraction $f_{S}$ vs Leggett bound $f_{S}^{L}$

If the periodic potential (2) is weak (small $V_{0}$), the density changes
induced by the perturbation can be calculated via linear response theory and
take the form

$\delta n(x)=V_{0}\chi(0,q)\cos(qx)\;,\\\ $ (6)

where $\chi(\omega=0,q)$ is the static response function, related to the
inverse energy weighted moment of the dynamic structure factor $S(\omega,q)$
by [9]:

$-\frac{1}{\bar{n}}\chi(0,q)=2m_{-1}(q)$ (7)

By expanding the Leggett’s integral up to terms quadratic in $V_{0}$, one
finds the following result for the Leggett’s upper bound:

$f^{L}_{S}=1-2V^{2}_{0}(m_{-1}(q))^{2}$ (8)

in terms of $V_{0}$ and $m_{-1}(q)$. Since Eq.(8) is based on linear response
theory, the moment $m_{-1}(q)$ in Eq.(8) has to be evaluated using the
unperturbed ($V_{0}=0$) uniform system.

In the following we will consider many body systems, characterized by the
absence of transverse gapless excitations. This includes, among others, single
component Bose and Fermi superfluids in the presence of an external periodic
perturbation of the form (2). In order to calculate the superfluid fraction
along the $x$-direction it is convenient to calculate the expectation value
$\langle P_{x}\rangle$ of the longitudinal momentum
$P_{x}=\sum_{k=1}^{N}p_{k,x}$ acquired by the system in the presence of the
longitudinal perturbation $-vP_{x}$ and use the definition

$f_{S}=1-\lim_{v\to 0}\frac{\langle P_{x}\rangle}{Nmv}$ (9)

By applying perturbation theory the superfluid fraction can then be written in
the form

$f_{S}=1-\frac{2}{Nm}\sum_{n\neq 0}\frac{|\langle
0|P_{x}|n\rangle|^{2}}{E_{n}-E_{0}}.$ (10)

If the Hamiltonian is translationally invariant, i.e. if $[H,P_{x}]=0$, then
$f_{S}=1$. If instead the Hamitonian contains the additional term (2) breaking
translational invariance, then the system will acquire a normal component.
Using the commutation relation

$[H,P_{x}]=-\frac{q}{2}V_{0}(\rho_{q}-\rho_{-q})\;,$ (11)

one can rewrite Eq.(10) in the form

$f_{S}=1-\frac{1}{Nm}V_{0}^{2}q^{2}\sum_{n\neq 0}\frac{|\langle
0|\rho_{q}|n\rangle|^{2}}{(E_{n}-E_{0})^{3}}=1-2V_{0}^{2}m_{1}(q)m_{-3}(q),$
(12)

showing that the superfluid fraction is controlled by the inverse cubic energy
weighted moment $m_{-3}(q)$. While Eq.(8) is valid only for a weak periodic
potential, Eq. (12) holds for arbitrary values of $V_{0}$. Of course for small
$V_{0}$ the moment $m_{-3}$ can be safely evaluated using the unperturbed
uniform system.

Eq.s (8) and (12) are the main results of this section. They show that, in the
considered limit of weak perturbations (small $V_{0}$), the Leggett’s bound
coincides with the exact value of the superfluid fraction only if the general
inequality

$\frac{m_{1}}{m_{-1}}\geq\frac{m_{-1}}{m_{-3}}\;,$ (13)

which follows from the positiveness of the domain of the dynamic structure
factor, reduces to an equality. The equality holds only in few notable cases:
i) for a dilute BEC gas, where the Bogoliubov excitation spectrum is
characterized by a single mode for all values of the wave vector $q$, with
frequency $\omega_{=}\sqrt{q^{2}\mu_{b}/m+q^{4}/(2m)^{2}}$, where $\mu_{b}$ is
the mean field chemical potential of the BEC gas. ii) In general for
superfluids in the phonon limit $q\to 0$, as a consequence of the fact that in
this limit the phonon mode exhausts the relevant sum rules entering the
inequality (13). This includes, in particular the case of the interacting
superfluid Fermi gas considered in this paper but also the case of superfluid
${}^{4}He$. In general the inequality (13) implies that $f_{S}^{L}\geq f_{S}$,
in agreement with Leggett’s prescription.

It is useful to stress that the existence of a finite $m_{-3}$ for a Fermi gas
requires the particle-hole excitations spectrum to be gapped (as it is indeed
the case in the superfluid state), otherwise $m_{-3}$ diverges due to the
presence of gapless particle-hole excitations at the Fermi surface [9].

###  II.2 B. Curvature effect in the excitation spectrum

Interestingly, the moments $m_{-1}$ and $m_{-3}$, controlling the behavior of
the Leggett bound $f_{S}^{L}$ and of the actual value $f_{S}$ of the
superfluid fraction in the presence of the periodic potential (2), can also be
used to characterize the dispersion relation of the Anderson-Bogoliubov mode
of a uniform system at small $q$

$\omega^{2}=c^{2}q^{2}(1+\gamma q^{2}),$ (14)

where $c$ is the sound velocity, related to the inverse compressibility
$\kappa^{-1}=n(\partial\mu/\partial n)$, through the thermodynamic relation
$c^{2}=\kappa^{-1}/m$, and $\gamma$ measures the curvature effect beyond the
linear dispersion. The coefficient $\gamma$ is positive in the case of a
dilute BEC gas $(\gamma=1/(4\mu_{b})$), while it is expected to become
negative in the BCS region of the BEC-BCS crossover. Its change of sign takes
place around unitarity and has been the object of several theoretical
calculations [10, 11, 12, 13] and of a recent measurement in a superfluid
atomic Fermi gas using Bragg spectroscopy [14].

To understand the connection between the results for the superfluid fraction
in the presence of the perturbation $V_{0}\cos(qx)$ and the curvature effect,
it is useful to explore the behavior of the contributions to the various sum
rules originating from the gapless branch and from the gapped excitations in
the small $q$ regime (see Table 1). The $q^{4}$ dependence of the strenght
contribution of the gapped states follows from general arguments based on
Galilean invariance (see for example [9]) and implies that the $m_{-3}$ and
$m_{-1}$ moments of the dynamic structure factor are exhausted by the gapless
mode not only in the lowest $q$ phonon regime, but also including the first
corrections in $q^{2}$.

Table 1: Low momentum behavior of the leading contributions to a given physical quantity (left column), coming from gapless collective modes (central column) and from gapped single-particle excitations (right column). Notice that the sum of the $q^{4}$ contributions of the gapless and gapped excitations to the energy weighted moment $m_{1}$ cancel out in order to ensure the fulfillment of the $f$-sum rule, $m_{1}(q)=q^{2}/(2m)$.  | Gapless | Gapped
---|---|---
$|\langle 0|\rho_{q}|n\rangle|^{2}$ | $q+q^{3}$ | $q^{4}$
$\omega(q)$ | $q+q^{3}$ | const
$m_{1}(q)$ | $q^{2}+q^{4}$ | $q^{4}$
$m_{-1}(q)$ | const + $q^{2}$ | $q^{4}$
$m_{-3}(q)$ |  $1/q^{2}$ \+ const | $q^{4}$

A non trivial consequence is that, at low $q$, the average square excitation
energy

$\omega^{2}=m_{-1}(q)/m_{-3}(q)$ (15)

exactly coincides with the expansion (14) of the square frequency, including
the curvature correction in $\gamma$. Viceversa, the most popular ratio
$m_{1}(q)/m_{-1}(q)$ correctly reproduces only the leading phonon term
$c^{2}q^{2}$. Using the small $q$ expansions
$m_{-1}(q)=(\kappa/2)(1-\beta_{-1}q^{2})$ and
$q^{2}m_{-3}(q)=(m\kappa^{2}/2)(1-\beta_{-3}q^{2})$ in Eq.(15) and comparing
the obtained expression with Eq.(14) one finds

$\gamma=\beta_{-3}-\beta_{-1}.$ (16)

By substituting the same expansions in Eq.s (8) and (12), we note that the
inequality $f_{S}\leq f_{S}^{L}$ yields $\beta_{-3}\leq 2\beta_{-1}$. When
$\gamma$ vanishes, corresponding to $\beta_{-1}=\beta_{-3}$, the $q^{2}$
corrections to the Leggett bound and to the exact superfluid fraction differ
exactly by a factor $2$. In the case of a superfluid Fermi gas this happens
very close to unitarity [14]. The situation is different for a dilute BEC gas,
where $\beta_{-3}=2\beta_{-1}$, since $f_{S}^{L}=f_{S}$ for all values of $q$.

It’s also interesting to compare the above results with the effective field
theory for the unitary Fermi gas developed in [10], which includes the leading
order corrections in an expansion in small gradients of the quantum
hydrodynamic Lagrangian density (see also the related discussion in Ref.[11]).
The corresponding predictions for the static response function and for the
dispersion of the collective mode at small $q$ are given by [10]

$\displaystyle\chi(0,q)$ $\displaystyle=$
$\displaystyle-n\kappa\left[1+A\left(c_{1}-\frac{9}{2}c_{2}\right)\frac{q^{2}}{k_{F}^{2}}\right],$
(17) $\displaystyle\omega^{2}$ $\displaystyle=$ $\displaystyle
c^{2}q^{2}\left[1-A\left(c_{1}+\frac{3}{2}c_{2}\right)\frac{q^{2}}{k_{F}^{2}}\right],$
(18)

where $c_{1}$ and $c_{2}$ are dimensionless parameters entering the Lagrangian
density, $k_{F}=(3\pi^{2}n)^{1/3}$ is the Fermi wave-vector,
$A=2\pi^{2}\sqrt{2(1+\beta)}$ and $\beta$ is the dimensionless Bertsch
parameter entering the equation of state of the homogeneous Fermi gas at
unitarity: $\mu=(1+\beta)k_{F}^{2}/(2m)$. The BdG equations yield
$\beta=-0.41$ (see Sec. III), while more accurate Monte Carlo calculations
predict $\beta=-0.58$ [15, 16], which is quite close to the value
$\beta=-0.63$ reported in the experiment [17].

By comparing Eqs. (17) and (18) with Eq.(7) and (15), respectively, we find
that the parameters of the effective field theory can be expressed as

$\displaystyle c_{1}$ $\displaystyle=$
$\displaystyle(2\beta_{-1}-3\beta_{-3})\frac{k_{F}^{2}}{4A},$ $\displaystyle
c_{2}$ $\displaystyle=$
$\displaystyle(2\beta_{-1}-\beta_{-3})\frac{k_{F}^{2}}{6A}.$ (19)

In Ref.[11] it was demonstrated that the coefficient $c_{2}$ must be positive.
From Eq.(II.2) we see that this is equivalent to the condition
$f_{S}<f_{S}^{L}$.

##  III III. Numerical results

###  III.1 A. Dynamical BCS theory and sum rules

In this section we evaluate the sum rule moments of a unitary Fermi gas with
the help of the dynamical BCS theory developed in Ref.[8] (see also an earlier
work by Minguzzi and collaborators [18]). In particular we will make use of
the general expression for the density response function

$\displaystyle\chi(\omega,q)$ $\displaystyle=$
$\displaystyle-\frac{1}{2\pi^{2}}\Big{(}I^{\prime\prime}(\omega,q)$ (20)
$\displaystyle-$
$\displaystyle\frac{F(\omega,q)}{I_{11}(\omega,q)I_{22}(\omega,q)-\omega^{2}I_{12}^{2}(\omega,q)}\Big{)},$

where $I^{\prime\prime},F$ and $I_{ij}$ are functions, whose definition in
terms of integrals over momenta can be found in Ref.[8]). Notice that, with
respect to the result quoted in Ref.[8] (see their Eq.(25)), the rhs of
Eq.(20) contains an extra factor of 2, because in our definition of the
density response function we account for the spin degeneracy [see Eq. (6)].
The first term in Eq.(20), proportional to $I^{\prime\prime}$, represents the
bare contribution from the linear response formalism applied to the BCS state,
while the second term accounts for the fluctuations of the order parameter
over its mean-field value, which are included via the random phase
approximation (RPA). In particular the vanishing of the denominator of the
second term,
$I_{11}(\omega,q)I_{22}(\omega,q)-\omega^{2}I_{12}^{2}(\omega,q)=0$, yields
the dispersion relation of the Anderson-Bogoliubov mode,
$\omega=\omega_{AB}(q)$, which has a phonon-like behavior at low momentum.

Although not explicitly written, all functions in Eq.(20) depend also on the
two mean-field parameters describing the BCS-BEC crossover in Fermi gases,
namely the energy gap $\Delta$ and the chemical potential $\mu$. Their values
can be obtained by inverting the BCS gap equation and the number equation
describing the BCS-BEC crossover for uniform systems, see for instance [8]. In
particular, for a unitary Fermi gas, corresponding to $1/(k_{F}a)=0$, one
finds $\mu=0.591E_{F}$ and $\Delta=0.686E_{F}$. The dynamic structure factor
can be calculated numerically from Eq.s (3) and (20). For a given momentum
$q$, there is a threshold frequency $\omega_{tr}(q)$, where the collective
mode merges with the continuum of single-particle excitations [8]

$\omega_{tr}(q)=\left\\{{2\Delta\;\;\;\;\;\;\;\;\;\;\textrm{for}\;\mu>0\;\text{and}\;q<2\sqrt{2m\mu},\atop
2\sqrt{(q^{2}/(8m)-\mu)^{2}+\Delta^{2}}\;\;\;\text{otherwise.}}\right.$ (21)

For $\omega<\omega_{tr}(q)$, all the functions $I^{\prime\prime},F$ and
$I_{ij}$ in Eq.(20) are real, due to the energy gap in the single-particle
energy spectrum. In this case the only contribution to the dynamic structure
factor originates from the collective mode, as the vanishing of the
denominator in the second term in Eq.(20) causes the appearance of an
imaginary part, via the formula $1/(x+i0^{+})=\textrm{P}1/x-i\pi\delta(x)$,
where P stands for the principal part. In contrast, for
$\omega>\omega_{tr}(q)$, the only contribution to $S(q,\omega)$ comes from the
continuum of single-particle excitations, with all the functions
$I^{\prime\prime},F$ and $I_{ij}$ in Eq.(20) acquiring an imaginary part by
themselves.

From the dynamic structure factor we extract the moments (5) and use them to
estimate the superfluid density and the associated Leggett bound in the
presence of a weak periodic potential, corresponding to $V_{0}\ll E_{F}$, with
$E_{F}=k_{F}^{2}/(2m)$ being the Fermi energy of the unperturbed system. Since
the potential induces a quadratic in $V_{0}$ corrections to both quantities,
we write $f_{S}=1-\alpha_{S}(V_{0}/E_{F})^{2}$ and
$f_{S}=1-\alpha_{S}^{L}(V_{0}/E_{F})^{2}$, where $\alpha_{S}$ and
$\alpha_{S}^{L}$ are dimensionless functions. From Eq.s (8) and (12) we find,
respectively, $\alpha_{S}^{L}(q)=2(m_{-1}(q))^{2}E_{F}^{2}$ and
$\alpha_{S}(q)=2m_{1}(q)m_{-3}(q)E_{F}^{2}$. In Fig.1 we plot such quantities
as a function of the ratio $q/k_{F}$ for a unitary Fermi gas. We see that in
general $\alpha_{S}^{L}(q)\neq\alpha_{S}(q)$, showing that the Leggett bound
is not exhausted in the presence of the external periodic potential. A notable
exception occurs for vanishing $q$, where the two quantities coincide and
reduce to $\kappa^{2}E_{F}^{2}/2=9E_{F}^{2}/(8\mu^{2})=3.221$. Moreover, we
see from Fig.1 that the corrections to the superfluid density become less and
less important as $q$ increases and vanish for infinite $q$, reflecting the
behavior of the static response function.

Figure 1: Lattice-induced perturbative corrections to the superfluid density,
$\alpha_{S}=(1-f_{S})E_{F}^{2}/V_{0}^{2}$, and to the associated Leggett’s
bound, $\alpha_{S}^{L}=(1-f_{S}^{L})E_{F}^{2}/V_{0}^{2}$, of a unitary Fermi
gas, calculated from the sum-rule approach (see Eq.s (12) and (8)).

In Fig.2 we plot the low-momentum estimates
$\omega(q)=\sqrt{m_{1}(q)/m_{-1}(q)}$ (green circles) and
$\omega(q)=\sqrt{m_{-1}(q)/m_{-3}(q)}$ (red diamonds) for the energy
dispersion of the Anderson-Bogoliubov collective mode of a unitary Fermi gas
based on the sum rule approach. The black solid line represents the full
numerical result, $\omega=\omega_{AB}(q)$, showing a positive, albeit small,
curvature [12]. At higher momentum, however, the correction to the phononic
dispersion becomes negative, as the frequency of the collective mode is
bounded above by the threshold frequency $\omega_{tr}$. We see that the sum
rule estimate $\omega(q)=\sqrt{m_{1}(q)/m_{-1}(q)}$ reproduces only the sound
velocity of the superfluid in the phononic regime, while missing the curvature
correction. In contrast, the sum-rule estimate
$\omega(q)=\sqrt{m_{-1}(q)/m_{-3}(q)}$ is more accurate as it includes the
first non linear correction.

Figure 2: Dispersion relation of the Anderson-Bogoliubov collective mode of
the unitary Fermi gas. The sum rule estimates (symbols) in the long wavelength
limit are compared with the full numerical result $\omega=\omega_{AB}(q)$
(black solid line), first discussed in Ref.[12] . Here $k_{F}$ is the Fermi
momentum and $v_{F}=k_{F}/m$ is the Fermi velocity. The sound velocity of the
system is $c=0.444v_{F}$.

###  III.2 B. BdG equations for systems in a periodic potential

So far we have considered a weak external periodic potential, so that its
effect on the transport properties of the unitary Fermi gas can be calculated
within perturbation theory. Below we compute the superfluid fraction and the
Leggett bound for a periodic potential of arbitrary strength, by solving
numerically the BdG equations. Since the external potential (2) is periodic,
the superfluid fraction defined in Eq. (9) coincides with the effective mass
ratio $m/m^{*}$ of the Fermi gas, calculated from the curvature of the lowest
energy band according to

$f_{S}=\frac{m}{m^{*}}\equiv\lim_{Q\to
0}\frac{1}{\bar{n}}\frac{\partial^{2}e(\bar{n},Q)}{\partial Q^{2}},$ (22)

where $e=E/V$ is the energy density and $Q$ is the wave-vector of the
supercurrent. To see this, we note that -$\langle P_{x}\rangle$ in Eq. (9)
represents the total momentum of the system along the x-axis measured with
respect to the moving frame, so that the supercurrent density $j$ in the lab
frame is given by

$j=\bar{n}v-\frac{\langle P_{x}\rangle}{mV}=f_{S}\bar{n}v+O(v^{2}).$ (23)

Since the supercurrent density is related to the wave-vector by $j=\partial
e/\partial Q$, from Eq. (23) we find $v=Q/m$ and

$e(\bar{n},Q)=e(\bar{n},0)+f_{S}\bar{n}\frac{Q^{2}}{2m}+O(Q^{3}),$ (24)

which is consistent with Eq. (22). In the absence of the lattice ($V_{0}=0$),
one has $e(\bar{n},Q)=e(\bar{n},0)+\bar{n}Q^{2}/(2m)$, and hence Eq.(22)
yields $m^{*}=m$ and $f_{S}=1$.

In order to calculate the superfluid fraction via Eq. (22), we proceed as in
Ref. [19, 20] by imposing an order parameter of the Bloch form
$\Delta(\mathbf{r})=e^{i2Qx}\tilde{\Delta}(x)$ in the BDG equations, where
$\tilde{\Delta}(x)$ is a periodic function with period $d$. This implies that
the corresponding eigenfunctions take the form
$u(\mathbf{r})=\tilde{u_{\mathbf{k}}}(x)e^{iQx}e^{i\mathbf{k}\cdot\mathbf{r}}/\sqrt{V}$
and
$v(\mathbf{r})=\tilde{v_{\mathbf{k}}}(x)e^{-iQx}e^{i\mathbf{k}\cdot\mathbf{r}}/\sqrt{V}$,
where $k_{x}$ lies in the first Brillouin zone and the amplitudes
$\tilde{u_{\mathbf{k}}}(x)$ and $\tilde{v_{\mathbf{k}}}(x)$ are periodic
functions with period $d$ obeying

$\left(\begin{array}[]{cc}H_{Q}(x)&\tilde{\Delta}(x)\\\
\tilde{\Delta}^{\ast}(x)&-H_{-Q}(x)\end{array}\right)\left(\begin{array}[]{c}\tilde{u_{\mathbf{k}}}(x)\\\
\tilde{v_{\mathbf{k}}}(x)\end{array}\right)=\epsilon_{\mathbf{k}}\left(\begin{array}[]{c}\tilde{u_{\mathbf{k}}}(x)\\\
\tilde{v_{\mathbf{k}}}(x)\end{array}\right)\;,$ (25)

where

$H_{Q}(x)\equiv\frac{1}{2m}\left[\mathbf{k}_{\perp}^{2}+\left(-i\partial_{x}+Q+k_{x}\right)^{2}\right]+V_{\rm
ext}(x)-\mu\,.$ (26)

For a given value of the momentum $\mathbf{k}$, the BdG equations (25)
represents an eigenvalue problem with the corresponding eigenfunctions being
normalized according to
$\int_{0}^{d}[\tilde{u}_{\mathbf{k}\alpha}^{*}(x)\tilde{u}_{\mathbf{k}\beta}(x)+\tilde{v}_{\mathbf{k}\alpha}^{*}(x)\tilde{v}_{\mathbf{k}\beta}(x)]dx=d\delta_{\alpha,\beta}$.

The local paring field $\tilde{\Delta}(x)$ and the chemical potential $\mu$,
appearing in Eq.s (25) and (26), are variational parameters to be determined
self-consistently. The first constraint is provided by the gap equation

$\tilde{\Delta}(x)=-g\frac{1}{V}\sum_{\mathbf{k}\alpha}\tilde{u}_{\mathbf{k}\alpha}(x)\tilde{v}_{\mathbf{k}\alpha}^{*}(x),$
(27)

where $g$ is the bare coupling constant of the contact interaction potential
and the sum is restricted over the solutions with
$\epsilon_{\mathbf{k}\alpha}\geq 0$. It is well known that the rhs of Eq.(27)
is actually plagued by a ultraviolet divergence and must be regularized.
Following the procedure suggested in Ref.[21], we introduce an explicit high-
energy cut-off in the sum over eigenenergies, and replace $g$ with a position-
dependent effective coupling constant. The second constraint to be satisfied
is that the spatial average of the density profile must yield the mean density
$\bar{n}$, that is $\bar{n}=\int_{0}^{d}n(x)dx/d$, where

$n(x)=\frac{1}{V}\sum_{\mathbf{k}\alpha}2|\tilde{v}_{\mathbf{k}\alpha}(x)|^{2}.$
(28)

Once the mean field parameters are known, we compute the energy density
$e=E/V$ of the system according to

$\displaystyle e$ $\displaystyle=$
$\displaystyle\frac{1}{V}\sum_{\mathbf{k}\alpha}\int_{0}^{d}\frac{dx}{d}\Big{[}2(\mu-\epsilon_{\mathbf{k}\alpha})|\tilde{v}_{\mathbf{k}\alpha}(x)|^{2}$
(29)
$\displaystyle+\tilde{\Delta}^{*}(x)\tilde{u}_{\mathbf{k}\alpha}(x)\tilde{v}_{\mathbf{k}\alpha}^{*}(x)\Big{]}.$

and use Eq.(22) to extract the superfluid fraction from the obtained result.
Notice that the two contributions in the rhs of Eq.(29) are separately
divergent but their sum is finite, implying that this equation is well-
defined.

For the Leggett’s bound we instead solve the BdG equations (25) at
equilibrium, i.e. with $Q=0$, and then extract $f_{S}^{L}$ from Eq.s (1) and
(28). In Fig.3 we display $f_{S}$ and $f_{S}^{L}$ as a function of the ratio
$V_{0}/\epsilon_{F}$ for the unitary Fermi gas, the two panels corresponding
to two different values of the dimensionless parameter $k_{F}d=4$ (a) and
$k_{F}d=10$ (b). We see that the Leggett bound is generally not saturated for
the unitary Fermi gas, but the difference tends to reduce as $k_{F}d$
increases. In particular in the limit $k_{F}d\gg 1$ the two values $f_{S}$ and
$f_{S}^{L}$ are expected to approach the prediction of the local density
approximation (LDA), where the density profile of the Fermi gas can be
calculated starting from the expression
$\mu_{0}=\mu_{\textrm{hom}}(n(x))+V_{0}\cos(qx)$ for the chemical potential,
with $\mu_{\textrm{hom}}(n)=(1+\beta)E_{F}=(1+\beta)(3\pi^{2}n)^{2/3}/(2m)$
the equation of state of the uniform Fermi gas at unitarity [22]. The density
profile in the LDA then takes the form

$n_{\textrm{LDA}}(x)=\frac{1}{3\pi^{2}}\left(\frac{2m}{1+\beta}\right)^{3/2}\Big{[}\mu_{0}-V_{0}\cos(2\pi
x/d)\Big{]}^{3/2}\;,$ (30)

yielding the prediction

$\frac{1}{f_{S}^{LDA}}=\langle
n\rangle\langle\frac{1}{n}\rangle=h\left(\frac{V_{0}}{\mu_{0}}\right)\int_{0}^{1}d\tilde{x}\Big{[}1-\frac{V_{0}}{\mu_{0}}\cos(2\pi\tilde{x})\Big{]}^{-3/2},$
(31)

for the inverse superfluid fraction, where $\tilde{x}=x/d$ and the function
$h$ is defined as
$h(V_{0}/\mu_{0})=\int_{0}^{1}d\tilde{x}[1-(V_{0}/\mu)\cos(2\pi\tilde{x})]^{3/2}$.
The chemical potential $\mu_{0}$ is related to the Fermi energy $E_{F}$ and
the potential strength $V_{0}$ via the normalization condition
$\bar{n}=\int_{0}^{d}n_{\textrm{LDA}}(x)dx/d$, yielding

$\mu_{0}=\frac{(1+\beta)E_{F}}{h(V_{0}/\mu_{0})^{2/3}}.$ (32)

Eq.s (31) and (30) show that in the LDA the superfluid density is independent
of the lattice period $d$ and only depends on the ratio $V_{0}/E_{F}$ (see
green solid line in Fig.3). In particular the superfluid fraction vanishes for
$V_{0}/E_{F}=(1+\beta)/h(1)^{2/3}=0.522$, corresponding to $V_{0}/\mu_{0}=1$,
where the integral in Eq.(31) diverges.

For small values of $V_{0}/E_{F}$ the LDA expression for the superfluid
density exhibits the quadratic dependence (see also [19] and Eq.(12) in the
$q\to 0$ limit)

$f_{S}^{LDA}=1-\frac{9}{8(1+\beta)^{2}}\frac{V_{0}^{2}}{E_{F}^{2}},$ (33)

It is interesting to compare the above expansion with the corresponding
expansion for the compressibility, in the same LDA regime ($q\to 0$ limit):

$\kappa^{LDA}=\frac{2}{3}(1+\beta)E_{F}\Big{[}1+\frac{1}{8(1+\beta)^{2}}\frac{V_{0}^{2}}{E_{F}^{2}}\Big{]},$
(34)

Eq.(34) is easily obtained by applying perturbation theory for the calculation
of the energy in the presence of the perturbation (2) and plays a crucial role
in the propagation of sound (see Sect IV below). Notice in particular that the
effect of the perturbation is smaller by a factor $9$, with respect to the
analogous expansion (33) for the superfluid fraction. It is also interesting
to compare the above results with the expansions holding in a dilute BEC gas
[7] where, in the LDA the $V_{0}^{2}$ correction to the compressibility is
exactly vanishing, as a consequence of the linear density dependence exhibited
by the equation of state.

In Fig.3 we also show the superfluid fraction (dotted line) and its Leggett’s
bound (dashed line) for a Fermi gas in the weakly interacting (WI) BCS regime,
corresponding to $\Delta\ll E_{F}$. In this limit the density profile and the
energy of the system can be calculated by assuming that the Fermi gas is non
interacting. Of course the non interacting Fermi gas is not superfluid because
of the occurrence of transverse gapless excitations. However, in the
superfluid phase the superfluid fraction $f_{S}$ coincides with the effective
mass (see (Eq.(22)). The evaluation of the effective mass in the non
interacting Fermi gas is expected to provide the actual value of $f_{S}^{WI}$,
if the condition $qk_{F}/m\gg\Delta$, corresponding to wave lengths of the
perturbation much smaller than the size of Cooper pairs, is satisfied. One can
then write [23]:

$f_{S}^{WI}\simeq\frac{1}{\bar{n}V}\sum_{\mathbf{k}j}\frac{\partial^{2}\epsilon_{j}(k_{x})}{\partial
k_{x}^{2}}2\Theta\Big{(}\mu-\epsilon_{j}(k_{x})-\frac{\mathbf{k}_{\perp}^{2}}{2m}\Big{)},$
(35)

where $\epsilon_{j}(k_{x})$ are the single-particle dispersion relations
induced by the 1D external potential, with $j$ being the band index, and the
integration over $k_{x}$ is restricted to the first Brillouin zone.

The Leggett bound in the WI regime can be obtained from Eq.(1) by using the
density profile of the non interacting Fermi gas

$n(x)=\frac{1}{\bar{n}V}\sum_{\mathbf{k}j}|\psi_{k_{x}j}(x)|^{2}2\Theta\Big{(}\mu-\epsilon_{j}(k_{x})-\frac{\mathbf{k}_{\perp}^{2}}{2m}\Big{)},$
(36)

where $\psi_{k_{x}j}(x)$ are the amplitudes of the Bloch states.

A first useful comparison between the superfluid fractions calculated at
unitarity and in the weekly interacting BCS gas concerns the behavior of
Leggett bound in the local density approximation (see Eq.(31)). Since the
chemical potential of the non interacting Fermi gas and of the gas at
unitarity have the same density dependence, except for the factor $(1+\beta)$
(see Eq.(32)), one finds that in the weakly interacting case $f_{S}^{LDA}$
always stays significantly above the corresponding value at unitarity. In
particular the LDA superfluid fraction vanishes for a larger value of the
ratio $V_{0}/E_{F}$, equal to $0.522/(1+\beta)=0.883$. The situation can be
however very different for the actual value $f_{S}$ of the superfluid
fraction. For example, for small values of $V_{0}$ the superfluid fraction of
the weakly interacting BCS gas always stays below the value at unitarity. For
large wave vectors (see Fig. 3a) the same happens also for larger values of
$V_{0}$, suggesting that interactions play a crucial role to facilitate
superfluid transport, as previously reported by Watanabe and Pethick with
important implications on the behavior of neutrons stars [20].

Another interesting feature emerging from Fig.3 (ab) is that in the WI regime
$f_{S}$ exhibits a linear dependence for small $V_{0}$. We have checked that
this behavior occurs only for $k_{F}d\geq\pi$, while for $k_{F}d<\pi$ we find
that $1-f_{S}\propto V_{0}^{2}$, as previously observed for the interacting
system. This effect can be simply understood by noting that for a
noninteracting uniform Fermi gas the moment $m_{-3}(q)$ in Eq.(12) diverges
for $q<2k_{F}$, because the dynamic structure factor $S(q,\omega)$ is linear
in $\omega$ at low frequency, reflecting the presence of gapless single-
particle excitations at the Fermi surface. In the superfluid Fermi gas these
excitations are instead gapped due to Cooper pairing. In the Appendix we
provide an explicit derivation of the behavior of $f_{S}$ for small $V_{0}$ in
the WI regime based on degenerate perturbation theory.

Figure 3: Superfluid fraction $f_{S}$ (blue triangles) and its Leggett’s bound
$f_{S}^{L}$ (red circles) for the unitary Fermi gas as a function of the ratio
$V_{0}/E_{F}$ for $k_{F}d=4$ (upper panel) and $k_{F}d=10$ (lower panel). The
green solid line refers to the local density approximation (LDA) at unitarity.
In each plot we also include the corresponding predictions for the weakly
interacting BCS regime for $f_{S}$ (black dotted line), $f_{S}^{L}$ (black
dashed line) and for LDA (green dashed line).

##  IV IV. Behaviour of the sound velocities

The results discussed in the previous Section on the superfluid density of the
interacting Fermi gas have important implications on the value of the sound
velocity and in particular on its anisotropy in the presence of the periodic
perturbation (2). According to the hydrodynamic theory of superfluids, the
square sound velocities along the longitudinal and transverse directions are
given by

$\displaystyle c^{2}_{x}$ $\displaystyle=$ $\displaystyle f_{S}\kappa^{-1}$
$\displaystyle c^{2}_{y}$ $\displaystyle=$ $\displaystyle
c^{2}_{z}=\kappa^{-1},$ (37)

emphasizing that their anisotropy, caused by the presence of the 1D periodic
potential, follows from the anisotropy of the superfluid fraction. From an
experimental point of view a natural strategy to determine the value of the
superfluid fraction, associated with the motion of the superfluid along the
$x$-direction, then consists of measuring the ratio

$\frac{c^{2}_{x}}{c^{2}_{y}}=f_{S}\;$ (38)

which is independent on the value of the compressibility. This procedure has
been recently successfully employed in the case of a dilute Bose gas confined
in a 2D box trap [7].

##  V V. CONCLUSIONS

In this work we have investigated the superfluid fraction of a Fermi gas at
unitarity, in the presence of a one-dimensional periodic potential. We have
shown that, in sharp contrast with the dilute BEC gas case, the Leggett’s
upper bound for the superfluid fraction can significantly overestimate its
actual value at zero temperature, unless the period $2\pi/q$ of the applied
lattice becomes much larger than the size of Cooper pairs, i.e. if
$q\ll\Delta/(mk_{F})$, where $\Delta$ is the pairing gap energy. We have
pointed out the occurrence of important deviations of Leggett’s bound from the
actual value $f_{S}$ of the superfluid fraction at unitarity and shown that
these deviations becomes more and more important as one approaches the weakly
interacting BCS regime of small and negative values of the $s$-wave scattering
length. Using sum rule arguments we have also shown that the deviations of the
Leggett bound from the superfluid fraction explain the change of sign of the
curvature of the Anderson-Bogoliubov mode dispersion taking place around
unitarity.

We thank Y. Castin, F. Dalfovo, S. Giorgini, A. Leggett and W. Zwerger for
fruitful discussions. This work was granted access to the HPC resources of
TGCC under the allocations AD010513635 and AD010513635R1 supplied by GENCI
(Grand Equipement National de Calcul Intensif). S.S. acknowledges continuous
financial support from PAT (Provincia Autonoma di Trento.)

##  VI APPENDIX

For completeness, in this Appendix we provide an analytical calculation of
$f_{S}$ for a weak periodic potential in the weakly interacting BCS regime,
under the assumption that $q\gg\Delta m/k_{F}$, which is equivalent to
$k_{F}d\ll E_{F}/\Delta$. As discussed in the main text, under this condition
one is allowed to identify the inverse effective mass, calculated using the
ideal Fermi gas, with the superfluid fraction of the weakling interacting BCS
superfluid. Below we will consider separately two cases, $k_{F}d<\pi$ and
$k_{F}d\geq\pi$.

###  VI.1 Case $k_{F}d<\pi$

We calculate $f_{S}^{WI}$ from Eq.(12) since $m_{-3}(q)$ is expected to be
finite for a uniform system. For $q>2k_{F}$ the dynamic structure factor of
the noninteracting Fermi gas is given by

$S(\omega,q)=\frac{m}{4\pi^{2}q}\left[k_{F}^{2}-\frac{m^{2}}{q^{2}}\left(\omega-\frac{q^{2}}{2m}\right)^{2}\right],$
(39)

if the rhs of Eq. (39) is positive, that is for
$|\omega-q^{2}/(2m)|<qk_{F}/m$, and zero otherwise.

$\displaystyle m_{-3}(q)$ $\displaystyle=$
$\displaystyle\frac{m}{4\pi^{2}q\bar{n}}\int_{q^{2}/(2m)-qk_{F}/m}^{q^{2}/(2m)+qk_{F}/m}\frac{d\omega}{\omega^{3}}$
(40)
$\displaystyle\times\left[k_{F}^{2}-\frac{m^{2}}{q^{2}}\left(\omega-\frac{q^{2}}{2m}\right)^{2}\right]$

To perform the integral, we introduce dimensionless variable
$y=[\omega-q^{2}/(2m)]/E_{F}$ and $p=q/k_{F}$, so that the integral becomes

$m_{-3}(q)=\frac{m^{3}}{\pi^{2}k_{F}^{3}\bar{n}p}\int_{-2p}^{2p}\left(1-\frac{y^{2}}{4p^{2}}\right)\frac{dy}{(y+p^{2})^{3}}$
(41)

By performing the integration in Eq.(41) and recalling that
$\bar{n}=k_{F}^{3}/(3\pi^{2})$, we obtain $m_{-3}(q)=F(q/k_{F})/E_{F}^{3}$,
where

$F(p)=\frac{3}{8p^{3}}\left[\ln\left(\frac{p-2}{p+2}\right)+\frac{4p}{p^{2}-4}\right].$
(42)

Since $m_{1}(q)=q^{2}/(2m)$, from Eq.(12) we find

$1-f_{S}^{WI}\simeq
2\left(\frac{q}{k_{F}}\right)^{2}F\left(\frac{q}{k_{F}}\right)\left(\frac{V_{0}}{E_{F}}\right)^{2}.$
(43)

Notice that the function $F(p)$ in Eq.(42) diverges at $p=2$ as
$F(p)\approx\frac{3}{128(p-2)}$, implying that for $q=2k_{F}$ the correction
to the superfluid fraction induced by the periodic potential are no longer
quadratic but linear in $V_{0}$ , as we shall see below.

###  VI.2 Case $k_{F}d\geq\pi$

We use degenerate perturbation theory applied to the external periodic
potential (2) to show that $1-f_{S}\propto V_{0}$. For simplicity, we will
prove the claim for the special case $k_{F}d=\pi$, so that only the lowest
energy band is populated. The proof can be easily adapted to the general case
where higher bands are also populated.

We consider two states with wave-vector $k_{x}$ and $k_{x}-G$, where $G$ is a
reciprocal lattice point. The two energy bands are obtained by the condition

$\begin{vmatrix}\frac{k_{x}^{2}}{2m}-E&\frac{V_{0}}{2}\\\
\frac{V_{0}}{2}&\frac{(k_{x}-G)^{2}}{2m}-E\end{vmatrix}=0,$ (44)

whose solutions are

$\displaystyle E_{\pm}(k_{x})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\frac{k_{x}^{2}}{2m}+\frac{(k_{x}-G)^{2}}{2m}\right)$
(45)
$\displaystyle\pm\left[\frac{1}{4}\left(\frac{k_{x}^{2}}{2m}-\frac{(k_{x}-G)^{2}}{2m}\right)^{2}+\frac{V_{0}^{2}}{4}\right]^{1/2}.$

We consider $0\leq k_{x}\leq\pi$ and choose $G=2\pi/d$ in Eq.(45) so that
$E_{\pm}(\pi/d)=E_{R}\pm V_{0}/2$. For $k_{F}=\pi/d$, the Fermi energy
$E_{F}=E_{R}$ sits between the two bands, so the upper band is completely
empty and according to in Eq.(35) it does not contribute to the superfluid
fraction. Let us introduce the dimensionless variables
$\tilde{E}_{-}=E_{-}/E_{R}$ and $x=k_{x}\pi/d$. Then

$\tilde{E}_{-}(x)=x^{2}+2-2x-[4(1-x)^{2}+V_{r}^{2}]^{1/2},$ (46)

where $V_{r}=V_{0}/(2E_{R})$ is the rescaled potential strength. From Eq. (46)
we find

$\frac{\partial^{2}\tilde{E}_{-}(x)}{\partial
x^{2}}=2-\frac{4V_{r}^{2}}{[4(1-x)^{2}+V_{r}^{2}]^{3/2}},$ (47)

showing that the curvature of the bands diverges as $2-4/V_{r}$ for vanishing
$V_{r}$. Substituting Eq.s (46) and (47) into Eq.(35), we find

$f_{S}^{WI}\simeq
1-\frac{4V_{r}^{2}\int_{0}^{1}\frac{[-(x-1)^{2}+\sqrt{4(1-x)^{2}+V_{r}^{2}}]}{[4(1-x)^{2}+V_{r}^{2}]^{3/2}}dx}{2\int_{0}^{1}[-(x-1)^{2}+\sqrt{4(1-x)^{2}+V_{r}^{2}}]dx}.$
(48)

Let us discuss the behavior of the two integrals in Eq.(48), starting from the
upper one. Since the major contribution to the integral comes from the region
around $x=1$, the term proportional to $-(x-1)^{2}$ in the numerator can be
neglected as it vanishes for $x=1$ (its contribution to the integral is of
order $V_{r}^{2}\ln V_{r}$). For the lower integral, we can safely set
$V_{r}=0$ as no singularity is present, so that Eq. (48) reduces to

$f_{S}^{WI}\simeq
1-\frac{4V_{r}^{2}\int_{0}^{1}[4(1-x)^{2}+V_{r}^{2}]^{-1}dx}{2\int_{0}^{1}[-(x-1)^{2}+2(1-x)]dx}.$
(49)

Since $V_{r}$ is small, in the upper integral we can safely replace $0$ with
$-\infty$, so that the integration becomes elementary, yielding

$f_{S}^{WI}\simeq 1-\frac{\pi
V_{r}}{4/3}=1-\frac{3\pi}{8}\frac{V_{0}}{E_{F}},$ (50)

showing that the correction to the superfluid fraction is linear in $V_{0}$,
which agrees with our numerics.

## References

  * Leggett [1970] A. J. Leggett, Can a solid be ”superfluid”?, Phys. Rev. Lett. 25, 1543 (1970).
  * [2] In general the superfluid density (and hence the superfluid fraction) is an anisotropic tensor and in the following we will adopt the notation $f_{s}$ to indicate the superfluid fraction along the $x$-axis.
  * Peshkov [1944] V. P. Peshkov, Second sound in helium II, J. Phys. (Moscow) 8, 381 (1944).
  * Sidorenkov  _et al._ [2013] L. A. Sidorenkov, M. K. Tey, R. Grimm, Y.-H. Hou, L. Pitaevskii, and S. Stringari, Second sound and the superfluid fraction in a fermi gas with resonant interactions, Nature 498, 78 (2013).
  * Christodoulou  _et al._ [2021] P. Christodoulou, M. Gałka, N. Dogra, R. Lopes, J. Schmitt, and Z. Hadzibabic, Observation of first and second sound in a bkt superfluid, Nature 594, 191 (2021).
  * Leggett [1998] A. J. Leggett, On the superfluid fraction of an arbitrary many-body system at $T=0$, J. Stat. Phys. 93, 927 (1998).
  * Chauveau  _et al._ [2023] G. Chauveau, C. Maury, F. Rabec, C. Heintze, G. Brochier, S. Nascimbene, J. Dalibard, J. Beugnon, S. M. Roccuzzo, and S. Stringari, Superfluid fraction in an interacting spatially modulated bose-einstein condensate, Phys. Rev. Lett. 130, 226003 (2023).
  * Combescot  _et al._ [2006] R. Combescot, M. Y. Kagan, and S. Stringari, Collective mode of homogeneous superfluid fermi gases in the bec-bcs crossover, Phys. Rev. A 74, 042717 (2006).
  * Pines and Nozières [1966] D. Pines and P. Nozières,  _The Theory of Quantum Liquids: Normal Fermi liquids_, Advanced book classics series (W.A. Benjamin, 1966).
  * Son and Wingate [2006] D. Son and M. Wingate, General coordinate invariance and conformal invariance in nonrelativistic physics: Unitary fermi gas, Annals of Physics 321, 197 (2006), january Special Issue.
  * Haussmann  _et al._ [2009] R. Haussmann, M. Punk, and W. Zwerger, Spectral functions and rf response of ultracold fermionic atoms, Phys. Rev. A 80, 063612 (2009).
  * Kurkjian  _et al._ [2016] H. Kurkjian, Y. Castin, and A. Sinatra, Concavity of the collective excitation branch of a fermi gas in the bec-bcs crossover, Phys. Rev. A 93, 013623 (2016).
  * Zou  _et al._ [2018] P. Zou, H. Hu, and X.-J. Liu, Low-momentum dynamic structure factor of a strongly interacting fermi gas at finite temperature: The goldstone phonon and its landau damping, Phys. Rev. A 98, 011602 (2018).
  * Biss  _et al._ [2022] H. Biss, L. Sobirey, N. Luick, M. Bohlen, J. J. Kinnunen, G. M. Bruun, T. Lompe, and H. Moritz, Excitation spectrum and superfluid gap of an ultracold fermi gas, Phys. Rev. Lett. 128, 100401 (2022).
  * Carlson  _et al._ [2003] J. Carlson, S.-Y. Chang, V. R. Pandharipande, and K. E. Schmidt, Superfluid fermi gases with large scattering length, Phys. Rev. Lett. 91, 050401 (2003).
  * Astrakharchik  _et al._ [2004] G. E. Astrakharchik, J. Boronat, J. Casulleras, and S. Giorgini, Equation of state of a fermi gas in the bec-bcs crossover: A quantum monte carlo study, Phys. Rev. Lett. 93, 200404 (2004).
  * Ku  _et al._ [2012] M. J. H. Ku, A. T. Sommer, L. W. Cheuk, and M. W. Zwierlein, Revealing the superfluid lambda transition in the universal thermodynamics of a unitary fermi gas, Science 335, 563 (2012).
  * Minguzzi, A.  _et al._ [2001] Minguzzi, A., Ferrari, G., and Castin, Y., Dynamic structure factor of a superfluid fermi gas, Eur. Phys. J. D 17, 49 (2001).
  * Watanabe  _et al._ [2008] G. Watanabe, G. Orso, F. Dalfovo, L. P. Pitaevskii, and S. Stringari, Equation of state and effective mass of the unitary fermi gas in a one-dimensional periodic potential, Phys. Rev. A 78, 063619 (2008).
  * Watanabe and Pethick [2017] G. Watanabe and C. J. Pethick, Superfluid density of neutrons in the inner crust of neutron stars: New life for pulsar glitch models, Phys. Rev. Lett. 119, 062701 (2017).
  * Bulgac and Yu [2002] A. Bulgac and Y. Yu, Renormalization of the hartree-fock-bogoliubov equations in the case of a zero range pairing interaction, Phys. Rev. Lett. 88, 042504 (2002).
  * Pitaevskii and Stringari [2016] L. Pitaevskii and S. Stringari,  _Bose-Einstein condensation and superfluidity_ (Oxford University Press, 2016).
  * Pitaevskii  _et al._ [2005] L. P. Pitaevskii, S. Stringari, and G. Orso, Sound propagation and oscillations of a superfluid fermi gas in the presence of a one-dimensional optical lattice, Phys. Rev. A 71, 053602 (2005).


* Spitler et al. (2016) Spitler, L. G., Scholz, P., Hessels, J. W. T., et al. 2016, Nature, 531, 202, doi: 10.1038/nature17168
  * Swat et al. (2016) Swat, M. J., Grenon, P., & Wimalaratne, S. 2016, Bioinformatics, 32, 2719, doi: 10.1093/bioinformatics/btw170
  * Tendulkar et al. (2021) Tendulkar, S. P., Gil de Paz, A., Kirichenko, A. Y., et al. 2021, ApJ, 908, L12, doi: 10.3847/2041-8213/abdb38
  * The CHIME Collaboration et al. (2022) The CHIME Collaboration, Amiri, M., Bandura, K., et al. 2022, arXiv e-prints, arXiv:2201.07869.  https://arxiv.org/abs/2201.07869
  * Thornton et al. (2013) Thornton, D., Stappers, B., Bailes, M., et al. 2013, Science, 341, 53, doi: 10.1126/science.1236789
  * Wilks (1938) Wilks, S. S. 1938, The Annals of Mathematical Statistics, 9, 60 , doi: 10.1214/aoms/1177732360
  * Yao et al. (2017) Yao, J. M., Manchester, R. N., & Wang, N. 2017, ApJ, 835, 29, doi: 10.3847/1538-4357/835/1/29
  * Zhang & Zhang (2022) Zhang, R. C., & Zhang, B. 2022, ApJ, 924, L14, doi: 10.3847/2041-8213/ac46ad
  * Zyla et al. (2020) Zyla, P. A., Barnett, R. M., Beringer, J., et al. 2020, Progress of Theoretical and Experimental Physics, 2020, 083C01, doi: 10.1093/ptep/ptaa104