\section{Introduction}

Recommender systems, which support many real-world applications, have become increasingly accurate in recent years.
Important improvements are brought about by the use of deep learning \cite{cheng2016wide,guo2017deepfm} and the leverage of user historical data~\cite{din,dien}. In these systems, embedding learning (i.e., learning vector representations) for users and items plays a critical role.

The standard method, \emph{full embedding}, learns a unique vector representation for each user and item. However, this approach faces important limitations in practice.
Firstly, for cold users and items that have insufficient historical interactions (which commonly exist in practice), it is hard for them to learn an embedding that generalizes well on their own.
Secondly, web-scale recommender systems usually involve an enormous number of users and items; the memory cost incurred by full embedding can be prohibitively high: 1 billion users, each assigned a $64$ dimensional embedding in $32$-bit floating point, require $238$ GB of memory.

To mitigate the problem of cold users and items, a common idea is to impose regularization, e.g., forcing all embeddings to follow a given clustering structure \cite{eTREE,jnkm}. However, using a fixed predefined clustering structure lacks flexibility; the given clustering structure can hardly be optimal. Besides, existing methods \cite{eTREE,jnkm} require training in full embedding, still suffering from high memory costs.

To mitigate the problem of high memory cost, hashing \cite{modulo,binaryhash,hashgnn,prehash} is typically adopted to map entities to a fixed-size embedding table. This can directly tackle the problem of memory cost. However, forcing a hashing that is essentially random \cite{modulo,binaryhash} leads to unsatisfactory performance; the performances of current learning-to-hash methods \cite{hashgnn,prehash} also have a large space for improvements. Besides, pretraining a hash function \cite{hashgnn} or binding each hash bucket with a warm user \cite{prehash} also makes these methods lack flexibility.

In chasing a more flexible and integrated solution to these two problems, while at the same time improving the performance, we propose in this paper Clustered Embedding Learning (CEL). CEL achieves dynamic and automatic clustering of users and items, which gives it the potential to have superior performance. In CEL, entities within the same cluster may use and jointly learn a shared embedding; thus it can significantly reduce the memory cost.

CEL performs clustering in a top-down (i.e., divisive) fashion, starting with a single cluster that includes all users (or items) and recursively splitting a cluster into two. Given that we want to share information among similar entities, such divisive clustering is a natural choice: sharing information of all entities at the beginning, then seeking to refine the clustering when possible.

Though being a natural thought, implementing divisive clustering for embedding learning is challenging. In the circumstance of \textit{shared embedding}, it appears to have no standard to differentiate entities within the same cluster. To tackle this problem, we propose to refer to their gradients with respect to the shared embedding, since the differences in the gradients indicate their different desired optimization directions of the embedding. Then naturally, we propose to split the cluster along the first principal component of the gradients.

After cluster split and associated embedding optimizations, the current cluster assignment might have room for improvement, e.g., the embedding of some other cluster might better fit the interactions of an entity. Thus, fine-grained automatic cluster optimization also requires enabling cluster reassignment. Given that the training loss appears to be the only referable metric, we propose to conduct cluster reassignment based on direct loss minimization.

With such a cluster split and reassignment paradigm, CEL is ideally a sustainable framework that can not only support starting with zero data, but also can automatically increase the number of clusters, as well as automatically optimize the cluster reassignments, as the interaction data increases and the number of users/items increases.

\begin{figure*}[ht]
    \centering
    \includegraphics[width=0.96\linewidth]{rankingf.png}
    \caption{The standard behavior prediction paradigm in recommender system.}
    \label{fig:rankframework}
\end{figure*}

We validate the proposed CEL both theoretically and empirically. Firstly, we prove the identifiability and the existence of a unique optimal number of clusters for CEL in the context of nonnegative matrix factorization. Then, we show that the accelerated version of CEL (CEL-Lite) has an optimal time complexity, which allows efficient online updates. Finally, we conduct extensive experiments on public real-world datasets as well as a private business dataset, integrated with various feature interaction models (including NMF, DIN, and DeepFM); the results show that CEL consistently outperforms current state-of-the-art methods. In particular, when incorporating
CEL into the business model, it brings an improvement of $+0.6\%$ in AUC, which translates into a significant revenue gain; meanwhile, the size of the embedding table gets 2650 times smaller.

\section{Related Works}

One line of research regularizes the embeddings to follow a given clustering structure to improve their qualities.
For example, Joint NMF and K-Means (JNKM) \cite{jnkm} forces the embeddings to follow a K-clustering structure; eTree \cite{eTREE} learns the embeddings under an implicit tree (i.e., hierarchical clustering). Prior works in this line also include {HSR} \cite{HSR}, {HieVH} \cite{sun2017exploiting}, etc.
The major distinction of our work is that CEL does not predefine the clustering structure. Instead, it automatically optimizes the number of clusters and the cluster assignment. Besides, these methods require training in full embedding size, while CEL is trained with reduced embedding size.

Another line of work uses the hash to share and jointly learn the embeddings. The plain version applies the modulo operation \cite{modulo} on the hash ID of each entity to obtain its hash code (i.e., the embedding index). Sharing embedding within each hash bucket, the embedding size can be reduced by the magnitude of the modulo divisor.
Binary Hash (BH) \cite{binaryhash,zhang2016discrete} adopts a code block strategy on binarized hash IDs to reduce direct collisions.
Adaptive Embedding (AE) \cite{DeepRec} explicitly avoids collisions for warm users who have sufficient historical data.
However, the hash assignments in these methods are essentially random and will not adapt/learn during learning.
There are also several works that tried to learn a hash function. For example, HashGNN \cite{hashgnn} learns a binary code hash function with graph neural networks. However, it requires pretraining a hash function.
Preference Hash (PreHash) \cite{prehash} learns to hash users with similar preferences to the same bucket. PreHash binds a warm user to each bucket, while the embedding of other users is a weighted sum of the embeddings of their top-$K$ related buckets. By contrast, CEL does not force bindings between users and clusters.

Some other works seek to reduce memory cost by reducing the embedding dimension or the number of bits for each dimension (i.e., quantization) \cite{Mixeddimension,Automatedembedding,amtl,zhang2016discrete,yang2019quantization}, which are orthogonal to our work in the perspective of reducing the memory cost.
The type of clustering algorithms most relevant to CEL is Divisive Hierarchical Clustering (DHC) \cite{dhcreview}.
Particularly, Principal Direction Divisive Partitioning (PDDP) \cite{pddp} splits clusters via principal coordinate analysis \cite{PCoA}, which is remotely related to our proposed GPCA (Section~\ref{sec:gpca}).

\begin{figure*}[ht]
    \centering
    \includegraphics[width=0.98\linewidth]{framework.png}
    \\

    \hspace{35mm}(a)\hspace{80mm}(b)

    \caption{(a) An overview of CEL: it alternatively performs embedding optimization and cluster optimization. (b) The illustration of GPCA: it splits clusters along the first principal component of item gradients.}
    \label{fig:celframework}
\end{figure*}
\section{The Proposed Framework} \label{sec:methodology}

In this section, we present Clustered Embedding Learning (CEL) as a plug-and-play embedding learning framework which can be integrated with any differentiable feature interaction model such as Nonnegative Matrix Factorization (NMF), NeuMF~\cite{neuralCF}, Wide\&Deep \cite{cheng2016wide}, DeepFM~\cite{guo2017deepfm}, MMoE~\cite{ma2018modeling}, DIN~\cite{din} and DIEN~\cite{dien}. We illustrate the role of embedding learning in the behavior prediction  that is commonly used in recommender systems in Figure~\ref{fig:rankframework}.

Assume we have a user-item interaction data matrix $\mathbf{X}\in\mathbb{R}^{N\times M}$ of $N$ users and $M$ items, where $\mathbf{X}(i,j)$ denotes the interaction (rating, click, etc.) of the $i^{th}$ user on the $j^{th}$ item. We aim to recover the data matrix with two low-rank matrices, i.e.,  $\mathbf{X}\approx Y(\mathbf{A},\mathbf{B})$, where $\mathbf{A}\in\mathbb{R}^{N\times R}$ and $\mathbf{B}\in\mathbb{R}^{M\times R}$ can be regarded as the user embedding matrix and the item embedding matrix, respectively. Each row of $\mathbf{A}$ (or $\mathbf{B}$) corresponds to the embedding of a user (or item). $R$ is the embedding dimension.
The prediction $Y(\mathbf{A},\mathbf{B})=[y_{ij}]=[y(\mathbf{A}_i,\mathbf{B}_j)]\in\mathbb{R}^{N\times M}$ is produced by a differentiable feature interaction model $y$ that takes the embedding of user $i$ and the embedding of item $j$ as inputs.
Generally, embedding learning can be formulated as an optimization problem:
\begin{equation}

    \argmin_{\mathbf{A},\mathbf{B}} \; \mathcal{L},
    \;
    \mathcal{L}=\big\|\mathbf{W}\odot(\mathbf{X}-Y(\mathbf{A},\mathbf{B}))\big\|_{F}^2 ,
\label{equ:nmf_objective}
\end{equation}
where $\mathbf{W} \in \{0, 1\}^{N\times M}$ contains 1s at the indices of observed entries in $\mathbf{X}$, and 0s otherwise. $F$ denotes the Frobenius norm and $\odot$ denotes element-wise multiplication.

CEL clustering can be similarly applied to users and items. For clarity, we will only discuss the clustering of items. To initiate the clustering of items, we decompose the item embedding matrix as a product following a clustering structure:
\begin{equation}
    \mathbf{B}=\mathbf{S}_q\mathbf{B}_q,\;\;q=\{0, 1, 2, \ldots\},
\label{equ:sqbq}
\end{equation}
where $\mathbf{S}_q\in \{0, 1\}^{M\times M_q}$ and $\mathbf{B}_q\in\mathbb{R}^{M_q\times R}$ denote the cluster assignment matrix and the cluster embedding matrix of items, respectively. $M_q$ denotes the current number of clusters and $q$ denotes the current number of splits. We adopt hard cluster assignment: each row of $\mathbf{S}_q$ has exactly one $1$, while the other entries are zero.

CEL solves Eq.~\eqref{equ:nmf_objective} by alternatively performing embedding optimization and cluster optimization (Figure~\ref{fig:celframework}).

\subsection{Embedding Optimization}

In embedding optimization, we optimize $\mathbf{A}$ and $\mathbf{B}_q$ subject to a fixed $\mathbf{S}_q$.
The optimization of $\mathbf{A}$ can be directly conducted with Eq. \eqref{equ:nmf_objective} and \eqref{equ:sqbq}. Since hard cluster assignment is adopted, the embedding of each cluster can be optimized separately.

The loss of the $k^{th}$ cluster embedding of $\mathbf{B}_q$ can be denoted as:
\begin{equation}
\begin{split}
    \mathcal{L}_k=\big\|\mathbf{W}\odot\big[&\mathbf{X}\cdot \text{diag}(\mathbf{S}_q(:,k))
    \\
    &-Y(\mathbf{A},\mathbf{S}_q(:,k)\mathbf{B}_q(k,:))\big]\big\|_F^2,
\end{split}
\end{equation}
where colon $(:)$ stands for the slicing operator that takes all values of the respective field, and $\text{diag}(\cdot)$ maps a vector into a diagonal matrix. It has
$\mathcal{L}=\sum_{k=1}^{M_q}\mathcal{L}_k$.

\subsection{Cluster Optimization}

In this section, we introduce the two proposed cluster optimization operations: cluster reassignment and cluster split. We assume that the reassignment occurs every $t_1$ steps of embedding optimization, while the split occurs every $t_2$ steps of embedding optimization.

\subsubsection{Cluster Reassignment}

An item assigned to a sub-optimal cluster can result in performance loss. When the embeddings have been optimized, we may perform cluster reassignment so that an item can be assigned to a better cluster. Reassignment is an important step toward automatically grouping similar items.

One problem CEL faces is how to measure the distance between an item and a cluster, since their representations are in different domains: items have interaction data, while clusters have cluster embeddings. A natural metric seems to be the fitness of the cluster embedding to the item. Thus, we seek to directly minimize the loss of each item w.r.t. their cluster assignments:
\begin{equation}
    \mathbf{S}_q =
    \argmin_{\mathbf{S}_q} \; \big\|\mathbf{W}\odot(\mathbf{X}-Y(\mathbf{A},\mathbf{S}_q\mathbf{B}_q))\big\|_{F}^2.
\label{equ:reassignment}
\end{equation}
With fixed embeddings, the reassignment of each item is independent, therefore Eq. \eqref{equ:reassignment} can also be equivalently formulated as a set of sub-problems, one for each item:
\[
\mathbf{S}_q(j,:) =
\argmin_{\mathbf{S}_q(j,:)} \; \big\|\mathbf{W}(:,j)\odot(\mathbf{X}(:,j)-Y(\mathbf{A},\mathbf{S}_q(j,:)\mathbf{B}_q))\big\|_{F}^2.
\]

\subsubsection{Cluster Split with GPCA} \label{sec:gpca}

With embeddings optimized, we may consider cluster split to achieve a finer clustering. Each time we select one cluster and split it into two. There are many heuristic criteria for choosing the cluster. For example, we may choose the cluster that has the largest: total loss, mean loss, number of associated items, or number of associated interactions ($\|\mathbf{W}\mathbf{S}_q(:,k)\|_0$). According to our experiments, they generally lead to similarly good results. This demonstrates the robustness of CEL upon cluster-choosing criteria. In this paper, we adopt the number of associated interactions as the default criterion.

The key challenge for CEL is how to split a cluster. With hard cluster reassignment and embedding sharing, embeddings within a cluster are the same. Thus we need a representation that can differentiate items in the same cluster.
To tackle this problem, we propose Gradient Principal Component Analysis (GPCA) that splits clusters according to items' gradients w.r.t. their shared embedding, which indicates their different desired optimization directions of the embedding.
Such a split strategy can help achieve finer-grained clustering of similar items while further minimizing the loss.

Specifically, we assemble the gradients of items of a given cluster into a matrix $\mathbf{G}=[\ldots, g_j,\ldots]^\top$, $\forall j\in\{j|\mathbf{S}_q(j,k)=1\}$, where $g_j=-\partial\mathcal{L}/\partial \mathbf{B}(j,:)$ and $k$ denotes the cluster to be split.
Then, we normalize $\mathbf{G}$ (as a standard procedure before performing PCA, without changing its notation) and compute its first principal component:
\begin{equation}
\mathbf{p}=\argmax_{\|\mathbf{p}\|_2=1}\;\mathbf{p}^\top\mathbf{G}^\top\mathbf{G}\mathbf{p}.
\label{equ:pcavector}
\end{equation}
After that, we split the cluster into two according to their first principal component scores $g_j^\top\mathbf{p}$:
\begin{equation}
\hspace{-2mm}
    \left\{
        \begin{array}{ll}
          \mathbf{S}_{q+1}(j,k)=1,
          \;\;\;\;\;\;\;\;\;\; \text{ if }g_j^\top\mathbf{p}< \delta; \\
          \mathbf{S}_{q+1}(j,M_q+1)=1,
          \;\; \text{ if }g_j^\top\mathbf{p}\geq \delta.
        \end{array}
    \right.
\label{equ:pca score}
\end{equation}
One set of items is still assigned to cluster $k$, while the other set is assigned to a new cluster $M_q+1$. We may adjust the threshold $\delta$ to control/balance the number of items in each cluster.

We present the pseudocode of CEL in Algorithm \ref{alg:cel}.

\subsection{Theoretical Analysis of CEL under NMF}
\label{sec:identifiability}

CEL can be integrated with any differentiable feature interaction model, however, most of them (e.g., when involving other features or integrated with neural networks) are hard to analyze. Even so, we are able to study the theoretical properties of CEL when it is integrated with Nonnegative Matrix Factorization (NMF). These analyses can to some extent demonstrate the validity of CEL.

In NMF, $Y(\mathbf{A}, \mathbf{B}) = \mathbf{A} \mathbf{B}^\top$ and $\mathbf{A}, \mathbf{B}\geq 0$, which factorizes the data matrix $\mathbf{X}$ into the product of two low-rank matrices. Despite that low-rank matrix factorization generally has no unique solution, nonnegativity helps establish the identifiability of NMF. NMF can thus produce essentially unique low-dimensional representations \cite{lee1999learning}. The NMF of $\mathbf{X} = \mathbf{A} \mathbf{B}^\top$ is regarded as \emph{essentially unique} if the nonnegative matrix $\mathbf{A}$ and $\mathbf{B}$ are identifiable up to a permutation and scaling. We show that the solution of CEL is identifiable:

\begin{theorem}
(Identifiability of CEL) Assume that the data matrix $\mathbf{X}$ follows $\mathbf{X} = \mathbf{A} \mathbf{B}_q^\top \mathbf{S}_q^\top$, where $\mathbf{A} \in\mathbb{R}^{N\times R}$, $\mathbf{B}_q\in \mathbb{R}^{M_q\times R}$, $\mathbf{S}_q \in \{0,1\}^{M\times M_q}$, $\|\mathbf{S}_q(j,:)\|_0 = 1,\forall j\in [1,M]$, $rank(\mathbf{X}) = rank(\mathbf{A}) = R$, $\mathbf{B}_q$ has no repeated rows and, without loss of generality, $M_q\geq R$. If $\mathbf{S}_q$ is of full-column rank and rows of $\mathbf{B}_q$ are sufficiently scattered, then $\mathbf{A}$, $\mathbf{B}_q$ and $\mathbf{S}_q$ are essentially unique.
\label{theorem:CELIdentifiability}
\end{theorem}

The formal definitions (e.g., sufficiently scattered) \cite{eTREE,jnkm,identifiability} and detailed proofs are provided in Appendix
B.

We assume $\mathbf{S}_q$ to be of full-column rank, which means there is no empty cluster. Enforcing such a condition in our implementation has improved the model performance.
The assumption of $\mathbf{B}_q$ being sufficiently scattered is common and is likely to be satisfied as shown by \cite{eTREE}.
Based on Theorem~\ref{theorem:CELIdentifiability}, we can further deduce the existence of a unique optimal number of clusters:

\begin{proposition}
(Optimal Number of Clusters of CEL) Under the assumptions of Theorem \ref{theorem:CELIdentifiability} and further assume that the data matrix can also be decomposed as $\mathbf{X} = \mathbf{A}' \mathbf{B}_{q+1}^\top \mathbf{S}_{q+1}^\top$, where $\mathbf{A}' \in\mathbb{R}^{N\times R}$, $\mathbf{B}_{q+1}\in \mathbb{R}^{M_{q+1}\times R}$, $\mathbf{S}_{q+1} \in \{0,1\}^{M\times M_{q+1}}$, $\|\mathbf{S}_{q+1}(j,:)\|_0 = 1,\forall j\in [1,M]$, $rank(\mathbf{A}') = R$, and $M_{q+1}>M_{q}\geq R$. If $\mathbf{S}_{q+1}$ is of full-column rank and the rows of $\mathbf{B}_{q+1}$ are sufficiently scattered, then $\mathbf{B}_{q+1}$ has only $M_{q}$ distinct rows, which can be obtained via permutation and scaling of the rows in $\mathbf{B}_q$.
\label{proposition:clusternumber}
\end{proposition}

The proof is built on the fact that the identifiability of CEL implies $\mathbf{S}_{q+1}\mathbf{B}_{q+1}$ is a permutation and scaling of $\mathbf{S}_{q}\mathbf{B}_{q}$. The implication of Proposition \ref{proposition:clusternumber} is twofold: 1) there exists a unique optimal number of clusters\footnote{A formal proof is provided in Appendix~B.}; 2) once the optimal number of clusters is reached, further split is unnecessary and will only introduce duplicated cluster embeddings.

Based on this observation and given that we do not have the exact formulation of the optimal number, we will keep splitting, while always initializing the embedding of the newly split cluster to be the same as the chosen cluster (this will retain the optimality of the previous solution). Besides, in online learning settings where we continually receive streaming data, the optimal number of clusters may keep increasing. In this sense, cluster split is well-motivated and makes CEL a sustainable learning framework.

\subsection{CEL-Lite for Efficient Online Update}
\label{sec:online}

\subsubsection{Online Learning.}

Modern recommender systems generally deal with a massive amount of data.
Such systems are thus usually updated in an online fashion: data arrive in sequence over time, and updates take place when the accumulated data reaches a predefined quantity. This framework closely follows the online learning paradigm and allows the model to rapidly adapt to the newest data.

We adapt CEL into CEL-Lite for the online learning setting by optimizing its efficiency in dealing with batch data. CEL-Lite has a theoretically optimal time complexity and inherits most of the advantages of CEL. CEL-Lite is also suitable when we have access to the full data but want to train the model in an online fashion.

\subsubsection{Data and Embedding Optimization.}

Under the online setting, each new batch $\mathbf{X}_b$ consists of a set of user-item interactions $\{\mathbf{X}(i,j)\}$. To enable the usage of interaction history as well as mitigate \emph{catastrophic forgetting} \cite{forgetting}, we buffer from online batches a small number of the latest interactions ($\leq\!n$) from each user and item for \emph{experience replay} \cite{replay}, which is commonly used in online systems. Since each interaction is represented by a triplet of $\langle$user Id, item Id, rating$\rangle$, the storage cost incurred by buffering is acceptable. Moreover, the interactions can be stored in distributed file systems, unlike the embeddings or model parameters which must be stored in memory for efficient updates.

In CEL-Lite, one step of embedding optimization is performed for each batch via stochastic gradient descent with related buffered interactions. The overall time complexity of embedding optimization is then at most $\mathcal{O}(nDR)$, where $D$ denotes the total number of interactions $\|\mathbf{W}\|_0$.

\subsubsection{Efficient Cluster Reassignment.}

To achieve fast cluster reassignment, we incorporate the following accelerations into CEL-Lite: 1) only considers $m$ additional randomly sampled clusters for each reassignment; 2) only items involved in the current batch are considered for reassignment.

With these two strategies, the cluster assignments are only improved locally. Nevertheless, with sufficient updates, it may be able to converge to the global optimal. Accordingly to our experiments, such accelerations significantly improve efficiency, while only incurring a small drop in performance. The time complexity of each reassignment is now only $\mathcal{O}(nmbR)$, where $b$ is the batch size. Given there are $D/b$ batches, there are at most $D/b$ cluster reassignments. Thus, the overall time complexity is $\mathcal{O}(nmDR)$.

\subsubsection{Efficient Cluster Split.}

To control the overall complexity of cluster split for CEL-Lite, we propose to constrain the number of associated interactions in each cluster and ensure clusters are approximately balanced. These two conditions limit the maximal number of split operations and the time complexity of each split. Concretely, we set up the following split strategies:
\begin{strategy}
If an item has more than $d$ interactions, it will directly split and form a cluster by itself.
\label{strategy:1}
\end{strategy}
\begin{strategy}
A cluster will split if and only if it has more than $2d$ associated interactions ($\|\mathbf{W}\mathbf{S}_q(:,k)\|_0$). The split should be balanced such that the difference between the number of associated interactions of resulting clusters is within $d$.
\label{strategy:2}
\end{strategy}

Strategy~\ref{strategy:1} handles extremely warm items. After executing strategy~\ref{strategy:1}, each clustered item has no more than $d$ interactions; strategy~\ref{strategy:2} then gets to easily balance the splits, which can be done by tuning\footnote{The tuning of $\delta$ can be done by \emph{Quickselect} in $\mathcal{O}(\|\mathbf{S}_q(:,k)\|_0)$.} the threshold $\delta$ in Eq.~\eqref{equ:pca score}.

With strategy \ref{strategy:1} and \ref{strategy:2}, all clusters have at least $d/2$ interactions. Therefore, there are at most $2D/d$ clusters and hence at most $2D/d$ GPCA splits in total. Calculating the gradient of each item with only their $n$ buffered interactions, the time complexity of a single GPCA split is approximately $\mathcal{O}(n\|\mathbf{S}_q(:,k)\|_0R)$, which is no more than $\mathcal{O}(2ndR)$. The overall time complexity of splits is thus $\mathcal{O}(4nDR)$.

Note that we use a constant $d$ in strategies \ref{strategy:1} and \ref{strategy:2} for simplicity and ease of use in practice. It can in fact generalize to any monotonically nondecreasing sublinear functions of $D$ without altering the overall time complexity.
These strategies allow the maximum number of clusters to grow according to $D$, which practices our theoretical results in Section \ref{sec:identifiability}.

\subsubsection{The Optimal Time Complexity}
\label{sec:time complexity}

Regarding $n,m$ as constants, each of the three optimizations of CEL-Lite is of a time complexity of $\mathcal{O}(DR)$. Thus, the overall complexity of CEL-Lite is also $\mathcal{O}(DR)$.
Note that once we need to process all the data and each invokes an update of the embedding, the time complexity is already $\mathcal{O}(DR)$. That is, the time complexity of $\mathcal{O}(DR)$ in a sense is the lower bound of embedding learning.

\section{Experiments}

\begin{table}[t]
\centering
\caption{Datasets summary.}
\resizebox{.42\textwidth}{!}{
\begin{tabular}{lrrr}
\hline
Dataset        & Users & Items & No. interactions \\ \hline

MovieLens-1m   & $6,040$  & $3,952$  & $1,000,209$ \\

Video Games  & $826,767$ & $50,210$ & $1,324,753$ \\
Electronics  & $4,201,696$ & $476,002$ & $7,824,482$ \\
One-day Sales  & $\approx 13$m & $\approx40$m & $\approx700$m \\\hline
\end{tabular}}
\label{tab:datasets}
\end{table}

In this section, we empirically analyze CEL and compare it with related embedding learning baselines, including: Vanilla full embedding; eTree \cite{eTREE}; JNKM \cite{jnkm}; Modulo \cite{modulo} (i.e., with plain hash); BH \cite{binaryhash}; AE \cite{DeepRec}; HashGNN \cite{hashgnn}; and PreHash \cite{prehash}.
They are tested with a variety of feature interaction models, including: NMF, DeepFM~\cite{guo2017deepfm}, DIN~\cite{din}, and a business model from an E-commerce platform.

The datasets considered include: MovieLens, a review dataset with user ratings ($1\!\sim\!5$) on movies; Video Games and Electronics, review datasets with user ratings ($1\!\sim\!5$) of products from Amazon; and One-day Sales, sub-sampled real service interaction data within a day of  the E-commerce platform Shopee. The key statistics of these datasets are given in Table~\ref{tab:datasets}.

The original eTree and JNKM do not share embedding within clusters. To make the comparison thorough, we adapt eTree such that in the test phase, leaf nodes (users/items) will share the embedding of their parent nodes. Similarly, we adapt JNKM such that clustered items share the embedding of their cluster in the test phase. PreHash has a hyperparameter $K$ that selects top-$K$ related clusters for each user or item. We report its results with different values of $K$.

We reuse $M_q$ to also denote: the number of parent nodes in eTree; the number of clusters in JNKM; and the number of hash buckets in hash-based methods. We refer to $M_q/M$ as the \textbf{compression ratio}.

\begin{table*}[htb]
\centering
\caption{The test Mean Square Error (MSE) of rating predictions on MovieLens-1m dataset (averaged over 5-fold cross validation). A compression ratio (defined as $M_q/M$) of $100\%$ means the result is obtained under full embedding, while other compression ratios indicate that the results are obtained with clustered/hashed/shared embedding. Results marked with (p) indicate the results obtained with personalization (Section~\ref{sec:personalization}).}
\resizebox{0.81\textwidth}{!}{
\begin{tabular}{c|ccc|cccc|cccc}
\hline
\multicolumn{4}{c|}{
Feature Interaction Model
} & \multicolumn{4}{c|}{NMF} & \multicolumn{4}{c}{DIN} \\ \hline
\multicolumn{4}{c|}{Compression Ratio}
             &$100\%$ &$5\%$ &$1\%$  &$0.5\%$ &$100\%$ &$5\%$ &$1\%$   &$0.5\%$       \\ \hline
 & LTC & LMC & SIT & & & & & & & &
\\ \hline
Vanilla full embedding & \checkmark & \ding{53} & \checkmark & $0.8357$ & - & -  & -
                  & $0.7903$ & - & -  & - \\
eTree & \ding{53} & \ding{53} & \ding{53}        & $0.7841$ & $0.8399$ & $0.8482$ & $0.8614$
             & $0.7640$ & $0.8306$ & $0.8379$ & $0.8470$ \\
JNKM& \ding{53} & \ding{53} & \ding{53}        & $0.7626$ & $0.8903$ & $0.8616$ & $0.8498$
             & $0.7486$ & $0.8868$ & $0.8560$ & $0.8466$ \\
Modulo & \checkmark & \checkmark & \checkmark    & - & $1.0227$ & $1.0688$  & $1.0732$
             & - & $1.0206$ & $1.0647$ & $1.0692$ \\
BH     & \checkmark & \checkmark & \checkmark     & - & $0.9753$ & $1.0170$  & $1.0526$
             & - & $0.9907$ & $1.0040$ & $1.0495$ \\
AE    & \checkmark & \checkmark & \checkmark       & - & $0.9959$ & $1.0343$  &  $1.0630$
             & - & $0.9895$ & $1.0280$ & $1.0624$ \\
HashGNN   & \ding{53} & \ding{53} & \ding{53}   & - & $0.8491$ & $0.8829$ & $0.8882$
             & - & $0.8467$ & $0.8663$ & $0.8723$ \\
PreHash ($K=4$)   & \checkmark & \checkmark & \ding{53} & - & $0.8169$ & $0.8832$ & $0.8987$
                    & - & $0.8057$ & $0.8678$ & $0.8897$ \\
PreHash ($K=M_q$) & \ding{53} & \checkmark & \ding{53} & - & $0.7958$ & $0.8047$ & $0.8101$
                    & - & $0.7871$ & $0.7893$ & $0.7942$ \\
CEL (Ours)  & \ding{53} & \checkmark & \checkmark     & $\textbf{\small{\textsf{0.7507}}}^p$ & $\textbf{\small{\textsf{0.7784}}}$ & $\textbf{\small{\textsf{0.7858}}}$ & $\textbf{\small{\textsf{0.7926}}}$ & $\textbf{\small{\textsf{0.7390}}}^p$ & $\textbf{\small{\textsf{0.7718}}}$ & $\textbf{\small{\textsf{0.7787}}}$ & $\textbf{\small{\textsf{0.7868}}}$ \\
CEL-Lite (Ours) & \checkmark & \checkmark & \checkmark  & $\mathbf{0.7519}^p$ & $\mathbf{0.7820}$ & $\mathbf{0.7901}$ & $\mathbf{0.8014}$ & $\mathbf{0.7421}^p$ & $\mathbf{0.7766}$ & $\mathbf{0.7789}$ & $\mathbf{0.7906}$ \\ \hline
\end{tabular}}
\label{tab:movielens ranting prediction}
\end{table*}

\subsection{Performance Versus Baselines}
\subsubsection{Rating Prediction}

In this experiment, we perform rating prediction on the MovieLens datasets. We adopt NMF as well as DIN (both with $R=64$) as our feature interaction models, respectively. We set $t_1=40$ and $t_2=10$ for CEL, while $b=2000$, $n=20$, $m=10$, $t_1=1$ for CEL-Lite.
Since some baselines apply their methods only on items or users, to ensure a fair comparison, for this experiment we only compress (clustering, hashing, etc.) the item embeddings while keeping full user embeddings for all methods.

The results are shown in Table~\ref{tab:movielens ranting prediction}.
It can be observed that CEL outperforms existing baselines under all settings. According to our results, training with predefine clustering structure, eTree and JNKM brings clear improvements over vanilla methods in the full embedding setting. However, they consistently underperform CEL, especially when using a compressed embedding table. For hash based methods, the improvements brought by BH and AE (over naive Modulo) are relatively limited, while HashGNN and PreHash achieve more clear performance gains. Among them, PreHash with $K=M_q$ achieves the best performance in the compressed settings. Nevertheless, CEL clearly outperforms PreHash.
Meanwhile, the results of this experiment demonstrate that CEL-Lite, being much more computationally efficient than CEL, has a similar performance to CEL and still clearly outperforms other baselines. The degradation in performance might be attributed to its subsampling accelerations.

In addition, we check in Table~\ref{tab:movielens ranting prediction} three important properties for each method: Low Time Consumption (\textbf{LTC}), which indicates whether it achieves a time complexity of $\mathcal{O}(DR)$ for both training and test. In particular, algorithms whose time complexity scales linearly with $M_q$ is not regarded as being LTC. Thus, eTree, JNKM, PreHash (with $K=M_q$), and CEL is not regarded as LTC. Neither does HashGNN, as it requires running GNN aggregation of neighborhoods, which has a worst case time complexity exponential in $D$; Low Memory Consumption (\textbf{LMC}), which indicates whether the method can significantly reduce the memory cost for both training and testing compared with full embedding. JNKM, eTree, and HashGNN require a full embedding table during training, thus do not satisfy LMC; Support Incremental Training (\textbf{SIT}), which indicates whether the method supports incremental training, is crucial for practical systems. As the data increases, the predefined tree/clustering structure of eTree/JNKM can be sub-optimal, the pretrained hash function of HashGNN may become outdated, and the pre-specified pivot users in PreHash may become less representative. As a result, these methods may require periodical retraining.

\begin{table*}[ht]
\centering
\caption{The test AUC (presented in \% for a better interpretation) on Video Games and Electronics datasets.}

\resizebox{0.99\textwidth}{!}{
\begin{tabular}{c|cccc|cccc|cccc|cccc}
\hline
Dataset      & \multicolumn{8}{c|}{Video Games}  & \multicolumn{8}{c}{Electronics} \\ \hline
Feature Interaction      & \multicolumn{4}{c|}{DeepFM}  & \multicolumn{4}{c|}{DIN} & \multicolumn{4}{c|}{DeepFM}  & \multicolumn{4}{c}{DIN} \\ \hline
Compression Ratio     & $100\%$    & $10\%$    & $5\%$    & $2\%$   & $100\%$    & $10\%$    & $5\%$    & $2\%$   & $100\%$    & $10\%$    & $5\%$    & $2\%$   & $100\%$    & $10\%$    & $5\%$& $2\%$   \\ \hline
Vanilla full embedding      & $68.19$ & - & - & - & $70.73$ & - & - & -
            & $66.95$ & - & - & - & $67.62$ & - & - & -  \\
Modulo      & - & $61.68$ & $58.81$ & $56.34$
            & - & $63.35$ & $61.94$ & $60.25$
            & - & $63.13$ & $54.54$ & $48.61$
            & - & $63.70$ & $61.69$ & $58.93$ \\
HashGNN     & - & $68.51$ & $64.43$ & $61.84$
            & - & $71.59$ & $71.11$ & $70.82$
            & - & $64.49$ & $61.41$ & $59.45$
            & - & $65.50$ & $63.49$ & $60.40$ \\
PreHash ($K=M_q$)       & - & $69.19$ & $66.21$ & $63.89$
                        & - & $72.72$ & $71.63$ & $71.46$
                        & - & $67.50$ & $67.18$ & $66.99$
                        & - & $70.30$ & $68.80$ & $68.77$ \\
CEL-Lite    & $\mathbf{70.85}^p$ & $\mathbf{70.33}$ & $\mathbf{70.21}$ & $\mathbf{69.82}$
            & $\mathbf{73.85}^p$ & $\mathbf{73.80}$ & $\mathbf{73.50}$ & $\mathbf{73.39}$
            & $\mathbf{68.98}^p$ & $\mathbf{68.90}$ & $\mathbf{68.69}$ & $\mathbf{68.03}$
            & $\mathbf{71.63}^p$ & $\mathbf{71.56}$ & $\mathbf{70.18}$ & $\mathbf{69.81}$ \\ \hline
\end{tabular}}
\label{tab:cr predictions}
\end{table*}
\subsubsection{Online Conversion Prediction}

\label{exp:Conversion Rate}

To promote desirable behaviors from users, such as clicks and purchases, conversion rate prediction is an important learning task in online E-commerce systems. In this section, we test the models with the conversion prediction task, i.e., predicting whether a user will perform a desirable behavior on an item. We use the Video Games and Electronics datasets for this experiment, treating ratings not less than $4$ as positive feedback (i.e., conversions). For performance evaluation, we use the standard Area under the ROC Curve (AUC) on the test sets, which can be interpreted as the probability that a random positive example gets a score higher than a random negative example.

The experiment is conducted in an online fashion, such that the data is received incrementally over time. We adopt HashGNN and PreHash, the top-2 performing hash methods in the rating prediction experiment (as well as the naive Modulo) as our baselines, though they do not support incremental training. The hash function of HashGNN requires pretraining with full data before starting the actual embedding learning, thus we allow it to access the full data before the online training. PreHash requires binding the warmest users to its buckets, we thus allow it to go through all data to find the warmest users before the online training. It is worth noting that, unlike these two methods, CEL(-Lite) can naturally start embedding learning from scratch without any extra information.

In this experiment, we compress both the user embeddings and the item embeddings. To further demonstrate that our method works well with different feature interaction models, we adopt DeepFM in this experiment, while keeping using DIN. We set $R=8$ and $b=256$. We set $n=20$, $m=50$, and $t_1=1$ for CEL-Lite.

The results are shown in Table~\ref{tab:cr predictions}. It can be observed that CEL-Lite still clearly outperforms the baselines.

\begin{table}[ht]
\centering
\caption{The test AUC  (\%) on One-day Sales dataset.}
\resizebox{.35\textwidth}{!}{
\begin{tabular}{c|cccc}
\hline
Dataset                & \multicolumn{4}{c}{One-day Sales} \\ \hline
\begin{tabular}[c]{@{}c@{}}Embedding size \\
\end{tabular} & Full          & 1m       & 100k        & 10k         \\

\hline
Business Model                 & $87.09$      & $86.80$     & $86.55$          & $86.46$          \\
CEL-Lite                    & -           & -          & $\mathbf{87.78}$     & $\mathbf{87.70}$     \\ \hline
\end{tabular}}
\label{tab:one day sales}
\end{table}

\subsubsection{Integrated with Business Model}

To further validate the performance of CEL-Lite, we conduct experiments on dataset from the E-commerce platform Shopee, that needs to handle billions of requests per day. We use the model deployed in that platform as our baseline, which is a fine-tuned business model that has a sophisticated network structure and heavy feature engineering.
It optionally uses full embedding, or uses Modulo to compress the size of embedding table. In this experiments, we integrate CEL-Lite into this business model, replacing its embedding learning module.

As shown in Table~\ref{tab:one day sales}, CEL-Lite with a small number of clusters (with 10k or 100k embedding entities for users and items, respectively) outperforms the business model with full embedding, and more significantly, the business model with compression.
Specifically, CEL-Lite achieves an improvement of $+0.61\%$ in AUC with the embedding table size being 2650 times smaller.
Note that a $+0.1\%$ gain in absolute AUC is regarded as significant for the conversion rate tasks in practice \cite{cheng2016wide,din}.

\begin{figure}[h]
    \centering
    \includegraphics[width=.45\textwidth]{plots/hist.pdf}
    \caption{Averaged test MSE of items with different number of interactions. PreHash (with $K=M_q$) is the best-performing baseline in Table~\ref{tab:movielens ranting prediction}. CEL significantly outperforms PreHash on cold items. Personalization degenerates the performance of cold items, while improves the performance of warm items.}
    \label{fig:loss vs frequency}
\end{figure}

\begin{figure*}[ht]
    \centering
    \hspace{-1mm}
    \includegraphics[width=57mm]{plots/dense_ml10m_twolines.pdf}
    \includegraphics[width=56mm]{plots/normal_ml10m_twolines.pdf}
    \includegraphics[width=57mm]{plots/sparse_ml10m_twolines.pdf}
    \hspace{6mm}(a)\hspace{52mm} (b)\hspace{53mm} (c)

    \caption{The test MSE of different (initial) compression ratios on (a) the \emph{dense} MovieLens-1m, (b) the \emph{original} MovieLens-1m, and (c) the \emph{sparse} MovieLens-1m. Standard errors over $5$ runs are indicated by the shades. CEL with personalization at a given compression ratio indicates the fully personalized result, initiated with CEL at that compression ratio. $100\%$ means gradually split to full embedding.}
    \label{fig:sparsity vs ratio}
\end{figure*}

\subsection{Personalization and Cold-start Problems}
\subsubsection{Personalization.} \label{sec:personalization}

When training with shared embeddings, we may choose to untie the cluster assignments at a specific point to learn a personalized embedding (with the shared embedding as initialization) for each user and item. After being untied, embeddings can be further optimized according to their own historical data. This makes sense if we want a full embedding, in case we have enough data and memory. We refer to this procedure as \emph{personalization}, which is similar to the personalization in federated learning settings \cite{tan2022towards}. There are various personalization techniques. We leave them for future investigations. Here, we simply add the following distance regularization to the original objective to guide the personalization:
\[\mathcal{L}_{p}=\frac{\lambda_p}{2}\sum_{j=1}^M\|\mathbf{B}(j,:)-\mathbf{S}_q(j,:)\mathbf{B}_q\|_2^2.\]
We set $\lambda_p=50$ in our experiment.
Table~\ref{tab:movielens ranting prediction}\&\ref{tab:cr predictions} shows CEL with personalization achieves the best results. In particular, it has outperformed all other full embedding methods.

\subsubsection{Mitigating Cold-start Problems.}

We investigate the test loss of items with respect to the number of associated interactions in the training set (on the MovieLens-1m dataset, with a compression ratio of $1\%$). As shown in Figure~\ref{fig:loss vs frequency}, CEL significantly outperforms PreHash on cold items. This could be because PreHash emphasizes the warm items, while the cold items are roughly represented as a weighted sum of the warm ones. This may also be why PreHash performs slightly better in the interval of $[800,2000]$. On the other hand, we can see that if we do personalization for CEL, the performance of cold items degenerates. This strongly evidences those cold items cannot be well optimized with only their own data, and CEL can help in this regard. Personalization brings benefits for the intervals of $>\!20$, hence can have a better overall performance. We further try only personalizing the warm items ($>\!20$), and find a $1.2\%$ improvement on test MSE.

\subsubsection{Sparsity and Compression Ratio.}

To investigate the impact of data sparsity and compression ratio on CEL, we create a dense MovieLens-1m dataset by filtering out the users and items with no more than $100$ interactions, and a sparse MovieLens-1m dataset by randomly dropping $50\%$ interactions (but still ensuring at least $1$ sample) for each user.

As shown in Figure~\ref{fig:sparsity vs ratio}, on the dense and original datasets, increasing the number of clusters consistently improves the overall performance of CEL; while on the sparse dataset, the optimal compression ratio occurs at $2\%$: further splits will degenerate its performance.

We can also see that personalization at different compression ratios consistently brings improvement over CEL on the dense and original datasets, which is reasonable as warm items may have dominated the overall performance in these datasets. By contrast, on the sparse dataset, personalization consistently leads to degeneration of performance.
More specifically, on the dense dataset, increasing the number of clusters consistently improves the performance of personalization, while on the original and sparse dataset, there is an optimal compression ratio of $1\%\sim2\%$ for personalization.

Overall, the number of clusters is not the more the better. Too many clusters on a sparse dataset or a dataset that has many cold items, may lead to performance degeneration: partially on cold items or even the overall performance.

\subsection{Ablation Studies on Clustering}
\subsubsection{Interpretation of Clustering.}

To see how and how well CEL clusters items, we visualize the learned clusters of CEL on MovieLens-1m at $M_q=10$ in Figure~\ref{fig:cel_interpretation}.

We find that at this stage, the clusters are mainly divided according to the ratings of items. The standard deviations of ratings within clusters are generally low. This result is reasonable because the task is to predict the ratings.

We can also see that movies within a cluster appear to have some inner correlation. For example, cluster $6$ mainly consists of action or science fiction action movies; cluster $5$ mainly consists of music or animation movies that are suitable for children; cluster $9$ mainly consists of crime movies that have gangster elements; and cluster $8$ contains $70\%$ movies directed by Stanley Kubrick. Note that such information is not fed to the model in any form, but learned by CEL automatically.
A calculation of the averaged entropy of the distribution of genres in each cluster further evidences that CEL can better distinguish movies of different genres than others.

Interestingly, we find movies of the same type may form multiple clusters of different average ratings. For example, clusters $4$ and $6$ are both mostly action, but with a clear difference in their ratings.

\begin{table}[htb]
\centering
\caption{The test MSE of rating predictions on MovieLens-1m. We compare CEL with rule-based clusterings, and test using rule-based clusterings as initialization for CEL.}
\scalebox{0.92}{
\begin{tabular}{cc|ccc}
\hline
\multicolumn{2}{c|}{Number of clusters ($M_q$)} & $18$ & $40$ & $200$      \\ \hline
Initialization & Continue Running & && \\ \hline
By Genre &  - & $1.0515$ & - & - \\
By Avr. Rating &  -  & $0.8397$ & $0.8390$ & $0.8380$ \\
By Genre       & CEL & $0.8066$ & $0.7859$ & $\mathbf{0.7776}$ \\
By Avr. Rating & CEL & $\mathbf{0.7959}$ & $0.7892$ & $0.7833$ \\
None           & CEL & ${0.7985}$ & $\mathbf{0.7858}$ & $0.7784$ \\ \hline
\end{tabular}}
\label{tab:by_rule}
\end{table}

\subsubsection{Rule-based Clusterings}

CEL achieves automatic clustering of entities, but it is still not clear how it works against rule-based clusterings. Given that we find in the initial stage CEL mainly clusters items according to the genres and ratings of movies on the MovieLens dataset. In this section, we compare CEL with rule-based clustering based on genres (partitioning movies according to their genres in the meta data; there are 18 genres in total) and ratings (partitioning movies evenly according to their averaged rating).

As shown in Table~\ref{tab:by_rule}, neither of these two rule-based clusterings works well.
We further try using these two rule-based clusterings as initialization for CEL, to see whether it can help improve the performance. As we can see, such initialization brings either no or insignificant improvements.

\begin{figure*}[ht]
    \centering
    \includegraphics[width=1\linewidth]{plots/cel_clusters7.pdf}
    \caption{Clusters learned by CEL, ordered by their averaged ratings. We provide a one-line description for each cluster by summarizing the movies contained inside, and list 5 movies with the most number of ratings, with outliers marked by *. The table on the right shows the averaged entropy of the genre distribution in the clusters, calculated with the meta-data available.}
    \label{fig:cel_interpretation}
\end{figure*}

\subsubsection{Split Methods}

\begin{table}[htb]
\centering
\caption{Test MSE of different split methods.}
\scalebox{0.92}{
\begin{tabular}{ccc}
\hline
MovieLens         & 100k & 1m \\
\hline
Random Split      &  $0.9392$ & $0.8516$   \\
User-Inspired Bisecting 2-Means &  $0.9245$ & $0.8355$   \\
Gradient Bisecting 2-Means &  $0.8847$ & $0.8186$   \\
Gradient Random Projection &  $0.9201$  & $0.8269$   \\
GPCA              &  $\mathbf{0.8707}$ & $\mathbf{0.7858}$\\
\hline
\end{tabular}}
\label{tab:split methods}
\end{table}

We test several different split methods, including: Random Split, which randomly splits items into two sets; Bisecting 2-Means, an advanced split method from DHC \cite{steinbach2000comparison}, which bisects items into 2 clusters according to the distance between item representations; Random Projection, which is similar to GPCA but projects the gradients in a random direction instead of the first principal component; and GPCA.
Bisecting 2-Means is performed either on the gradient matrix or on the user-inspired representations (i.e., representing an item by averaging the embeddings of its associated users). For all methods, we split till a compression ratio of $1\%$.

The results are shown in Table~\ref{tab:split methods}. We can see that random split yields the worst result, which indicates the importance of split methods; Bisecting 2-Means works better on gradients than user-inspired representations, which suggests the gradients are more informative; Bisecting 2-Means and random project of gradient works less than GPCA, which suggests the rationality of split according to the first principal component.

\section{Conclusions}

In this paper, we proposed Clustered Embedding Learning (CEL), a novel embedding learning framework that can be plugged into any differentiable feature interaction model. It is capable of achieving improved performance with reduced memory costs. We have shown the theoretical identifiability of CEL and the optimal time complexity of CEL-Lite. Our extensive experiments suggested that CEL can be a powerful component for practical recommender systems.

We have mainly considered the embeddings of users and items in this paper. Nevertheless, CEL can be similarly applied to other ID features or categorical features, such as shop ID or seller ID in E-commerce, search words in search engines, object categories, and geographic locations. We leave these for future work.

\section*{Acknowledgments}
This research is supported in part by the National Natural Science Foundation of China (U22B2020); the National Research Foundation Singapore and DSO National Laboratories under the AI Singapore Programme (AISG2-RP-2020-019); the RIE 2020 Advanced Manufacturing and Engineering Programmatic Fund (A20G8b0102), Singapore; the Nanyang Assistant Professorship; the Shanghai University of Finance and Economics Research Fund (2021110066).

\bibliographystyle{ACM-Reference-Format}


\section{Introduction}

With the recent discovery of gravitational waves from binary black hole mergers \cite{LIGOScientific:2016aoc} and observation of the shadow of supermassive black hole candidates \cite{EventHorizonTelescope:2019dse,EventHorizonTelescope:2019ggy,EventHorizonTelescope:2020qrl,EventHorizonTelescope:2022wkp} we are entering an era in which proposed theories of gravity alternative to General Relativity (GR) may be tested via observations \cite{LIGOScientific:2016lio,Will:2014kxa,Yunes:2013dva}. The last few years have seen an enormous increase in the amount of theoretical works aimed at showing how to test the nature of black hole candidates and possible deviations from GR and the observations that are currently providing the most data are:
\begin{itemize}
    \item[(i)] Gravitational waves. So far, LIGO and Virgo have detected about one hundred gravitational waves from merger events of stellar mass black holes \cite{GWTC1,GWTC2,GWTC3,GWTC3-test}.
    \item[(ii)] Shadow. With  present technology the observation of the shadow is possible only for two supermassive black hole candidates, the one at the center of the galaxy M87 \cite{EventHorizonTelescope:2019dse,EventHorizonTelescope:2020qrl} and the one at the center of the Milky Way \cite{EventHorizonTelescope:2022wkp}.
    \item[(iii)] Accretion disk spectra from Active galactic Nuclei (AGNs). There are thousands of known spectra from quasars and AGNs, which bear information about the properties of the central objects.
    \item[(iv)] Accretion disk spectra from x-ray binaries. Stellar mass black hole candidates in the Milky Way have been detected also from the observation of their accretion disk spectra when they are in a binary system. Again the spectra provide information on the properties of the central object \cite{Bambi:2016sac,Cao:2017kdq,Bambi:2020jpe}.
    \item[(v)] Particle motion in the 'vicinity' of black hole candidates. The motion of stars orbiting a supermassive black hole may be described within the test particle approximation and such orbits have been observed only for the compact object at the Milky Way's center \cite{Ghez:1998ph,Ghez:2003qj,Ghez:2008ms,GRAVITY:2020gka}.
\end{itemize}
All of these observations may be used to determine some of the properties of black hole candidates. Under the assumption that the central object is a black hole as described in GR it may then be possible to use the observations to measure its mass and spin. On the other hand if the object is described by a black hole mimicker or a black hole in a theory different from GR the same observations may infer different values for the mass, spin and other relevant parameters of the object.
Many possibilities have been considered in the literature, from wormholes (see for example \cite{Li:2014coa,Nedkova:2013msa}) to modified black hole solutions (see \cite{Abdujabbarov:2016hnw,Cunha:2016wzk,Tsukamoto:2017fxq}), to boson stars and other black hole mimickers (see \cite{Cardoso:2019rvt,Vincent:2015xta,Abdikamalov:2019ztb,Li:2022eue}).

The general picture that is emerging is that one may not be able to rule out alternative theories of gravity and/or black hole mimickers from a single set of observations. Independent measurements of the same property, coming from independent observations are necessary. For this reason it is important to look at other potential observational tools to constrain the geometry in the vicinity of black hole candidates
\cite{Berti:2015itd, Cardoso:2019rvt}.

There have been several theoretical studies on the effects of curved spacetime on neutrino oscillations \cite{Wudka:1991tg,Grossman:1996eh,Cardall:1996cd,Piriz:1996mu,Fornengo:1996ef,Ahluwalia:1996ev,Bhattacharya:1999na,Pereira:2000kq,Crocker:2003cw,Lambiase:2005gt,Godunov:2009ce,Ren:2010yf,Geralico:2012zt,Chakraborty:2013ywa,Visinelli:2014xsa,Chakraborty:2015vla,Zhang:2016deq, Alexandre:2018crg, Blasone:2019jtj, Buoninfante:2019der, Boshkayev:2020igc, Mandal:2021dxk, Koutsoumbas:2019fkn, Swami:2020qdi, Capolupo:2020wlx}. A common conclusion reached in most of these works is the increased effective neutrino oscillation length due to gravitational redshift in curved spacetime backgrounds.
Inspired from this, one more experimental test that has recently been proposed to potentially break the aforementioned degeneracies from future measurements is the observation of gravitational lensing of neutrinos by massive sources and the effect that it has on neutrino oscillations \cite{Chakrabarty:2021bpr}. Of course, this kind of test is still out of reach from our present technological capabilities. Nevertheless it is important to explore the potential advantages that it may bear over the more 'traditional' lensing of photons.

It's a very intriguing proposal to test modified theories of gravity using neutrino experiments. We are particularly interested in theoretical models with extra dimensions that were introduced to solve the hierarchy problem between the electroweak and Planck scales as well as inspired by string theory \cite{Arkani-Hamed:1998jmv, Antoniadis:1998ig, Arkani-Hamed:1998sfv}. Indeed, the impact of extra dimensions on neutrino oscillation has been thoroughly investigated in the high energy physics community for over the past two decades. However, most of these investigations are based on the Large Extra Dimensions (LED) model and also make the assumption that sterile neutrinos are present and propagating in a higher-dimensional bulk \cite{Dienes:1998sb, Arkani-Hamed:1998wuz, Barbieri:2000mg} (see \cite{Basto-Gonzalez:2021aus} for recent studies). Here, we provide a different approach that is model independent and does not require the existence of sterile neutrinos to investigate the fingerprints of extra dimensions on the 4-dimensional curved spacetime.

In the present article we consider neutrino oscillations in the geometry of a Kaluza-Klein (KK) black hole \cite{Gibbons:1985ac,Rasheed:1995zv,Larsen:1999pp,Kudoh:2003ki,Horowitz:2011cq}.
Kaluza-Klein theory is a 5-dimensional extension of GR which includes Maxwell's electrodynamics by considering an additional field in the extra spatial dimension. By compactifying the extra dimension one obtains a 4-dimensional action for which the field equations reduce to those of Einstein and Maxwell.
Interestingly the KK black hole can also be obtained from a scalar-tensor theory of gravity coupled to Maxwell's electrodynamics with a dilaton field \cite{Frolov:1987rj,Larsen:1999pp}.
The shadow of a KK black hole was studied in \cite{Amarilla:2013sj,Mirzaev:2022xpz} while binary mergers and estimates for the spin and quasinormal modes were considered in \cite{Hirschmann:2017psw,Jai-akson:2017ldo}. More recently, the stability of orbits in the KK black hole was studied in \cite{Blaga:2023exi} while constraints on the theory's parameters from x-ray spectroscopy were obtained in \cite{Tripathi:2021rwb}.
Here we show that the probability of oscillation of neutrinos lensed by a KK black hole depends on the geometry via the value of the boost parameter $v$ and thus in principle may allow to distinguish such black holes from the ones in GR.

The article is organised as follows: A brief review of the KK black hole and neutrino oscillations are provided in sections \ref{KK} and \ref{neutrino}, respectively. In section \ref{4} we develop the formalism to describe neutrino oscillations in the KK metric and in section \ref{5} we obtain the oscillation probabilities for neutrinos lensed by stellar mass KK black holes.
The implications for future observations of astrophysical black hole candidates are discussed in section \ref{7}.
Throughout the article we make use natural units setting $G=c=1$, with the exception of the plots of the oscillation probabilities where we use GeV units in order to provide more realistic data as expected from lensing by a stellar mass black hole.

\section{The Kaluza-Klein metric}\label{KK}

We consider here the case of the Einstein-Maxwell-dilaton theory arising naturally as a low energy limit in string theory. The action of this theory is essentially that of gravity coupled to Maxwell and a dilaton field \cite{Horne:1992zy}
\begin{equation}
    S=\int d^4x \sqrt{-g}\left[-R+2(\nabla \varphi)^2 + e^{-2\alpha\varphi}F^2\right] .
\end{equation}
It has been shown that for $\alpha=\sqrt{3}$ the black hole solution obtained from the above action can be identified with a 5-dimensional KK black hole \cite{Frolov:1987rj,Larsen:1999pp}.
In fact, given a 4-dimensional vacuum solution of Einstein's equation in GR and taking its product with $\mathbb{R}$ one
obtains a 5-dimensional solution which is translation invariant in the extra dimension. Then the new solution is obtained by boosting, with boost parameter $v$, the 5-dimensional solution in the new direction. Therefore, to obtain a rotating black hole solution, one starts with the Kerr solution and applies the above procedure.
When reduced to 4D, the black hole solution obtained can be viewed as a `charged' solution with a non-trivial dilaton field \cite{Dobiasch:1981vh,Chodos:1980df,Gibbons:1985ac}. The resulting Kaluza-Klein line element in Boyer-Lindquist coordinates is given by \cite{Frolov:1987rj}
\begin{equation}\label{metric}
    \begin{aligned}
        ds^2 &= - \frac{1 - Z}{B} dt^2 - \frac{2 a Z \sin^2 \theta}{B \sqrt{1 - v^2}} dt d\phi + \\
        &+ \left[ B (r^2 + a^2) + a^2 \sin^2 \theta \frac{Z}{B} \right] \sin ^2 \theta d\phi^2 + \\
        &+ B \Sigma \left( \frac{dr^2}{\Delta} + d\theta^2 \right),
    \end{aligned}
\end{equation}
where
\begin{equation}
    \begin{aligned}
        B &= \left( 1 + \frac{v^2 Z}{1 - v^2} \right)^{1/2}, \\
        Z &= \frac{2mr}{\Sigma}, \\
        \Delta &= r^2 + a^2 - 2mr, \\
        \Sigma &= r^2 + a^2 \cos^2 \theta.
    \end{aligned}
\end{equation}
Here $ m $ and $ a $ are mass and spin parameters of the original Kerr solution and $ v $ is the boost velocity. The dilaton field is given by
\begin{equation}
    \varphi = -\frac{\sqrt{3}}{2}\log B.
\end{equation}
In terms of $m$, $v$, and $a$, the physical mass $M$, charge $Q$, and angular momentum $J$ as measured by observers at infinity are written as
\begin{eqnarray}\label{met-MJQ}
    M &=& m \left( 1 + \frac{v^2}{2(1 - v^2)} \right), \nonumber \\
    Q &=& \frac{mv}{1 - v^2}, \\
    J &=& \frac{ma}{\sqrt{1 - v^2}} \nonumber,
\end{eqnarray}
where the physical mass and charge are not independent anymore since they are related by the boost $v$ which is the result of the reduction of the boosted BH solutions in 5-dimensional gravity into 4-dimensional KK BHs. The locations of the outer and inner horizon of this KK black hole coincides with that of the Kerr solution.

\section{Neutrino Oscillations}\label{neutrino}

Neutrinos are produced and detected in different flavor eigenstates given by $\ket{\nu_\alpha}$, where $\alpha=e,\mu,\tau$. These flavor eigenstates are a superposition of three mass eigenstates denoted by $\ket{\nu_i}$,
\begin{equation}
     \ket{\nu_\alpha} = \sum_{i=1}^3 U^*_{\alpha i}\ket{\nu_i},
\end{equation}
where $i=1,2,3$ and $U^*_{\alpha i}$ is a $3\times3$ unitary mixing matrix \cite{Pontecorvo:1957qd,Pontecorvo:1967fh,Maki:1962mu}.
Let us assume that the neutrinos have a plane-wave wavefunction and it propagates from source $S$ at $(t_S,x_S)$ to detector $D$ at $(t_D,x_D)$ (see Figure \ref{fig1}). So the evolution of the wavefunction at the detector can be given by
\begin{equation}
    \ket{\nu_i(t_D,x_D)} = \exp(-i\Phi_i)\ket{\nu_i(t_S,x_S)},
\end{equation}
where $\Phi_i$ is the phase of oscillation. Neutrinos produced at the source $ S $ in a flavor eigenstate $ \ket{\nu_\alpha} $ travel to the detector location $ D $ and in this case, the probability of the neutrino changing flavor from the state $\nu_\alpha$ to $\nu_\beta$ is given by
\begin{equation}
    \begin{aligned}
        \mathcal{P}_{\alpha\beta} &= |\braket{\nu_\beta | \nu_\alpha(t_D,x_D)}|^2 = \\
        &= \sum_{i,j}U_{\beta i}U_{\beta j}^*U_{\alpha j}U_{\beta i}^* \exp\left[ -i(\Phi_i - \Phi_j  ) \right].
    \end{aligned}
\end{equation}
In flat spacetime, the phase $ \Phi_i $ has the usual form
\begin{equation}\label{phaseflat}
    \Phi_i = E_i(t_D - t_S) - p_i(\vec{x}_D-\vec{x}_S).
\end{equation}
In a typical scenario, we assume that the mass eigenstates in a flavour eigenstate initially produced at the source have equal momentum or energy \cite{Bilenky:1978nj}. This assumption, together with $ (t_D-t_S) \simeq |(\vec{x}_D-\vec{x}_S)| $ for relativistic neutrinos ($ E_i \gg m_i $), leads to the phase difference
\begin{equation}
    \Delta \Phi_{ij} \equiv \Phi_i - \Phi_j \simeq \frac{\Delta m_{ij}^2}{2E_0}(\vec{x}_D-\vec{x}_S),
\end{equation}
where $ \Delta m_{ij}^2 = m_i^2 - m_j^2 $, and $ E_0 $ is the average energy of the relativistic neutrinos produced at the source. Therefore, the phenomena of neutrino oscillations in vacuum depend only on the squared mass differences. However, in curved spacetime, the phase and hence the probability depends also on the sum of mass squared of individual mass eigenstates \cite{Swami:2020qdi}.

In curved spacetime, the expression of the phase $ \Phi_i $ of neutrino propagation can be generalized to a covariant form \cite{Fornengo:1996ef}
\begin{equation}
    \Phi_i = \int^D_S p_\mu^{(i)}dx^\mu,
\end{equation}
where
\begin{equation}
    p_\mu^{(i)} = m_i g_{\mu\nu}\frac{dx^\nu}{ds}
\end{equation}
is the canonical conjugate momentum to the coordinates $ x^\mu $ and $ g_{\mu\nu} $ and $ ds^2 $ are the metric tensor and line element of the curved spacetime, respectively.

\section{Neutrino Oscillation in KK metric}\label{4}

\begin{figure*}
    \begin{center}
        \includegraphics[width=14cm]{scheme.pdf}
    \end{center}
    \vspace{-0.5cm}
    \caption{ Schematic diagram for weak lensing of neutrinos in the KK black hole spacetime. Neutrinos propagate from the source $S$ to detector $D$ in the exterior of a static and non-rotating massive object described by the KK-metric. \label{fig1}}
\end{figure*}

In this section, we shall calculate the oscillation phases of neutrinos in the KK metric considering the simplest case of a non-rotating object, i.e. setting $ a=0 $. Here we shall follow the notations introduced in \cite{Fornengo:1996ef,Swami:2020qdi,Chakrabarty:2021bpr}. In the equatorial plane $\theta=\pi/2$, the metric components \eqref{metric} become
\begin{equation}\label{gmunu}
    \begin{aligned}
        g_{tt} &= -\left(1-\frac{2m}{r}\right)\left(1+\frac{2mv^2}{r(1-v^2)}\right)^{-1/2}, \\
        g_{rr} &= \left(1-\frac{2m}{r}\right)^{-1}\left(1+\frac{2mv^2}{r(1-v^2)}\right)^{1/2}, \\
        g_{\phi\phi} &= r^2\left(1+\frac{2mv^2}{r(1-v^2)}\right)^{1/2}.
    \end{aligned}
\end{equation}
We also define the components of canonical momenta $p_\mu^{(k)}$ where $k$ corresponds to the mass eigenstate,
\begin{equation}\label{can-mom}
    \begin{aligned}
        p_t^{(k)} &= m_k g_{tt}\frac{dt}{ds}, \\
        p_r^{(k)} &= m_k g_{rr}\frac{dr}{ds}, \\
        p_\phi^{(k)} &= m_k g_{\phi\phi}\frac{d\phi}{ds}.
    \end{aligned}
\end{equation}
The mass shell relation is given by
\begin{equation}
    -m^2 = g^{tt}p_t^2+g^{rr}p_r^2+g^{\phi\phi}p_\phi^2,
\end{equation}
where we have omitted the index $k$ in order to de-clutter the calculations.

\subsection{Radial case}

Let us first check the neutrino oscillation phase in the radial case, i.e. neutrinos are produced in the gravitational potential and travel outward radially. In this case, $d\phi=0$. Therefore, from Eq.\eqref{can-mom}, we have
\begin{equation}\label{dtds_drds}
    \frac{dt}{ds}=\frac{p_t}{m_k g_{tt}}, \ \ \ \ \ \ \frac{dr}{ds}=\frac{p_r}{m_k g_{rr}}.
\end{equation}
Since the spacetime is static, we know that the energy $E$ of test particles is conserved, and we define $ p_t = E_k $, $ p_0 = -E_0 $ and $p_r = p_k(r)$, where the subscript $k$ refers to the neutrino flavor. Now, the phase of a neutrino propagating radially in a null trajectory becomes
\begin{equation}\label{phik1}
    \Phi_k = \int_{r_S}^{r_D} \left[-E_k\left(\frac{dt}{dr}\right)_0+p_k(r)\right]dr,
\end{equation}
where the suffix $ 0 $ indicates the null trajectory. Now from Eq.\eqref{dtds_drds}
\begin{equation}\label{dtdr0}
    \left(\frac{dt}{dr}\right)_0 = -\frac{E_0}{p_0(r)}\frac{g_{rr}}{g_{tt}},
\end{equation}
on a null trajectory, where $E_0$ is the energy of the massless particle at infinity. The mass-shell relation now gives
\begin{equation}\label{p0pr}
    p_0 = \pm E_0\sqrt{-\frac{g_{rr}}{g_{tt}}}, \ \ \ \ \ p_k(r) = \pm \sqrt{-\frac{g_{rr}E_k^2}{g_{tt}}-m_k^2 g_{rr}}.
\end{equation}
Now using Eq.\eqref{dtdr0} and Eq.\eqref{p0pr}, we can write the phase in the null trajectory Eq.\eqref{phik1} as
\begin{equation}
    \Phi_k = \pm\int_{r_S}^{r_D}E_k \sqrt{-\frac{g_{rr}}{g_{tt}}}\left[\sqrt{1-\frac{m_k^2g_{tt}}{E_k^2}}-1\right]dr.
\end{equation}
Since $ 0 < g_{tt} < 1 $, we can expand the square root inside the bracket to get
\begin{equation}\label{phik2}
    \Phi_k = \pm \int_{r_S}^{r_D}\sqrt{-g_{tt}g_{rr}}E_k\frac{m_k^2}{2E_k^2}dr.
\end{equation}
Now, in relativistic approximation ($m_k<<E_k$), the following expressions hold
\begin{equation}
    \begin{aligned}
        &E_k \simeq E_0 +\mathcal{O}\left(\frac{m_k^2}{2E_0}\right), \\
        &E_k\frac{m_k^2}{2E_k^2} \simeq E_0 \frac{m_k^2}{2E_0^2}.
    \end{aligned}
\end{equation}
So using these two relations, Eq.\eqref{phik2} becomes
\begin{equation}
    \Phi_k = \pm \frac{m_k^2}{2E_0} \int_{r_S}^{r_D} \sqrt{-g_{tt}g_{rr}}dr.
\end{equation}
Now for the metric coefficients in Eq.\eqref{gmunu}, we have $ g_{tt}g_{rr} = -1 $. So the phase becomes
\begin{equation}
    \Phi_k = \pm \frac{m_k^2}{2E_0} (r_D - r_S),
\end{equation}
which has the same form as that in a flat spacetime. Hence, the oscillation probability of neutrinos travelling radially outward in the equatorial plane of a weak gravitational field will be independent of the gravitational effects.

\subsection{Non-radial case}

Let us now investigate what happens in the non-radial case. For a neutrino travelling along a null trajectory, the phase can be written using Eq.\eqref{can-mom} and Eq.\eqref{p0pr} as
\begin{equation}
    \Phi_k = \int_{r_S}^{r_D}\left[ -E_k \left( \frac{dt}{dr} \right)_0 + p_r + J_k \left( \frac{d\phi}{dr} \right)_0 \right]dr,
\end{equation}
where we use conservation of energy and angular momentum $ p_t^{(k)}=-E_k $ and $ p_\phi^{(k)} = J_k $ and $ J $ is the angular momentum of the test particle. Now, again from Eq.\eqref{can-mom}, we can write
\begin{equation}
    \frac{dt}{dr} = \frac{-E_k}{p_r}\frac{g_{rr}}{g_{tt}}, \ \ \ \ \ \ \frac{d\phi}{dr} = \frac{J_k}{p_r}\frac{g_{rr}}{g_{\phi\phi}}.
\end{equation}
Along the null trajectory, these equations take the form
\begin{equation}\label{nr-dtdr}
    \left(\frac{dt}{dr}\right)_0 = -\frac{E_0}{p_0(r)}\frac{g_{rr}}{g_{tt}}, \ \ \ \ \ \ \ \left(\frac{d\phi}{dr}\right)_0 = \frac{J_0}{p_0(r)}\frac{g_{rr}}{g_{\phi\phi}}.
\end{equation}
Now, we would like to express the angular momentum $J_k$ as a function of the energy $E_k$, the impact parameter $b$ and the velocity of the test particle at infinity $v_k^{(\infty)}$ as
\begin{equation}\label{nr-jk}
    J_k = E_k b v_k^{(\infty)}.
\end{equation}
The metric is asymptotically flat, so using relativistic approximation, we can express the velocity and the angular momentum as
\begin{equation}\label{nr-jkvk}
    \begin{aligned}
        &v_k^{\infty} = \frac{\sqrt{E_k^2-m_k^2}}{E_k} \simeq 1 - \frac{m_k^2}{2E_k^2}, \\
        &J_k  \simeq E_k b \left( 1 - \frac{m_k^2}{2E_k^2} \right),
    \end{aligned}
\end{equation}
and for a massless particle
\begin{equation}\label{nr-j0}
    J_0 = E_0 b.
\end{equation}
Using Eqs. \eqref{nr-dtdr}, \eqref{nr-jk}, \eqref{nr-jkvk} and \eqref{nr-j0}, we can write the phase for the non-radial case as
\begin{equation}
    \begin{aligned}
        \Phi_k = \int_{r_S}^{r^D}\frac{E_0 E_k g_{rr}}{p_0(r)} &\Bigg[ \frac{1}{g_{tt}} + \frac{p_0 p_k(r)}{E_0 E_k g_{rr}} + \\
        &+\frac{b^2}{g_{\phi\phi}}\left(1 - \frac{m_k^2}{2E_k^2}\right) \Bigg]dr.
    \end{aligned}
\end{equation}
Now, we shall simplify this expression using the mass-shell relation. The following two equations can be obtained from the mass-shell relation
\begin{equation}\label{nr-msrels}
    \begin{aligned}
        &\frac{p_0(r)}{E_0 g_{rr}} = \pm \sqrt{-\frac{1}{g_{tt}g_{rr}} - \frac{b^2}{g_{rr}g_{\phi\phi}}}, \\
        &\frac{p_0(r)p_k(r)}{E_0 E_k B} = -\left( \frac{1}{g_{tt}} + \frac{b^2}{g_{\phi\phi}} + \frac{m_k^2}{2E_k^2} \right),
    \end{aligned}
\end{equation}
And from these, we can express the phase as
\begin{equation}
    \Phi_k = -\int_{r_S}^{r_D} \frac{E_0 g_{rr}}{p_0(r)} E_k \frac{m_k^2}{2E_k^2}dr.
\end{equation}
In the relativistic approximation, we get
\begin{equation}
    \begin{aligned}
        \Phi_k &= \frac{m_k^2}{2E_0}\int_{r_S}^{r_D}\frac{E_0 g_{rr}}{p_0(r)}dr = \\
        &= \frac{m_k^2}{2E_0}\int_{r_S}^{r_D}\sqrt{-g_{tt}g_{rr}}\left(1-\frac{b^2|g_{tt}|}{g_{\phi\phi}}\right)^{-1/2}dr.
    \end{aligned}
\end{equation}
We can see that for $ b = 0 $, we recover the usual phase expression flat background spacetime. However, for $b\neq 0$ we see that the geometry now plays a role.

Now we shall evaluate the integral in two cases.
\begin{itemize}
    \item[(i)] Neutrinos are produced in the gravitational potential and travel outwards non-radially towards the detector.
    \item[(ii)] Neutrinos coming from a distant source are lensed by the gravitational potential before reaching the detector.
\end{itemize}
Let us consider the integrand in the first case. Replacing the metric coefficients, we get
\begin{equation}
    \begin{aligned}
        &\sqrt{-g_{tt}g_{rr}}\left( 1 - \frac{b^2|g_{tt}|}{g_{\phi\phi}} \right)^{-1/2} = \\
        = &\left[ \frac{b^2}{r^2}\left( \frac{2m}{v^2(2m-r)+r} - 1 \right) + 1 \right]^{-1/2} = \\
        = &\left[ \frac{b^2}{r^2}\left( \frac{4M}{v^2(4M-r)+2r} - 1 \right) + 1 \right]^{-1/2}.
    \end{aligned}
\end{equation}
Here, in the second equality, we have used Eq.\eqref{met-MJQ} to express the integrand in terms of the physical mass $M$. We now apply the weak field limit $ M/r << 1 $ to obtain
\begin{equation}
    \begin{aligned}
        &\left[ \frac{b^2}{r^2}\left( \frac{4M}{v^2(4M-r)+2r} - 1 \right) + 1 \right]^{-1/2} = \\
        = &\left( \frac{r^2}{r^2-b^2} \right)^{1/2} + \frac{2 b^2 M}{(v^2 - 2)(r^2 - b^2)^{3/2}} + \mathcal{O}(M^2).
    \end{aligned}
\end{equation}
Now we perform the integral to obtain
\begin{equation}
    \begin{aligned}
        &\int_{r_S}^{r_D}\sqrt{-g_{tt}g_{rr}}\left(1-\frac{b^2|g_{tt}|}{g_{\phi\phi}}\right)^{-1/2}dr = \\
        = &\sqrt{r^2-b^2} + \frac{Mr}{\sqrt{r^2-b^2}} -\frac{M r v^2}{(v^2 - 2)\sqrt{r^2-b^2}}.
    \end{aligned}
\end{equation}
The first two terms here are the Schwarzschild contribution and the third term arises because of the non-zero boost velocity parameter. This term can also be written in terms of the charge as
\begin{equation}
    -\frac{M r v^2}{(v^2 - 2)\sqrt{r^2-b^2}} = \frac{Q r v}{2\sqrt{r^2 - b^2}}.
\end{equation}
Now by solving the relation $ Q/M = 2v/(2-v^2) $, we can find that the parameter $ v $ has two real solutions
\begin{equation}
    v = \frac{\pm \sqrt{M^2 + 2Q^2} - M}{Q}.
\end{equation}
For small physical charge, $ Q/M << 1 $, the positive part yields $ |v| \sim \left| Q/M \right| << 1 $. On the other hand, the negative part yields $ |v| \sim 2\left| M/Q \right| >> 1 $, which is nonphysical. Going forward, we shall only consider the positive part of the solution. This gives
\begin{equation}
    \begin{aligned}
        &\int_{r_S}^{r_D}\sqrt{-g_{tt}g_{rr}}\left(1-\frac{b^2|g_{tt}|}{g_{\phi\phi}}\right)^{-1/2}dr = \\
        = &\sqrt{r^2-b^2} + \frac{Mr}{\sqrt{r^2-b^2}} + \frac{Q^2 r}{2M\sqrt{r^2 - b^2}} = \\
        = &\sqrt{r^2-b^2} + M\left( 1 + \frac{Q^2}{2M^2} \right)\frac{r}{\sqrt{r^2 - b^2}}.
    \end{aligned}
\end{equation}
Then, the phase can be written as
\begin{equation}
    \begin{aligned}
        \Phi_k &= \frac{m_k^2}{2E_0}\Bigg[ \sqrt{r_D^2 - b^2} - \sqrt{r_S^2 - b^2} + \\
        &+ M \left( 1 + \frac{Q^2}{2M^2} \right) \left( \frac{r_D}{\sqrt{r_D^2 - b^2}} - \frac{r_S}{\sqrt{r_S^2 - b^2}} \right)  \Bigg].
    \end{aligned}
\end{equation}
In the limit $ v \rightarrow 0 $, the above expression reduces to that of the Schwarzschild case \cite{Fornengo:1996ef,Swami:2020qdi}.

Now, let us consider the second case. Here, neutrinos are produced from a distant source, travel towards the lensing object, passing at closest approach at $ r = r_0 $ and subsequently reaching the detector. The phase, in this case, can be written in two parts taking into account the sign of the momentum and the symmetry of the trajectory about $r_0$. Then
\begin{equation}
    \begin{aligned}
        \Phi_{k}(r_S \rightarrow r_0 \rightarrow r_D) &= \frac{m_k^2}{2E_0}\int_{r_0}^{r_S}\sqrt{\frac{-g_{tt}g_{rr}}{1-\frac{b^2|g_{tt}|}{g_{\phi\phi}}}}dr + \\
        &+ \frac{m_k^2}{2E_0}\int_{r_0}^{r_D}\sqrt{\frac{-g_{tt}g_{rr}}{1-\frac{b^2|g_{tt}|}{g_{\phi\phi}}}}dr.
    \end{aligned}
\end{equation}
In the weak-field approximation, the point of closest approach can be determined by solving the equation
\begin{equation}
    \left( \frac{dr}{d\phi} \right)_0 = \frac{p_0(r_0)g_{\phi\phi}}{J_0 g_{rr}} = 0.
\end{equation}
The solution of this equation takes a simple form and is given by
\begin{equation}
    r_0 \simeq b - \frac{2M}{2-v^2} \simeq b - M\left( 1 + \frac{Q^2}{2M^2} \right).
\end{equation}
In the second equality, we considered a small charge $ v \sim Q/M << 1 $. Then we express the integrand of the phase in terms of $ r_0 $ in the weak field limit as
\begin{equation}
    \begin{aligned}
        &\sqrt{\frac{-g_{tt}g_{rr}}{1-b^2|g_{tt}|/g_{\phi\phi}}} = \frac{r}{\sqrt{r^2 - r_0^2}} + \\
        &- \frac{2 M r_0}{(v^2 - 2)(r+r_0)\sqrt{r^2 - r_0^2}} + \mathcal{O}(M^2/r^2).
    \end{aligned}
\end{equation}
Then we perform the integral to get
\begin{equation}
    \begin{aligned}
        &\int_{r_S}^{r_D}\sqrt{-g_{tt}g_{rr}}\left(1-\frac{b^2|g_{tt}|}{g_{\phi\phi}}\right)^{-1/2}dr = \\
        &= \sqrt{r^2-r_0^2} - \frac{2M\sqrt{r^2-r_0^2}}{(r+r_0)(v^2-2)} = \\
        &= \sqrt{r^2-r_0^2} + M\left( 1 + \frac{Q^2}{2M^2} \sqrt{\frac{r-r_0}{r+r_0}} \right),
    \end{aligned}
\end{equation}
where in the second equality, we have considered the small charge approximation $ v \sim Q/M << 1 $.
Finally, we obtain
\begin{equation}
    \begin{aligned}
        &\Phi_k (r_S \rightarrow r_0 \rightarrow r_D) = \frac{m_k^2}{2E_0}\Bigg[ \sqrt{r_D^2 - r_0^2} +\sqrt{r_S^2 - r_0^2} + \\
        &+ M\left( 1 + \frac{Q^2}{2M^2} \right) \left( \sqrt{\frac{r_D - r_0}{r_D + r_0}} + \sqrt{\frac{r_S - r_0}{r_S + r_0}} \right) \Bigg], \\
    \end{aligned}
\end{equation}
Or
\begin{equation}
    \begin{aligned}
        &\Phi_k = \frac{m_k^2}{2E_0}\Bigg[ \sqrt{r_D^2 - b^2} +\sqrt{r_S^2 - b^2} + \\
        &+ M\left( 1 + \frac{Q^2}{2M^2} \right) \Bigg( \frac{b}{\sqrt{r_D^2 - b^2}} + \\
        &+ \frac{b}{\sqrt{r_S^2 - b^2}} +\sqrt{\frac{r_D-b}{r_D + b}} + \sqrt{\frac{r_S-b}{r_S + b}}\Bigg)\Bigg].
    \end{aligned}
\end{equation}
Now, for $ b << r_{S,D} $, we expand the above expression over $ b/r_{S,D} $ and keep terms up to $ \mathcal{O}\left(b^2/r_{S,D}^2\right) $ to get
\begin{equation}
    \begin{aligned}
        \Phi_k &= \frac{m_k^2}{2E_0}(r_D + r_S) \Bigg[ 1 - \frac{b^2}{2 r_D r_S} + \\
        &+ \left( 1 + \frac{Q^2}{2M^2} \right)\frac{2M}{r_D + r_S} \Bigg].
    \end{aligned}
\end{equation}
So, we can see that the phase depends on the boost parameter through the charge $ Q $. When $ v \rightarrow 0 \ \ (Q \rightarrow 0) $ we recover the expression for the phase in the Schwarzschild spacetime \cite{Fornengo:1996ef}.

\begin{figure}
    \begin{center}
        \includegraphics[width=8.5cm]{pemu_h_v00.png}
    \end{center}
    \caption{ Neutrino oscillation probability for $ v = 0 $, i.e. in the Schwarzschild case, when the mixing angle is $ \pi/6 $. The solid and the dashed curves represent normal hierarchy and inverted hierarchy respectively. Values of the other parameters are as follows: $ M = 4M_\odot $, $ \left| \Delta m^2 \right| = 10^{-3} \rm{eV}^2 $, and the lightest neutrino is massless. \label{fig2}}
\end{figure}

\section{Neutrino oscillation probability}\label{5}

\begin{figure*}
    \begin{center}
        \includegraphics[width=8.5cm]{pemu_h_v02.png}
        \includegraphics[width=8.5cm]{pemu_h_v04.png}
    \end{center}
    \caption{ Left Panel: Neutrino oscillation probability for $ v = 0.2 $ when the mixing angle is $ \pi/6 $. Right panel: Neutrino oscillation probability for $ v = 0.4 $ when the mixing angle is $ \pi/6 $. The solid and the dashed curves represent normal hierarchy and inverted hierarchy respectively. Values of the other parameters are as follows: $ M = 4M_\odot $, $ \left| \Delta m^2 \right| = 10^{-3} \rm{eV}^2 $, and the lightest neutrino is massless.\label{fig3}}
\end{figure*}

\begin{figure*}
    \begin{center}
        \includegraphics[width=8.5cm]{pemu_v_nh.png}
        \includegraphics[width=8.5cm]{pemu_v_ih.png}
    \end{center}
    \caption{ Left Panel: Neutrino oscillation probability for $ v = 0 $ (solid curve), $ v = 0.2 $ (dashed curve) and $ v = 0.4 $ (dotted-dashed curve) for normal hierarchy. Right panel: Neutrino oscillation probability for $ v = 0 $ (solid curve), $ v = 0.2 $ (dashed curve) and $ v = 0.4 $ (dotted-dashed curve) for inverted hierarchy.

    Values of the other parameters are as follows: $ M = 4M_\odot $, $ \left| \Delta m^2 \right| = 10^{-3} \rm{eV}^2 $, $ \alpha = \pi/6 $ and the lightest neutrino is massless. \label{fig4}}
\end{figure*}

\begin{figure*}
    \begin{center}
        \includegraphics[width=8.5cm]{pemu_M_v_nh.png}
        \includegraphics[width=8.5cm]{pemu_M_v_ih.png}
    \end{center}
    \caption{ Left Panel: Neutrino oscillation probability for $ v = 0, M = M_\odot $ (solid curve), $ v = 0.4, M = M_\odot $ (dashed curve), $ v = 0, M = 4M_\odot $ (dotted-dashed curve) and $ v = 0.4, M = 4M_\odot $ (dotted curve) for normal hierarchy. Right panel: Neutrino oscillation probability for $ v = 0, M = M_\odot $ (solid curve), $ v = 0.4, M = M_\odot $ (dashed curve), $ v = 0, M = 4M_\odot $ (dotted-dashed curve) and $ v = 0.4, M = 4M_\odot $ (dotted curve) for inverted hierarchy.

    Values of the other parameters are as follows: $ \left| \Delta m^2 \right| = 10^{-3} \rm{eV}^2 $, $ \alpha = \pi/6 $ and the lightest neutrino is massless.  \label{fig5}}
\end{figure*}

In this section, we shall calculate the oscillation probability of neutrinos that are lensed by a KK black hole. We shall consider neutrinos with mass eigenstates $ \ket{\nu_i} $ travelling in the KK black hole geometry from a source $ S $ to a detector $ D $ through different classical paths. A flavour eigenstate $ \ket{\nu_\alpha} $ propagated from $ S $ to $ D $ through a path denoted by $ p $ is given by
\begin{equation}
    \ket{\nu_\alpha(t_D,x_D)} = N \sum_i U_{\alpha i}^* \sum_p \exp \left( -i\Phi_i^p \ket{\nu_i (t_S,x_S)} \right),
\end{equation}
where $ \Phi_i^p $ is the phase of neutrino oscillation with the impact parameter $ b_p $ now being a path-dependent parameter and $ N $ is a normalization factor. We suppose that a neutrino is produced in a flavor eigenstate $ \alpha $ at the source $ S $ and detected in a flavor eigenstate $ \beta $ at the detector $ D $. Then the probability of oscillation is given by
\begin{equation}\label{prob-1}
    \begin{aligned}
        \mathcal{P}_{\alpha\beta}^{lens} &= |\braket{\nu_\beta |\nu_\alpha(t_D,x_D)}|^2 = \\
        &= |N|^2 \sum_{i,j}U_{\beta i}U_{\beta j}^*U_{\alpha j}U_{\alpha i}^* \sum_{p,q} \exp \left(
        -i\Delta \Phi_{ij}^{pq} \right),
    \end{aligned}
\end{equation}
where $N$ is given by
\begin{equation}
    |N|^2 = \left( \sum_i |U_{\alpha i}|^2 \sum_{p,q} \left( -i\Delta \Phi_{ii}^{pq} \right) \right)^{-1}.
\end{equation}
$ \Delta \Phi_{ij}^{pq} $ is the phase difference between the paths $p$ and $q$ which is given by
\begin{equation}\label{phiAB}
    \Delta\Phi_{ij}^{pq} = \Phi_p^i - \Phi_q^j = \Delta m_{ij}^2 A_{pq} + \Delta b_{pq}^2 B_{ij},
\end{equation}
with
\begin{equation}\label{aijbpq}
    \begin{aligned}
        &A_{pq} = \frac{r_S + r_D}{2E_0}\left[ 1 + \left( 1 + \frac{Q^2}{2M^2} \right)\frac{2M}{r_S+r_D} - \frac{\sum b_{pq}^2}{4r_Sr_D} \right], \\
        &B_{ij} = - \frac{\sum m_{ij}^2}{8E_0}\left( \frac{1}{r_S} + \frac{1}{r_D} \right).
    \end{aligned}
\end{equation}
In the above set of equations, the quantities $ \Delta b_{pq}^2 $, $ \sum b_{pq}^2 $, $ \Delta m_{ij}^2 $ and $ \sum m_{ij}^2 $ are given by
\begin{equation}
    \begin{aligned}
        &\Delta b_{pq}^2 = b_p^2 - b_q^2, \\
        &\sum b_{pq}^2 = b_p^2 + b_q^2, \\
        &\Delta m_{ij}^2 = m_i^2 - m_j^2, \\
        &\sum m_{ij}^2 = m_i^2 + m_j^2.
    \end{aligned}
\end{equation}
We can clearly notice from the above expressions that $ A_{pq} $ and $ B_{ij} $ are invariant under the change of their respective indices, i.e. $ A_{pq} = A_{qp} $ and $ B_{ij} = B_{ji} $. This suggests that under a simultaneous change of both of its indices, the phase difference changes sign $ \Delta\Phi_{pq}^{ij} = -\Delta\Phi_{qp}^{ji} $. For those neutrino paths for which $ p = q $, we have that $ \Delta b_{pq}^2 $ vanishes and the oscillation probability only depends on the difference of squared masses of neutrinos. However, for general paths, $ p \neq q $, the oscillation probability also depends on the sum of squared masses. Therefore, in principle, gravitationally lensed neutrinos can provide information about the individual neutrino masses if the spacetime properties are known from other observations \cite{Swami:2020qdi}.

Now we would like to consider a simple toy model of two neutrino flavors to understand how the boost velocity $ v $ affects the oscillation probability. Specifically, we shall evaluate
the two-flavor neutrino oscillation probability at a generic point in a plane connecting the source $ S $, the lens and the detector $ D $ in the weak field limit. We substitute Eq.\eqref{phiAB} in Eq.\eqref{prob-1} to get
\begin{equation}
    \begin{aligned}
        \mathcal{P}_{\alpha\beta}^{lens} &= |N|^2\Bigg[ \sum_{i,j}U_{\beta i}U_{\beta j}^*U_{\alpha j}U_{\alpha i}^*\Big( \sum_{p=q} \exp \left( -i\Delta m_{ij}^2 A_{pp} \right) + \\
        & + 2\sum_{p>q} \cos \left( \Delta b_{pq}^2 B_{ij} \right) \exp \left( -i\Delta m_{ij}^2A_{pq} \right) \Big) \Bigg],
    \end{aligned}
\end{equation}
where
\begin{equation}
    |N|^2 = \left( N_{path} + \sum_{i}|U_{\alpha i}|^2 \sum_{q>p} 2 \cos \left( \Delta b_{pq}^2 B_{ii} \right) \right)^{-1}.
\end{equation}
For simplicity, we consider the neutrinos to be propagating in the equatorial plane ($ \theta = \pi/2 $) and, in this case, $ N_{path} = 2 $. To understand the probabilities qualitatively and quantitatively, we consider a transition from electron to muon neutrino, $ \nu_e \rightarrow \nu_\mu $. For this process, the probability can be written as
\begin{equation}
    \begin{aligned}
        \mathcal{P}_{e\mu}^{lens} &= |N|^2 \sin^2 2\alpha \Bigg[ \sin^2\left( \Delta m^2 \frac{A_{11}}{2}\right) + \sin^2\left( \Delta m^2 \frac{A_{22}}{2}\right) + \\
        &- \cos \left( \Delta b^2 B_{12} \right)\cos \left( \Delta m^2 A_{12} \right) + \frac{1}{2} \cos \left( \Delta b^2 B_{11} \right) + \\
        &+ \frac{1}{2}\cos \left( \Delta b^2 B_{22} \right)   \Bigg],
    \end{aligned}
\end{equation}
where
\begin{equation}
    |N|^2 = \left[ 2\left( 1 + \cos^2\alpha \cos(\Delta b^2 B_{11}) + \sin^2 \alpha \cos (\Delta b^2 B_{22}) \right) \right]^{-1}.
\end{equation}
In these expressions, $ \Delta b^2 = \Delta b_{12}^2 $, $ \Delta m^2 = \Delta m_{21}^2 $ and $ A_{pq} $ and $ B_{ij} $ are given by Eq.\eqref{aijbpq}. We would like to express the path-dependent impact parameter $b_p$ in terms of the geometric quantities of the system. Figure \ref{fig1} shows the schematic representation of the geometric system of lensing phenomena. Neutrinos are produced at the source $ S $, get lensed by a KK black hole metric and are later detected at $ D $. The physical distances from the source to the lens and from the lens to the detector are $ r_S $ and $ r_D $, respectively in the $ (x,y) $ coordinate system. Let us consider another coordinate system $ (x',y') $ obtained by rotating the $ (x,y) $ system by an angle $ \varphi $. These two systems then would be related by the relations $ x' = x \cos \varphi + y \sin \varphi  $ and $ y' = -x \sin \varphi + y \cos \varphi $. Now, the angle of deflection $\delta$ of neutrinos in the rotated frame can be given by
\begin{equation}
    \begin{aligned}
        \delta = \frac{y_D' - b}{x_D'} &= -\frac{4m}{b}\left( 1 + \frac{v^2}{2(1-v^2)} \right) = \\
        &= -\frac{4M}{b} = -\frac{2R_x}{b},
    \end{aligned}
\end{equation}
where $ R_x $ is just $ 2 M $ and $ M $ is the physical mass of the source. Here $ (x_D', y_D') $ is the location of the detector. Now from Figure \ref{fig1}, we use the identity $ \sin \varphi = b/r_S $ and the above equation becomes
\begin{equation}\label{polynom}
    \left( 2R_x x_D + b y_D \right)\sqrt{1-\frac{b^2}{r_S^2}} = b^2\left( \frac{x_D}{r_S} + 1  \right) - \frac{2R_x b y_D}{r_S}.
\end{equation}
Assuming a nearly-circular trajectory of the detector around the lens, we have $ x_D = r_D \cos\phi $ and $ y_D = r_D \sin \phi $. Now, Eq.\eqref{polynom} is a quartic polynomial and it can be solved to obtain the impact parameters in terms of $ r_S $, $ R_x $ and the lensing location ($ x_D, y_D $). For a quantitative understanding, we consider a Sun-Earth system with the Sun as the lens and the Earth as a detector. In this example, we shall use the typical values of the geometrical quantities in the solar-system, assuming that the gravitational field of the Sun is represented by the exterior of the KK metric in the weak field limit. We assume the source is behind the Sun and it emits high-energy neutrinos with $ E_0 = 10 \rm{MeV} $. We consider $ r_D = 10^8 \rm{km} $ and $ r_S = 10^5r_D $ and numerically solve Eq.\eqref{polynom} and obtain two real roots $b_1$ and $b_2$ for each $ \phi $. Physical mass of the lens is considered to be $ M = 1-5 M_\odot $ and neutrino mass squared difference is $ |\Delta m^2| = 10^{-3} \rm{eV}^2 $.

\begin{figure*}
    \begin{center}
        \includegraphics[width=8cm]{pemu_degen_M_v-15.pdf}
        \hspace{0.5cm}
        \includegraphics[width=8cm]{pemu_degen_M_v-22.pdf}
    \end{center}
    \caption{The degeneracy between the determination of the mass parameter $M$ and the boost parameter $v$ for a given probability $\mathcal{P}_{e\mu}$ at a fixed angle $\phi$ is illustrated by the 2D contour plot of the implicit function $M(v)$ obtained from $\mathcal{P}_{e\mu} (M,v)=\rm{const}.$, namely the plot of the values of the boost parameter $ v $ and the lens mass $ M $ which give the same probability. Each curve corresponds to a fixed value of $\mathcal{P}_{e\mu}$ as shown on the line itself. The left panel shows the contour plots for $\phi=0.0015$, while the right panel shows the contour plots $M(v)$ for $\phi=0.0022$. Keep in mind that the probability also changes with the angle, therefore for the same source the curves $M(v)$ would correspond to different $\mathcal{P}_{e\mu}$ at different angles. \label{fig6}}
\end{figure*}

The oscillation probability for the two flavor toy model of neutrinos is shown in Figures \ref{fig2}, \ref{fig3}, \ref{fig4} and \ref{fig5} as a function of the the azimuth angle $\phi$ that signifies the detector's position in the orbit. Figure \ref{fig2} shows the oscillation probability of gravitationally lensed neutrinos in a Schwarzschild spacetime, ($ v = 0 $). The solid and the dashed line corresponds to normal hierarchy ($ \Delta m^2 > 0 $) and inverted hierarchy ($ \Delta m^2 < 0 $) of neutrinos respectively. Verifying the results of \cite{Swami:2020qdi,Chakrabarty:2021bpr}, we can see that the probability in the case of inverted hierarchy is quite large for particular values of $\phi$.

We see similar behavior also in case of non-zero $v$, i.e. in a KK black hole spacetime, in Figure \ref{fig3} and \ref{fig4}.
From Figure \ref{fig4} we can see that for sufficiently large values of the boost parameter $ v $, the probability deviates significantly from that of the Schwarzschild case, suggesting that the detection of a sufficiently large number of neutrinos could provide a way to test the geometry. The solid line there shows the probability of oscillation for $v=0$ (Schwarzschild case), and the dashed and the dotted-dashed line shows the oscillation probability for $ v = 0.2 $ and $ v = 0.4 $ respectively. This is true for both normal and inverted hierarchy.

Next, we investigate the effect of the mass of the lens and the boost parameter together in Figure \ref{fig5}. The left and right panels of Figure \ref{fig5} show the probability of oscillation of neutrinos for normal ($\Delta m^2 > 0$) and inverted hierarchy ($\Delta m^2 < 0$) for two different masses of the lens. The solid and dashed line corresponds to a solar mass lens with $ v = 0 $ and $ v = 0.4 $, respectively. The dotted and dotted-dashed line corresponds to a $ 4M_\odot $ lens with $
v = 0 $ and $ v = 0.4 $, respectively. We can see that a more massive lens can create a bigger shift in the oscillation probability. In certain cases, the shift in probability for a $ 4M_\odot $ gravitational lens is almost double that of a $ M_\odot $ lens.

Finally we investigated the degeneracy for the measurements of the lens mass and boost parameter simultaneously. In Figure \ref{fig6} we plotted the implicit function $M(v)$ obtained at a given angle $\phi$ for a fixed value of the probability $\mathcal{P}_{e\mu}$. We see that for any given value of $\mathcal{P}_{e\mu}$ there are several values of $M$ and $v$ that produce the same probability. This shows that, similarly to what happens with other kinds of observations, one needs more measurements in order to break the degeneracy and determine $M$ and $v$. Nevertheless, it is worth mentioning that the implicit plot changes with $\phi$, which is the detector's azimuthal angle, which in turn, for Earth, depends on the time of the year. Therefore by tracking the implict plots for $M(v)$ at different angles which correspond to the probability $\mathcal{P}_{e\mu}$ at that angle one may be able to constrain the range of allowed values of $M$ and $v$ for the source. In principle, by overlapping the obtained range with another range of possible values of $M$ for the same source, as obtained through a different method (for example, x-ray spectroscopy), one could be able to further constrain $M$ and therefore $v$, thus testing the nature of the geometry.

\section{Discussion}\label{7}

We considered gravitational lensing of neutrinos by a non rotating KK black hole as a possible tool to estimate the properties of the geometry and potentially distinguish such hypothetical objects from black holes in classical GR. We showed how the boost parameter $v$ of the KK black hole affects the probability of neutrino oscillations and showed that there exist a degeneracy between in the measurement of the boost parameter and gravitational mass of the object using this method. As a consequence, the observation of a black hole candidate with this method would provide only a range for the allowed values of the object's mass and boost parameter and would not be able to rule out the validity of the KK metric for the source.
It is well known that similar degeneracies arise for the measurement of the properties of black hole candidates also from other methods such as gravitational waves, spectroscopy or shadow (see for example \cite{Joshi:2013dva,Bambi:2013hza,Abdikamalov:2019ztb} for specific examples or \cite{Berti:2015itd, Cardoso:2016ryw} for a more general discussion of tests of GR). We argued that an independent measurement of the object's mass with two separate methods would allow to reduce the allowed range of values for $v$ and thus constrain the validity of the KK metric as the exterior field of a compact object.

As of now almost all observations of any black hole candidate have been done using only one tool and the precision of the data is in most cases not sufficient to provide stringent constraints on the metric parameters.
Similarly, the detection of neutrinos from galactic and extra-galactic astrophysical sources is still in its infancy \cite{IceCube:2018cha,IceCube:2018dnn} and therefore the methods discussed here are not experimentally applicable today. However, we believe that in the future, when more data and increased precision will be available, neutrino oscillations may become an important tool for our understanding of the nature of extreme compact objects.

\begin{acknowledgments}
 This work was supported in part by the NSRF via the Program Management Unit for Human Resources and Institutional Development, Research and Innovation [grant number B05F650021]. A.C. is also supported in part by Thailand Science research and Innovation Fund Chulalongkorn University (IND66230009). T.T. was supported by King Mongkut's University of Technology Thonburi's Post-doctoral Fellowship. DM and HC acknowledge support from Nazarbayev University Faculty Development Competitive Research Grant No. 11022021FD2926. HC would also like to thank AC for the kind hospitality during the 10th Bangkok Workshop on High-energy Theory.
\end{acknowledgments}\

\bibliographystyle{apsrev}
\bibstyle{apsrev}