11institutetext: Department of Mechanical, Electrical and Information
Engineering, Shandong University, Weihai, China. 22institutetext: School of
Management, Hangzhou Dianzi University, Hangzhou, China. 33institutetext:
College of Media Engineering, Communication University of Zhejiang, Hangzhou,
China. 44institutetext: School of Computing and Mathematical Sciences,
University of Leicester, Leicester, United Kingdom. 55institutetext:
Department of Radiation Oncology, Sichuan Cancer Hospital and Institution,
Sichuan Cancer Center, School of Medicine, Radiation Oncology Key Laboratory
of Sichuan Province, University of Electronic Science and Technology of China,
Chengdu, China. 66institutetext: School of Cyberspace, Hangzhou Dianzi
University, Hangzhou, China 77institutetext: Suzhou Research Institute of
Shandong University, Suzhou, China.
77email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>††footnotetext: $\dagger$ These authors contributed
equally to this work.

# MMFusion: Multi-modality Diffusion Model for Lymph Node Metastasis Diagnosis
in Esophageal Cancer

Chengyu Wu 11†† 0009-0009-0450-8649 Chengkai Wang 22†† 0009-0000-2690-9608
Yaqi Wang 33 0000-0002-4627-3392 Huiyu Zhou 44 0000-0003-1634-9840 Yatao Zhang
11 0000-0002-6152-0806 Qifeng Wang 55 0000-0001-8544-9007 Shuai Wang 6677
0000-0003-3730-6401

###### Abstract

Esophageal cancer is one of the most common types of cancer worldwide and
ranks sixth in cancer-related mortality. Accurate computer-assisted diagnosis
of cancer progression can help physicians effectively customize personalized
treatment plans. Currently, CT-based cancer diagnosis methods have received
much attention for their comprehensive ability to examine patients’
conditions. However, multi-modal based methods may likely introduce
information redundancy, leading to underperformance. In addition, efficient
and effective interactions between multi-modal representations need to be
further explored, lacking insightful exploration of prognostic correlation in
multi-modality features. In this work, we introduce a multi-modal
heterogeneous graph-based conditional feature-guided diffusion model for lymph
node metastasis diagnosis based on CT images as well as clinical measurements
and radiomics data. To explore the intricate relationships between multi-modal
features, we construct a heterogeneous graph. Following this, a conditional
feature-guided diffusion approach is applied to eliminate information
redundancy. Moreover, we propose a masked relational representation learning
strategy, aiming to uncover the latent prognostic correlations and priorities
of primary tumor and lymph node image representations. Various experimental
results validate the effectiveness of our proposed method. The code is
available at https://github.com/wuchengyu123/MMFusion.

###### Keywords:

Esophageal cancer Feature-guided Diffusion model Multi-modality Lymph node
metastasis

##  1 Introduction

Esophageal cancer, a leading cause of cancer-related mortality globally,
presents significant challenges in diagnosis and treatment planning [23].
Physicians rely on a multifaceted diagnostic approach, including blood tests,
radiomics, comorbidities, tumor markers, and imaging of the primary tumor and
lymph nodes, to formulate a comprehensive assessment of tumor progression [16,
19]. Integrating diverse data sources to identify lymph node metastasis (LNM)
indicators is challenging due to varying diagnostic methodologies among
clinicians, often influenced by personal experience. Therefore, establishing a
reliable system for tumor progression staging is essential to ensure
consistent and objective clinical decision-making and to customize treatment
strategies.

The necessity for utilizing multi-modal data in tumor progression diagnosis
stems from its potential to offer a holistic view of the patient’s condition.
CT imaging serves as a cornerstone for assessing the progression of esophageal
cancer, with prior studies demonstrating that models based on CT images of
tumors or lymph nodes have achieved favorable outcomes in this task [13, 3,
33, 6]. Additionally, to delve into deeper details in medical images, some
researchers extract radiographic features from specific areas and combine
these features with clinical data to enhance the predictive capability for LNM
[11, 29, 28, 21, 2]. However, these studies have generally simply combined
features from different modalities into a new vector [8], a practice that
might neglect the intrinsic medical prior feature relationships between
different type of features [8, 5, 4] and the complex interactions between data
patterns [9].

In the actual diagnostic process for cancer metastasis, clinicians often
assess the condition of the primary tumor and the proximity of nearby lymph
nodes to estimate the current state of metastasis and plan for surgery or
radiation therapy [18]. Our Multi-issue Masked Relational Representation
Learning (MMRL) strategy specifically addresses this clinical approach by
utilizing CT data of tumors with focal lesions and lymph nodes without focal
lesions, recognizing the priority of clinical importance between these
datasets. Thus, the MMRL strategy is deployed to capture the latent feature
relationships between these two types of data, addressing a gap previously
overlooked in research.

Following the MMRL strategy, our Heterogeneous Graph Aggregation (HGA) module
integrates complex cross-modal relationships via self-attention to enhance
multi-modal data fusion. This prepares for our Conditional Feature-guided
Diffusion (CFD) process, which refines the diagnostic process by utilizing
features aggregated by HGA, including CT images of tumors and lymph nodes,
clinical, hematological, and radiomics data to guide the training process. The
CFD process simulates the data generation process, progressively
reconstructing noise distributions to accurately capture the complete multi-
modal conditional distribution of these inputs [22, 7, 15, 30]. Also, the
randomness inherent in the diffusion process minimizes modality information
redundancy [30, 26], improving the model’s predictive accuracy and reliability
[12].

The core contributions of our study are as follows: 1) We collect a multi-
modal dataset of 1,354 cases of ESCC, utilizing various modalities for lymph
node metastasis diagnosis. 2) We employ a CFD process based on a multi-modal
heterogeneous graph to explore data relations, effectively reducing redundancy
in multi-modal features and leading to optimal model performance. 3) We
introduce MMRL strategy, leveraging intra-tissue self-attention and cross-
tissue relational mask self-attention for representational learning to explore
the priority modeling of inter-tissue relationships and the interactive
learning of prognostic information.

##  2 Methods

###  2.1 Architecture Overview

Figure 1: Overview framework of our proposed MMFusion. Initially, the MMRL
strategy is applied to extract correlation from image representations.
Subsequently, we employ HGA to facilitate multi-modal feature interaction.
Finally, the CFD method is utilized for feature redundancy elimination.

In this section, the Multi-modal lymph node Metastasis difFusion (MMFusion)
model is outlined for LNM diagnosis. Initially, as illustrated in Fig. 1, the
gross target volumes (GTV) for the tumor (GTV-T) and lymph node (GTV-N) are
derived by applying corresponding masks to the original CT scans.
Subsequently, latent imaging representations of GTV-N
$\bm{x}_{N}^{i},i\in\\{1,2,3\\}$ and GTV-T $\bm{x}_{T}$ are extracted using a
pretrained ResNet [1] and processed through the MMRL strategy to discern
disease-related information and prognostic correlation priorities between
multi-tissue latent representations. These features, along with data from
hematology, clinical and radiomics sources, are integrated into a
heterogeneous graph. An HGA process, leveraging a Graph Attention Network
(GAT) [25], then identifies potential multimodal feature interactions.
Finally, the CFD approach is applied to eliminate multi-modal feature
redundancy.

###  2.2 Multi-tissue Masked Relational Learning (MMRL) Strategy

Figure 2: Architecture of our proposed Multi-tissue Masked Relational
Representation Learning (MMRL) strategy.

The proposed MMRL strategy explores and learns the prognostic-related
relationship priorities and interactive information among multi-tissues
inspired by [32]. As shown in Fig. 2, commencing the exposition with the GTV-N
branch as a case study, the methodology initiates with the application of
intra-tissue Multi-head Self Attention (MSA) [24] to the concatenated GTV-N
representation, resulting in the generation of an intra-tissue representation
$\bm{x}_{N}^{S}$. Subsequently, a cross-tissue MSA is conducted. In this
phase, the masked relation matrix $\bm{R}_{cross}^{M}$ is derived by applying
random masking to the original relation matrix $\bm{R}_{cross}$. This process
results in the acquisition of both unmasked and masked cross-tissue
representations, designated as $\bm{x}_{N}^{cu}$ and $\bm{x}_{N}^{cm}$,
respectively. For practical application, the masked cross-tissue
representation $\bm{x}_{N}^{cm}$ is merged with $\bm{x}_{N}^{S}$ during the
training phase, whereas $\bm{x}_{N}^{cu}$ is combined with $\bm{x}_{N}^{S}$
during the inference phase. To ensure the integrity of feature preservation,
an alignment loss is implemented between $\bm{x}_{N}^{cu}$ and
$\bm{x}_{N}^{cm}$ to enforce a relation modeling constraint.

###  2.3 Multi-modal Heterogeneous graph-based Conditional Feature-guided
Diffusion model

####  2.3.1 Heterogeneous Graph Aggregation (HGA) Module

In graph theory, a heterogeneous graph is generally denoted as
$\mathcal{G}=(\mathcal{V},\mathcal{E},\mathcal{A},\mathcal{R})$, with vertex
feature matrices $\bm{F}_{a_{i}}\in$
$\mathbb{R}^{\left|\mathcal{V}_{a_{i}}\right|\times d_{a_{i}}}$ for vertex
type $a_{i}\in\mathcal{A}$, and a task to learn the $d$-dimensional vertex
representations $\bm{h}_{v}\in\mathbb{R}^{d}$ for all $v\in\mathcal{V}$ that
can capture structural and relational information involved in $\mathcal{G}$.
$\mathcal{E}$ $\in$ $\mathcal{R}$ denotes set of edges.

Our HGA module focuses on aggregating information on other modalities to
capture potential interactions between multi-modal features. Thus, we define a
heterogeneous graph with only one vertex of each multi-modal data type,
$\forall a_{i}\in\mathcal{A}:\left|\mathcal{V}_{a_{i}}\right|=1$, with GAT
employed to discern the inter-modal interactions. To ensure all possible
relationships between different modalities are considered, we utilize a fully
connected graph. At each layer, the node feature $\bm{F}_{a_{i}}$ will be
updated by its attentional neighbor modality to $\bm{F}_{a_{i}}^{\prime}$ as
follows[35]:

$\displaystyle\bm{F}_{a_{i}}^{\prime}=\sigma\left(\sum_{j\in\\{i\\}\cup\mathcal{N}_{i}}\bm{\alpha}_{ij}\bm{W}\bm{F}_{j}\right).$
(1)

####  2.3.2 Conditional Feature-guided Diffusion (CFD) Approach

For the purpose of eliminating feature redundancy and modeling accurate muti-
modality, a novel diffusion model is developed based on CARD [7]. In the
forward process, the Gaussian noise is added to the ground truth $\bm{y}_{0}$
conditioned by features output from HGA with an arbitrary time step
$t\in[1,T]$, which can be defined as:

$\displaystyle
q\left(\bm{y}_{t}\mid\bm{y}_{0},\bm{f}_{\phi}(\bm{x})\right)=\mathcal{N}\left(\bm{y}_{t};\sqrt{\bar{\alpha}_{t}}\bm{y}_{0}+\left(1-\sqrt{\bar{\alpha}_{t}}\right)\bm{f}_{\phi}(\bm{x}),\left(1-\bar{\alpha}_{t}\right)\mathbb{I}\right).$
(2)

where $\bar{\alpha}_{t}=\prod_{t}\alpha_{t}$, $\alpha_{t}=1-\beta_{t}$ with a
linear noise schedule $\\{\beta_{t}\\}_{t=1\dots T}\in(0,1)$.
$\mathcal{N}(\cdot,\cdot)$ denotes the Gaussian distribution,
$\bm{f}_{\phi}(\bm{x})$ denotes the output feature from HGA, and $\mathbb{I}$
is the identity matrix.

In the reverse process, the denoising UNet $\bm{\epsilon}_{\theta}$ [7] is
applied to learn the noise distribution and generate the final prediction
$\tilde{\bm{y}_{0}}$ during denoising steps from a noisy distribution guided
by the output features from HGA, which can be defined as:

$\displaystyle
p_{\theta}\left(\bm{y}_{0:T-1}\mid\bm{y}_{T},\bm{f}_{\phi}(\bm{x})\right)=\prod_{t=1}^{T}p_{\theta}\left(\bm{y}_{t-1}\mid\bm{y}_{t},\bm{f}_{\phi}(\bm{x})\right).$
(3)

where $\theta$ donates as the parameters in denoising UNet. Concretely, the
reverse process in the inference stage can be formalized as:

$\displaystyle\bm{y}_{t-1}=\frac{1}{\sqrt{\bar{\alpha}_{t}}}\left(\bm{y}_{t}-\left(1-\sqrt{\bar{\alpha}_{t}}\right)\bm{f}_{\phi}(\bm{x})-\sqrt{1-\bar{\alpha}_{t}}\bm{\epsilon}_{\theta}\left(\bm{y}_{t},\bm{f}_{\phi}(\bm{x}),t\right)\right).$
(4)

where $\bar{\alpha}_{t}$, $\bm{f}_{\phi}(\bm{x})$, and
$\bm{\epsilon}_{\theta}$ stand for the same as the previous definition,
$\bm{y}_{t}$ stands for the denoising output from CTD at timestep $t$. By
denoising the process under the guidance of the output features
$\bm{f}_{\phi}(\bm{x})$, the redundancy of multi-modal representations can be
effectively eliminated.

###  2.4 Optimization Method

For model optimization, there are two main parts of loss functions. For the
non-diffusion part, we utilize Binary Cross-Entropy (BCE) as the primary loss
for classification. Also, the Mean Square Error (MSE) loss is utilized for
masked and unmasked representational alignment in the MMRL strategy. The loss
of the non-diffusion part can be defined as:

$\displaystyle\mathcal{L}_{\text{non-
diff}}=\mathcal{L}_{\text{BCE}}(\hat{\bm{y}},\bm{y}_{0})+\mathcal{L}_{\text{MSE}}^{N}(\bm{x}_{N}^{c_{u}},\bm{x}_{N}^{c_{m}})+\mathcal{L}_{\text{MSE}}^{T}(\bm{x}_{T}^{c_{u}},\bm{x}_{T}^{c_{m}}).$
(5)

where $\hat{\bm{y}},\bm{y}_{0}$ are output predictions from the HGA module and
the ground truth. $\bm{x}_{N}^{c_{u}},\bm{x}_{N}^{c_{m}}$ are unmasked and
masked intermediate representations of GTV-N in MMRL strategy. For the
diffusion part, the denoising UNet $\bm{\epsilon}_{\theta}$ is utilized to
predict the noisy distribution, which is added in the forward process. The
loss function of the diffusion part can be defined as:

$\displaystyle\mathcal{L}_{\text{diff}}=\left\|\bm{\epsilon}-\bm{\epsilon}_{\theta}\left(\sqrt{\bar{\alpha}_{t}}\bm{y}_{0}+\sqrt{1-\bar{\alpha}_{t}}\bm{\epsilon}+\left(1-\sqrt{\bar{\alpha}_{t}}\right)\bm{f}_{\phi}(\bm{x}),\bm{f}_{\phi}(\bm{x}),t\right)\right\|^{2}.$
(6)

where $\bm{\epsilon}$ stands for Gaussian noise. $\bar{\alpha}_{t}$,
$\bm{f}_{\phi}(\bm{x})$, $\bm{y}_{0}$, and $\bm{\epsilon}_{\theta}$ stand for
the same as the previous definition. In total, the final loss function can be
defined as:

$\displaystyle\mathcal{L}_{\text{final}}=\mathcal{L}_{\text{non-
diff}}+\mathcal{L}_{\text{diff}}.$ (7)

##  3 Experiments

###  3.1 Dataset Description

A total of 1,354 ESCC patients are enrolled in this study at Sichuan Cancer
Hospital111This study was approved by the Institutional Ethics Committee of
Sichuan Cancer Hospital (No.SCCHEC-02-2020-015).. Corresponding data includes
preoperative contrast-enhanced CT, clinicopathological parameters, and follow-
up information. Contrast-enhanced CT images are retrospectively acquired from
imaging systems manufactured by GE Medical Systems or Philips. To further
assess the model’s robustness, we execute a three-fold cross-validation. The
baseline characteristics of one of the folds can be seen in Table S1.

###  3.2 Main Results

####  3.2.1 Comparison Study

We compare our proposed MMFusion against the state-of-the-art (SOTA) LNM
diagnostic methods to verify the effectiveness of our method. Models are
trained on the training set, selected based on the performance of the
validation set, and final results are reported on the test set. As shown in
Table 1, our proposed model achieves the best performance on all metrics
compared with other methods. Moreover, a paired t-test is conducted, and the
p-value is obtained. Results show that our model significantly outperforms
most of the compared methods, which also demonstrates the excellent
application prospects of our model in the diagnosis of LNM.

Table 1: Comparison study of the proposed method and other SOTA methods.
$\ast$ stands for significant difference (p$<$0.05)

Methods | Accuracy (%) |  Precision (%) |  F1 score (%) |  Recall (%)
---|---|---|---|---
Xie et al. [31] | 86.35 $\pm$ 1.09∗ | 76.30 $\pm$ 11.95 | 51.92 $\pm$ 3.25∗ | 52.83 $\pm$ 1.64∗
Jin et al. [14] | 87.33 $\pm$ 1.36 | 75.91 $\pm$ 3.24∗ | 62.84 $\pm$ 2.74∗ | 55.14 $\pm$ 3.21∗
Wan et al. [27] | 86.22 $\pm$ 0.92∗ | 71.23 $\pm$ 7.98∗ | 63.18 $\pm$ 5.87 | 57.98 $\pm$ 7.82
DensePriNet [34] | 86.22 $\pm$ 0.92∗ | 87.05 $\pm$ 4.12 | 65.21 $\pm$ 3.81 | 54.33 $\pm$ 1.55∗
Ou et al. [20] | 86.84 $\pm$ 0.17∗ | 71.85 $\pm$ 4.60∗ | 61.61 $\pm$ 2.61∗ | 59.20 $\pm$ 2.02
Huang et al. [10] | 87.21 $\pm$ 0.92 | 82.31 $\pm$ 6.31∗ | 56.51 $\pm$ 0.78∗ | 55.89 $\pm$ 0.36∗
MMFusion | 89.14 $\pm$ 0.99 | 93.37 $\pm$ 2.23 | 71.10 $\pm$ 3.73 | 65.97 $\pm$ 3.58

####  3.2.2 Ablation Study

To evaluate the effectiveness of our proposed sub-modules in the MMFusion, we
construct three baseline networks. Base1 removes all components from our
model, where only vanilla ResNet and MLP are included. Incorporating our MMRL
strategy into Base1 yields Base2. Then, by integrating the HGA module into
Base2, we obtain Base3. The results, presented in the table 2, reveal
progressive performance improvements across all metrics with the sequential
addition of each proposed module to Base1. This gradual enhancement
substantiates the individual contribution of each sub-module to the
overarching efficacy of our model.

Table 2: Ablation study of proposed methods. $\ast$ stands for significant
difference (p$<$0.05)

Methods | MMRL | HGA | CFD |  Accuracy (%) |  Precision (%) |  F1 score (%) |  Recall (%)
---|---|---|---|---|---|---|---
Base1 | - | - | - | 86.46 $\pm$ 2.11 | 77.59 $\pm$ 1.60∗ | 64.14 $\pm$ 1.13∗ | 61.16 $\pm$ 1.16
Base2 | $\checkmark$ | - | - | 87.94 $\pm$ 0.70 | 81.49 $\pm$ 9.16 | 66.90 $\pm$ 4.85 | 65.57 $\pm$ 4.22
Base3 | $\checkmark$ | $\checkmark$ | - | 88.65 $\pm$ 0.35 | 84.92 $\pm$ 5.30 | 66.95 $\pm$ 2.66 | 62.46 $\pm$ 2.36
MMFusion | $\checkmark$ | $\checkmark$ | $\checkmark$ | 89.14 $\pm$ 0.99 | 93.37 $\pm$ 2.23 | 71.10 $\pm$ 3.73 | 65.97 $\pm$ 3.58

####  3.2.3 Backbone and Relation Masking Replacement Study

We further verify the performance of our framework under different backbones
with and without relation masking. As Table 3 illustrates, integrating
relation masking significantly enhances model performance, underscoring the
MMRL strategy’s ability to capture latent prognostic-related priorities
information in multiple tissues. Table 3 also demonstrates that different
backbone networks have little impact on the overall performance. This
observation underscores that although backbone architectures lay the
groundwork, the actual performance boost comes from the strategic integration
and optimization of feature fusion techniques.

Table 3: Backbone and relation masking replacement study of proposed methods. CNN | GNN | MMRL |  Accuracy (%) |  Precision (%) |  F1 score (%) |  Recall (%)
---|---|---|---|---|---|---
DenseNet | GCN | w/o mask | 88.12 $\pm$ 0.58 | 88.01 $\pm$ 0.95 | 69.52 $\pm$ 6.91 | 63.87 $\pm$ 1.47
w/ mask | 89.01 $\pm$ 0.15 | 88.24 $\pm$ 3.01 | 70.10 $\pm$ 1.01 | 64.21 $\pm$ 0.72
| GAT | w/o mask | 88.05 $\pm$ 0.52 | 87.53 $\pm$ 3.52 | 66.25 $\pm$ 5.12 | 66.02 $\pm$ 1.05
| w/ mask | 88.65 $\pm$ 1.23 | 91.24 $\pm$ 2.11 | 68.15 $\pm$ 1.05 | 65.55 $\pm$ 2.12
ResNet | GCN | w/o mask | 87.15 $\pm$ 1.01 | 90.11 $\pm$ 0.07 | 68.01 $\pm$ 3.70 | 62.71 $\pm$ 2.15
w/ mask | 87.95 $\pm$ 1.76 | 83.28 $\pm$ 6.84 | 69.12 $\pm$ 3.74 | 66.50 $\pm$ 3.78
| GAT | w/o mask | 88.79 $\pm$ 0.34 | 91.02 $\pm$ 2.77 | 68.11 $\pm$ 5.02 | 63.67 $\pm$ 4.38
| w/ mask | 89.14 $\pm$ 0.99 | 93.37 $\pm$ 2.23 | 71.10 $\pm$ 3.73 | 65.97 $\pm$ 3.58

####  3.2.4 Visualization of Our Diffusion Procedure

To illustrate the reverse process of our CFD, the t-SNE [17] is applied to
visualize the denoised feature embeddings at consecutive time steps. As seen
in Fig. 3, the CFD gradually removes feature redundancy from the multi-modal
feature representation output produced by HGA with time step advanced,
resulting in a more apparent distribution of classes from the Gaussian
distribution. Moreover, we quantify the redundancy elimination effect of the
diffusion model using the DB score. As shown in Fig. 3, the DB score of
clustering diminishes over time steps, indicating the progressive elimination
of multi-modal feature redundancy and leading to enhanced accuracy in LNM
diagnostics.

Figure 3: Denoised multi-modal feature embedding output from HGA using t-SNE.
Red and blue stand for metastasis and non-metastasis, respectively. As the
time step encoding process advances, feature redundancy is gradually
eliminated, resulting in a clear distribution and a lower DB score, which
indicates that our model could effectively perform LNM diagnosis.

##  4 Conclusion

In this work, we introduce MMFusion, a novel framework that integrates primary
tumor and lymph node CT images with clinical, hematology, and radiomics data
for diagnosing lymph node metastasis in esophageal squamous cell cancer.
Specifically, an innovative multi-modal heterogeneous graph-based conditional
feature-guided diffusion model is designed to capture the complex interaction
within multi-modal representations and to minimize feature redundancy across
multi-modalities. Furthermore, we develop a masked relational representation
learning strategy to extract prognostic-related information and establish
correlation priorities among multiple tissues. Sufficient experiments are
conducted and show the effectiveness and superior performance of our proposed
method, highlighting its substantial potential to drive advancements in the
field.

####  4.0.1 Acknowledgments

This work was supported by National Natural Science Foundation of China (No.
62201323 & No. 62206242), Natural Science Foundation of Jiangsu Province (No.
BK20220266), Zhejiang Provincial Natural Science Foundation of China (No.
LDT23F01015F01), National Fund Joint Fund for Regional Innovation and
Development (NO. U20A20386), and Science and Technology Department of Sichuan
Province (No. 2023YFS0488 & No. 2023YFQ0055).

## References

  * [1] Chen, S., Ma, K., Zheng, Y.: Med3d: Transfer learning for 3d medical image analysis. arXiv preprint arXiv:1904.00625 (2019)
  * [2] Cong, M., et al.: Development of a predictive radiomics model for lymph node metastases in pre-surgical ct-based stage ia non-small cell lung cancer. Lung Cancer 139, 73–79 (2020)
  * [3] Elsherif, S.B., et al.: Role of precision imaging in esophageal cancer. Journal of Thoracic Disease 12(9), 5159 (2020)
  * [4] Fass, L.: Imaging and cancer: A review. Molecular Oncology 2 (2008). https://doi.org/10.1016/j.molonc.2008.04.001
  * [5] xia Feng, Q., et al.: An intelligent clinical decision support system for preoperative prediction of lymph node metastasis in gastric cancer. Journal of the American College of Radiology : JACR (2019). https://doi.org/10.1016/j.jacr.2018.12.017
  * [6] Guo, J., et al.: A prospective analysis of the diagnostic accuracy of 3t mri, ct and endoscopic ultrasound for preoperative t staging of potentially resectable esophageal cancer. Cancer Imaging 20 (2020). https://doi.org/10.1186/s40644-020-00343-w
  * [7] Han, X., Zheng, H., Zhou, M.: Card: Classification and regression diffusion models. Advances in Neural Information Processing Systems 35, 18100–18115 (2022)
  * [8] Hu, D., Li, S., Wu, N., Lu, X.: A multi-modal heterogeneous graph forest to predict lymph node metastasis of non-small cell lung cancer. IEEE Journal of Biomedical and Health Informatics 27(3), 1216–1224 (2023)
  * [9] Hu, Z., et al.: Enhancing the accuracy of lymph-node-metastasis prediction in gynecologic malignancies using multimodal federated learning: Integrating ct, mri, and pet/ct. Cancers 15 (2023). https://doi.org/10.3390/cancers15215281
  * [10] Huang, C., et al.: Development and validation of a deep learning model to predict survival of patients with esophageal cancer. Frontiers in Oncology 12, 971190 (2022)
  * [11] Huang, Y.q., et al.: Development and validation of a radiomics nomogram for preoperative prediction of lymph node metastasis in colorectal cancer. Journal of clinical oncology 34(18), 2157–2164 (2016)
  * [12] Hunter, J., Craig, P., Phillips, H.: On the use of random walk models with spatially variable diffusivity. Journal of Computational Physics 106, 366–376 (1993). https://doi.org/10.1016/S0021-9991(83)71114-9
  * [13] Jayaprakasam, V.S., et al.: Role of imaging in esophageal cancer management in 2020: update for radiologists. American Journal of Roentgenology 215(5), 1072–1084 (2020)
  * [14] Jin, C., et al.: Deep learning analysis of the primary tumour and the prediction of lymph node metastases in gastric cancer. British Journal of Surgery 108(5), 542–549 (2021)
  * [15] Kawar, B., Elad, M., Ermon, S., Song, J.: Denoising diffusion restoration models. Advances in Neural Information Processing Systems 35, 23593–23606 (2022)
  * [16] Kim, T., Kim, H.Y., Lee, K., Kim, M.S.: Multimodality assessment of esophageal cancer: preoperative staging and monitoring of response to therapy. Radiographics : a review publication of the Radiological Society of North America, Inc 29 2, 403–21 (2009). https://doi.org/10.1148/rg.292085106
  * [17] Van der Maaten, L., Hinton, G.: Visualizing data using t-sne. Journal of machine learning research 9(11) (2008)
  * [18] Marino, M., Avendaño, D., Zapata, P., Riedl, C., Pinker, K.: Lymph node imaging in patients with primary breast cancer: Concurrent diagnostic tools. The Oncologist 25, e231 – e242 (2019). https://doi.org/10.1634/theoncologist.2019-0427
  * [19] Ojiri, H.: Diagnostic imaging of the esophageal cancer. Esophageal Squamous Cell Carcinoma (2020). https://doi.org/10.1007/978-4-431-54977-2_3
  * [20] Ou, J., et al.: Ct radiomics features to predict lymph node metastasis in advanced esophageal squamous cell carcinoma and to discriminate between regional and non-regional lymph node metastasis: a case control study. Quantitative Imaging in Medicine and Surgery 11(2), 628 (2021)
  * [21] Sha, X., Gong, G., Qiu, Q., Duan, J., Li, D., Yin, Y.: Discrimination of mediastinal metastatic lymph nodes in nsclc based on radiomic features in different phases of ct imaging. BMC medical imaging 20(1), 1–8 (2020)
  * [22] Song, J., Meng, C., Ermon, S.: Denoising diffusion implicit models. arXiv preprint arXiv:2010.02502 (2020)
  * [23] Sung, H., et al.: Global cancer statistics 2020: Globocan estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA: a cancer journal for clinicians 71(3), 209–249 (2021)
  * [24] Vaswani, A., et al.: Attention is all you need. Advances in neural information processing systems 30 (2017)
  * [25] Veličković, P., et al.: Graph attention networks. arXiv preprint arXiv:1710.10903 (2017)
  * [26] Veraart, J., Novikov, D., Christiaens, D., Ades-aron, B., Sijbers, J., Fieremans, E.: Denoising of diffusion mri using random matrix theory. NeuroImage 142, 394–406 (2016). https://doi.org/10.1016/j.neuroimage.2016.08.016
  * [27] Wan, L., et al.: Prediction of lymph node metastasis in stage t1–2 rectal cancers with mri-based deep learning. European Radiology 33(5), 3638–3646 (2023)
  * [28] Wang, X., et al.: Can peritumoral radiomics increase the efficiency of the prediction for lymph node metastasis in clinical stage t1 lung adenocarcinoma on ct? European radiology 29, 6049–6058 (2019)
  * [29] Wu, S., et al.: A radiomics nomogram for the preoperative prediction of lymph node metastasis in bladder cancer. Clinical Cancer Research 23(22), 6904–6911 (2017)
  * [30] Xiao, Z., Kreis, K., Vahdat, A.: Tackling the generative learning trilemma with denoising diffusion gans. arXiv preprint arXiv:2112.07804 (2021)
  * [31] Xie, C., et al.: Prediction of individual lymph node metastatic status in esophageal squamous cell carcinoma using routine computed tomography imaging: comparison of size-based measurements and radiomics-based models. Annals of Surgical Oncology 29(13), 8117–8126 (2022)
  * [32] Yang, Q., Li, W., Li, B., Yuan, Y.: Mrm: Masked relation modeling for medical image pre-training with genetics. In: Proceedings of the IEEE/CVF International Conference on Computer Vision. pp. 21452–21462 (2023)
  * [33] Yao, J., Ye, X., Xia, Y., Zhou, J., Shi, Y., Yan, K., Wang, F., Lin, L., Yu, H., Hua, X.S., Lu, L., Jin, D., Zhang, L.: Effective opportunistic esophageal cancer screening using noncontrast ct imaging. In: Wang, L., Dou, Q., Fletcher, P.T., Speidel, S., Li, S. (eds.) Medical Image Computing and Computer Assisted Intervention – MICCAI 2022. pp. 344–354. Springer Nature Switzerland, Cham (2022)
  * [34] Zhao, X., et al.: A cross-modal 3d deep learning for accurate lymph node metastasis prediction in clinical stage t1 lung adenocarcinoma. Lung Cancer 145, 10–17 (2020)
  * [35] Zhu, S., Pan, S., Zhou, C., Wu, J., Cao, Y., Wang, B.: Graph geometry interaction learning. Advances in Neural Information Processing Systems 33, 7548–7558 (2020)


# Gravitational orbital Hall effect of vortex photons in Lense-Thirring metric

Wei-Si Qiu  Dan-Dan Lian<EMAIL_ADDRESS>Peng-Ming Zhang  School of
Physics and Astronomy, Sun Yat-sen University, 519082 Zhuhai, China

###### Abstract

Vortex photons, possessing an intrinsic orbital angular momentum (OAM) aligned
with the direction of propagation, are described using vortex electromagnetic
wave packets. Similar to the gravitational spin Hall effect (SHE), these
vortex photons are expected to exhibit intrinsic OAM-dependent trajectories
and separations when propagating through a gravitational field, a phenomenon
termed the gravitational orbital Hall effect (OHE). In this work, we construct
a vortex Laguerre-Gaussian electromagnetic wave packet and analyze its motion
by solving covariant Maxwell equations within the Lense-Thirring metric. Our
findings reveal that vortex photons with different intrinsic OAM not only
separate perpendicular to the null geodesic plane but also within it. This
behavior contrasts with the gravitational SHE, where photons of opposite spins
separate primarily perpendicular to the null geodesic plane. Moreover, the
relationship between the separation and intrinsic OAM differs significantly
from that between the separation and spin. These results suggest a unique
interaction between intrinsic OAM and gravity, distinct from the spin-gravity
coupling, indicating that the gravitational OHE might not be precisely
predicted by merely substituting spin with intrinsic OAM in the gravitational
SHE.

##  I Introduction

The spin Hall effect (SHE), characterized by spin-dependent trajectory
deviations of particles due to spin-orbit coupling in specific media or
external fields, is a well-documented phenomenon in light propagation through
gradient-index media and various optical materials. This effect has been
thoroughly discussed, revealing both its detectability and foundational
mechanisms Bliokh (2006); Bliokh  _et al._ (2008); Bliokh and Desyatnikov
(2009); Zhou  _et al._ (2013); Liu  _et al._ (2015); Bliokh  _et al._ (2015);
Ling  _et al._ (2023). Recent advancements have investigated the ability to
control the SHE through crystal birefringence, highlighting its utility in the
precise manipulation of micro-particles Fu  _et al._ (2019) and the behaviors
in optical fibers with curved step-index profiles Mieling and Oancea (2023).
Authoritative reviews have effectively summarized the core concepts and varied
applications of the optical SHE, significantly enhancing our comprehension of
this field Ling  _et al._ (2017); Liu  _et al._ (2022); Sheng  _et al._
(2023).

In gravitational fields, particles also exhibit trajectories and deviations
that depend on their spin, due to gravitational spin-orbit coupling, an effect
known as the gravitational SHE. This phenomenon was initially identified by
Papapetrou and Corinaldesi using the Mathisson-Papapetrou-Dixon (MPD)
equations Papapetrou (1951); Corinaldesi and Papapetrou (1951); Dixon (1970a,
b, 1974, 2015). Subsequent studies have extensively investigated the
gravitational SHE across various metrics employing these equations Carmeli
_et al._ (1977); Plyatsko  _et al._ (2011); Hackmann  _et al._ (2014);
Antoniou  _et al._ (2019); Duval  _et al._ (2019); Zhang  _et al._ (2020).
Nonetheless, the derivation of unique trajectories from the MPD equations
necessitates the application of artificial constraints to the angular momentum
$S^{\mu\nu}$. These constraints lead to different trajectories, highlighting
the influence of the chosen constraints Costa and Natário (2015).

To avoid this issue, researchers have employed alternative methods to study
the gravitational SHE. These methods include applying the Wentzel-Kramers-
Brillouin (WKB) approximation to Maxwell’s equations Frolov and Shoom (2012);
Oancea  _et al._ (2020); Deriglazov (2021); Harte and Oancea (2022); Andersson
and Oancea (2023); Oancea and Harko (2024), utilizing the Energy-Momentum
Tensor (EMT) Lian and Chen (2022); Lian and Zhang (2023), and using Foldy-
Wouthuysen transformation to derive non-relativistic Hamiltonians Obukhov
(2001); Silenko and Teryaev (2005); Obukhov  _et al._ (2009, 2013); Wang
(2024). These studies not only deepen our understanding of the interactions
between matter and gravity but also provide numerous predictions that can be
used to test gravitational theories.

Beyond possessing spin angular momentum (SAM), a single particle can also
exist in a vortex state, carrying an intrinsic orbital angular momentum (OAM)
aligned with its average momentum. The vortex state, described by the wave
function $\psi\propto e^{i\ell\phi}$, is an eigenstate of $z$-component of the
OAM operator, $\hat{L}_{z}=-i\hbar\frac{\partial}{\partial\phi}$, with
$\ell\hbar$ as its eigenvalue Landau and Lifshitz (2013); Cohen and Tannoudji
(1977). Here, $\phi$ is the azimuthal angle. These particles, when in such
vortex states, are referred to as ”vortex” particles, highlighting their
intrinsic OAM. Vortex photons were first discoveried by Allen et al. in 1992
Beijersbergen  _et al._ (1993). Since then, vortex photons have been generated
in a variety of optical experiments and used in industrial laser production
Bliokh  _et al._ (2015, 2017), demonstrating their applicability in numerous
areas of fundamental physics. These applications include optical trapping of
particles Gahagan and Swartzlander (1996); Garcés-Chávez  _et al._ (2003),
quantum information and communications Mair  _et al._ (2001); Leach  _et al._
(2010), astronomy and astrophysics Foo  _et al._ (2005); Harwit (2003), and
optical solitons Kruglov and Vlasov (1985); Swartzlander Jr and Law (1992).

Similar to the SHE, the dynamics of particles are also affected by intrinsic
OAM when passing through materials, a phenomenon termed the orbital Hall
effect (OHE). The OHE has been explored in photons moving along spatial curves
Lai  _et al._ (2018) and observed in electrons in titanium metals Choi  _et
al._ (2023). These findings lead to the expectation that vortex photons,
possessing intrinsic OAM, will similarly exhibit altered dynamics in
gravitational fields, an effect we term the gravitational OHE. However, the
MPD equations and WKB approximation, typically used for the gravitational SHE
analysis, are less effective for the gravitational OHE due to their inadequate
description of vortex states and inability to distinguish between SAM and
intrinsic OAM.

In this study, we focus on exploring the gravitational OHE for vortex photons
and its potential distinctions from the gravitational SHE. Our investigation
employs an alternative methodology comprising three main steps. First, we
construct an electromagnetic wave packet to model the initial state of a
vortex photon. Next, we examine the photon’s evolution in a gravitational
field by solving the covariant Maxwell equations. Finally, we assess the
motion of the vortex photon by analyzing the center of its energy-momentum
tensor (EMT). This approach allows us to compare the trajectories of vortex
photons with varying intrinsic OAM, thereby shedding light on the
gravitational OHE specific to vortex photons. This methodology has been
effectively utilized in previous studies to analyze the gravitational SHE of
photons within the Schwarzschild metric Lian and Chen (2022) and the
gravitational OHE of scalar vortex particles in a stellar gravitaional field
Lian and Zhang (2023).

The structure of the paper is organized as follows: Section II introduces the
vortex photon as a Laguerre-Gaussian electromagnetic wave packet and explores
its evolution by numerically solving the Maxwell equations within the Lense-
Thirring metric. In Section III, we examine the gravitational OHE for vortex
photons in this metric, using the center of the EMT. This section also
compares the gravitational OHE with the gravitational SHE and seek out the
differences between these two phenomena. Section IV offers further discussion
on our findings. Throughout this paper, we adopt the metric signature
$(-,+,+,+)$ and normalize units such that $c=\hbar=1$.

##  II Dynamics of vortex photons in the Lense-Thirring metric

The dynamics of vortex photons, represented as electromagnetic wave packets,
are governed by the Maxwell equations. Within a gravitational field, these
equations are expressed as

$\nabla_{\mu}F^{\mu\nu}=0.$ (1)

When considering the motion of these photons outside a gravitational source,
where the Ricci tensor $R_{\alpha\beta}$ vanishes, employing the covariant
Lorenz gauge condition, $\nabla_{\mu}A^{\mu}=0$, simplifies Eq. (1) to

$\nabla^{\mu}\nabla_{\mu}A^{\nu}=0.$ (2)

The rotation of celestial bodies significantly influences their gravitational
fields. In the weak field limit, the gravitational field around a slowly
rotating star is well-described by the Lense-Thirring metric Adler (2015);
Wald (2010). In Cartesian coordinates, the proper time $\text{d}s$ is defined
as:

$\mathrm{d}s^{2}=-\left(1-\frac{2GM}{r}\right)\mathrm{d}t^{2}+\left(1+\frac{2GM}{r}\right)(\mathrm{d}x^{2}+\mathrm{d}y^{2}+\mathrm{d}z^{2})-\frac{4GMa}{r^{3}}(x\mathrm{d}y\mathrm{d}t-y\mathrm{d}x\mathrm{d}t),$
(3)

where $M$ represents the mass of the gravitational source,
$r=\sqrt{x^{2}+y^{2}+z^{2}}$ the distance from the source’s center, and
$a=|\vec{J}|/M$ the angular momentum per unit mass, with the source’s angular
momentum $\vec{J}$ aligned along the $z$-axis. The Lense-Thirring metric is
thus approximated by:

$g_{\mu\nu}\simeq\eta_{\mu\nu}+h_{\mu\nu},$ (4)

with $\eta_{\mu\nu}=\text{diag}(-1,1,1,1)$ and the perturbation $h_{\mu\nu}$
as

$h_{\mu\nu}=\begin{pmatrix}\frac{2GM}{r}&\frac{2GMay}{r^{3}}&-\frac{2GMax}{r^{3}}&0\\\
\frac{2GMay}{r^{3}}&\frac{2GM}{r}&0&0\\\
-\frac{2GMax}{r^{3}}&0&\frac{2GM}{r}&0\\\ 0&0&0&\frac{2GM}{r}\\\
\end{pmatrix}$ (5)

Within the framework of the weak field approximation, where $GM/r\ll 1$, it is
practical to consider the gravitational impact on the electromagnetic field
$A^{\nu}$ as a perturbative effect:

$A^{\nu}=\bar{A}^{\nu}+\widetilde{A}^{\nu},$ (6)

with $\bar{A}^{\nu}$ representing the zeroth-order term in $GM/r$, satisfying
the equation $\eta^{\rho\mu}\partial_{\rho}\partial_{\mu}\bar{A}^{\nu}=0$.
According to Eq. (2), the perturbation $\widetilde{A}^{\nu}$, being first
order in $GM/r$, is approximately described by

$\partial_{\lambda}\partial^{\lambda}\widetilde{A}^{\nu}\simeq
h^{\lambda\alpha}\partial_{\lambda}\partial_{\alpha}\bar{A}^{\nu}+\eta^{\lambda\alpha}\widetilde{\Gamma}^{\beta}_{\lambda\alpha}\partial_{\beta}\bar{A}^{\nu}-\partial^{\alpha}\widetilde{\Gamma}^{\nu}_{\alpha\beta}\bar{A}^{\beta}-2\widetilde{\Gamma}^{\nu}_{\alpha\beta}\partial^{\alpha}\bar{A}^{\beta},$
(7)

where
$\widetilde{\Gamma}_{\lambda\alpha}^{\beta}=\frac{1}{2}\eta^{\beta\gamma}(\partial_{\alpha}h_{\lambda\gamma}+\partial_{\lambda}h_{\alpha\gamma}-\partial_{\gamma}h_{\lambda\alpha})$
denotes the first-order terms in $GM/r$ of the affine connection.

###  II.1 Physical system and vortex photons’ initial states

The physical system under consideration is depicted in Fig. 1, featuring a
gravitational source located at $(-b,0,0)$ that rotates around the $z$-axis.
Initially positioned at the coordinate origin, a vortex photon begins its free
fall with momentum $\vec{p}$ directed along the $z$-axis. Both its SAM
$\vec{S}$ and intrinsic OAM $\vec{L}$ are aligned with the $z$-axis.

Figure 1: Configurations of the physical system. The yellow sphere depicts the
rotating gravitational source situated at $(-b,0,0)$. A vortex photon,
symbolized by the orange sphere, initially positioned at the coordinate
origin, starts its free fall along the $z$-axis. The dashed line $x(\tau)$
symbolizes the vortex photon’s trajectory.

In the absence of gravity within Minkowski spacetime, a vortex photon can be
characterized by a vortex Laguerre-Gaussian electromagnetic wave packet Liu
and Ivanov (2023). The expression for this wave packet in momentum space is
given by

$\bar{A}_{f}^{i}(\vec{k},t)=\frac{\sqrt{2}N(\sigma_{\perp}k_{\perp})^{\ell}}{\omega^{3/2}}\epsilon^{i}\exp\left(-\frac{k_{\perp}^{2}\sigma_{\perp}^{2}}{2}-\frac{(k_{z}-p_{z})^{2}\sigma_{z}^{2}}{2}+i\ell\phi_{k}-i\omega
t\right),$ (8)

where $\vec{\epsilon}=\\{\epsilon^{i}\\}=(k_{z}^{2}+k_{y}^{2}-i\sigma
k_{x}k_{y},i\sigma k_{z}^{2}+i\sigma k_{x}^{2}-k_{x}k_{y},-i\sigma
k_{y}k_{z}-k_{x}k_{z})^{\text{T}}$, with $\sigma=0,\pm 1$, represents the
polarization vector. The wave packet adheres to the Coulomb Gauge,
$\partial_{i}\bar{A}^{i}=0$, with $\vec{A}^{0}=0$. The wave vector’s
transverse component is denoted as $k_{\perp}=\sqrt{k_{x}^{2}+k_{y}^{2}}$. The
angular frequency, $\omega$, is defined as
$\omega=\sqrt{k_{x}^{2}+k_{y}^{2}+k_{z}^{2}}$, and the azimuthal angle,
$\phi_{k}$, is given by $\phi_{k}=\arctan(\frac{k_{y}}{k_{x}})$. The parameter
$p_{z}$ refers to the average initial momentum of the wave packet along the
$z$-axis, and $N$ is the normalization factor.

It is verified that the wave packet, as defined by Eq. (8), carries a SAM of
$\sigma$ and an intrinsic OAM of $\ell$. Additionally, it is established that
the initial momentum of the wave packet is directed along the $z$-axis and its
initial position is at the coordinate origin. Thus, this wave packet aligns
with the configurations of the physical system considered in this study.

For a realistic vortex photon, the radius of its wave packet is expected to be
significantly smaller than its distance, $b$, from the center of the
gravitational source. For the wave packet described in Eq. (8), its radius,
$R$, can be approximated by $R\sim\sigma_{\perp}\sim\sigma_{z}$. Hence, we
establish the following relationship:

$\sigma_{\perp}\sim\sigma_{z}\ll b.$ (9)

This indicates that the wave packet is confined to a significantly small
volume in three-dimensional space, with its radius being much smaller than the
distance to the gravitational source. Consequently, the gravitational field
within the wave packet can be approximated as follows:

$\displaystyle h_{00}\simeq\frac{2GM}{b}\left(1-\frac{x}{b}\right),\quad
h_{ij}\simeq\frac{2GM}{b}\left(1-\frac{x}{b}\right)\delta_{ij},$
$\displaystyle h_{01}=h_{10}\simeq\frac{2GMay}{b^{3}},\quad
h_{02}=h_{20}\simeq-\frac{2GMa}{b^{2}}\left(1-\frac{2x}{b}\right),$ (10)

where only terms of the first order in $x^{i}/b$ have been considered.

While the gravitaional perturbation $h_{\mu\nu}$ has been significantly
simplified compared to its form in Eq. (5), Eq. (7) remains a second-order
inhomogeneous partial differential equation. This complexity presents
substantial challenges for direct solution. Additionally, the vortex photon’s
wave packet is articulated in momentum space, necessitating the transformation
of Eq. (7) into momentum space through a Fourier Transform across the three-
dimensional space. Given the wave packet $\bar{A}_{f}^{i}$ as depicted in Eq.
(8), the momentum space representation of Eq. (7) can be expressed as:

$\left(\partial_{t}^{2}+\omega^{2}\right)\tilde{A}_{f}^{\rho}=a_{1}^{\rho}(\vec{k})e^{-i\omega
t}+a_{2}^{\rho}(\vec{k})e^{-i\omega t}t,$ (11)

where
$\widetilde{A}^{\rho}_{f}=\frac{1}{2\pi^{3}}\int\widetilde{A}^{\rho}\exp(-i\vec{k}\cdot\vec{x})\text{d}^{3}x$
represents the perturbation term in momentum space. The coefficient
$a_{1}^{\rho}(\vec{k})$ encompasses all factors associated with terms that
include the time-dependent factor $\exp(-i\omega t)$. Similarly, the
coefficient $a_{2}^{\rho}(\vec{k})$ encapsulates all factors related to terms
that incorporate the time-dependent expression $t\exp(-i\omega t)$. Detailed
expressions of these equations are available in Appendix A.

Given the initial conditions:

$\left.\tilde{A}_{f}^{\rho}\right|_{t=0}=0\quad\text{and}\quad\left.\partial_{t}\tilde{A}_{f}^{\rho}\right|_{t=0}=0,$
(12)

Eq. (11) can be analytically solved. The perturbation term
$\widetilde{A}^{\rho}_{f}$ is expressed as:

$\tilde{A}_{f}^{\rho}=\frac{a_{1}^{\rho}(\vec{k})(1-e^{2i{\omega}t}+2i{\omega}t)e^{-i{\omega}t}}{4{\omega}^{2}}+\frac{a_{2}^{\rho}(\vec{k})(2{\omega}t-i+2i{\omega}^{2}t^{2}+ie^{2i{\omega}t})e^{-i{\omega}t}}{8{\omega}^{3}}.$
(13)

Accordingly, the wave packet of the free-falling vortex photons in coordinate
space is determined through the inverse Fourier Transform:

$A^{\rho}=\int{A}^{\rho}_{f}\exp(i\vec{k}\cdot\vec{x})\,\text{d}^{3}k,\quad
A^{\rho}_{f}=\bar{A}^{\rho}_{f}+\widetilde{A}^{\rho}_{f}.$ (14)

Due to the complexity of the wave packet $\bar{A}^{i}_{f}$, the inverse
Fourier Transform cannot be analytically performed. Therefore, we employ the
numerical Fast-Fourier-Transform (FFT) technique to calculate this
transformation and numerically obtain the wave packet $A^{\rho}$ in coordinate
space. Following these numerical computations, the energy-momentum tensor
(EMT) and the center of the vortex photon can be straightforwardly calculated.

##  III Gravitational Orbital Hall Effect for vortex photons

The dynamics of vortex photons are described through electromagnetic wave
packets, which preclude their simplification to point-like particles.
Consequently, traditional trajectories applicable to point-like particles
might not be suitable for vortex photons. In this study, we characterize the
motion of vortex photons by employing EMT, thereby defining their trajectories
via the EMT:

$\langle
x^{i}\rangle^{\mu\nu}=\frac{\int\sqrt{-g}x^{i}T^{\mu\nu}\text{d}x^{3}}{\int\sqrt{-g}T^{\mu\nu}\text{d}x^{3}},$
(15)

where the superscripts $\mu$ and $\nu$ denote the components of the EMT used
to define the vortex photons’ center. Here, $g=\det(g_{\mu\nu})$ represents
the determinant of the metric tensor $g_{\mu\nu}$.

Under the weak field approximation, where $h_{\mu\nu}\ll 1$, the gravitational
influence on the evolution of vortex photons is captured by the perturbation
term $\widetilde{A}^{\rho}$, which is first-order in $h_{\mu\nu}$.
Consequently, the EMT for vortex photons can be approximated as
$T^{\mu\nu}\simeq\bar{T}^{\mu\nu}+\widetilde{T}^{\mu\nu}$. Here,
$\bar{T}^{\mu\nu}$, the zeroth-order term in $h_{\mu\nu}$, represents the free
motion of vortex photons in the absence of gravitational effects. In contrast,
$\widetilde{T}^{\mu\nu}$, the first-order term in $h_{\mu\nu}$, accounts for
the gravitational influence on their motion. This allows for a simplified
expression for the center of vortex photons as:

$\langle x^{i}\rangle^{\mu\nu}\simeq\frac{\int
x^{i}\widetilde{T}^{\mu\nu}\text{d}x^{3}}{\int\bar{T}^{\mu\nu}\text{d}x^{3}}.$
(16)

In this study, we focus on utilizing the energy density of the symmetric EMT
to explore the gravitational OHE. Previous studies have indicated that the
choice of components and formulations of the EMT significantly influences the
position of the electromagnetic wave packet’s center in gravitational fields,
affecting the gravitational SHE as well Lian and Chen (2022). These findings
suggest that the gravitational OHE might be similarly affected by the choice
of the EMT components and formulations. Accordingly, we also compute the
gravitational OHE using different EMT components and expressions, facilitating
a comparative analysis with the gravitational SHE.

###  III.1 Separations of vortex photons with different intrinsic OAM

Figure 2: Trajectories and separations of vortex photons with intrinsic OAM
$\ell$. The upper panels illustrate the separations of vortex photons within
the null geodesic plane, defined as $\Delta x_{\ell}=\langle
x\rangle_{\ell}-\langle x\rangle_{\ell=0}$ and $\Delta z_{\ell}=\langle
z\rangle_{\ell}-\langle z\rangle_{\ell=0}$. The lower left panel presents the
trajectories of vortex photons in the transverse direction, perpendicular to
the null geodesic plane. The lower right panel shows the transverse trajectory
deviations between rotating ($a\neq 0$) and non-rotating ($a=0$) scenarios,
denoted by $\Delta y_{a}=\langle y\rangle_{a}-\langle y\rangle_{a=0}$. In
these illustrations, parameters are nondimensionalized using $b$ and specified
as: $GM/b=10^{-4}$, $p_{z}b=10^{7}$, and
$\sigma_{x}b=\sigma_{y}b=\sigma_{z}b=5\times 10^{3}$.

The results depicted in Fig. 2 demonstrate that the rotation of the
gravitational source affects the transverse trajectories of photons indeed.
However, for vortex photons moving within the same rotating gravitational
field, the alterations to their transverse trajectories by the source’s
rotation are consistent. As a result, the separations between vortex photons
with different intrinsic OAM are essentially unchanged. This indicates that
the gravitational OHE exhibits limited sensitivity to the gravitational
source’s rotation. Given this insensitivity, the discussion will not
extensively cover the impact of rotation but will instead concentrate on
investigating the differences between the gravitational spin and orbital Hall
effects.

Unlike the gravitational SHE, where photons of different spins primarily
diverge along the transverse direction perpendicular to the null geodesic
plane, vortex photons demonstrate separations both transversely and within the
null geodesic plane itself. This distinction highlights a fundamental
difference between the gravitational OHE and the gravitational SHE, with the
former affecting motion in a broader context.

To compare these two Hall effects in detail, we examine the transverse motion
of vortex photons. For photons whose travel time $t$ is much less than their
initial distance from the gravitational source ($t\ll b$), their transverse
trajectory can be approximated as:

$\langle y\rangle_{\ell}\simeq-\frac{GM\lambda\ell}{2\pi b^{2}}t,$ (17)

where $\lambda=2\pi/p_{z}$ represents the photon’s wavelength. In contrast,
the gravitational SHE influences the transverse trajectory as described by
Lian et al. Lian and Chen (2022):

$\langle y\rangle_{\sigma}\simeq\frac{2GM\lambda\sigma}{2\pi b^{2}}t.$ (18)

Eqs. (17) and (18) underline the difference between the gravitational spin and
orbital Hall effects. Specifically, while the gravitational SHE induces a
transverse movement given by $2GM\sigma\lambda t/b^{2}$, the gravitational OHE
dictates a transverse movement following $-GM\ell\lambda t/b^{2}$. Therefore,
a photon characterized by both a spin $\sigma$ and an intrinsic OAM $\ell$
will have its transverse trajectory determined by:

$\langle y\rangle_{\sigma,\ell}\simeq\frac{(2\sigma-\ell)GM\lambda t}{2\pi
b^{2}}.$ (19)

As previously discussed, the gravitational spin and orbital Hall effects are
significantly influenced by the selection of components and expressions of the
Energy-Momentum Tensor (EMT). Consequently, the transverse trajectory of a
vortex photon, carrying a total angular momentum $\sigma+\ell$, can be
generalized as follows:

$\langle y\rangle_{\sigma,\ell}\simeq\frac{(\alpha\sigma+\beta\ell)GM\lambda
t}{2\pi b^{2}}.$ (20)

Here, $\alpha$ and $\beta$ are coefficients that depend on the specific choice
of EMT component and expression used to define the vortex photon’s center. The
precise values of $\alpha$ and $\beta$ corresponding to these different
choices are detailed in Table 1.

Definitions of the center | $\alpha$ | $\beta$
---|---|---
Energy density of symmetric EMT | 2.0 (2.0) | -1.0 (-1.0)
Energy flux density of symmetric EMT | 1.0 (1.0) | -2.0 (-2.0)
Energy density of canonical EMT | 1.5 (1.5) | -1.0 (-1.0)
Energy flux density of canonical EMT | 0.5 (0.5) | -2.0 (-2.0)
Table 1: Values of $\alpha$ and $\beta$ for photons in vortex Laguerre-
Gaussian and Bessel states. Values in bold parentheses correspond to photons
in vortex Bessel states.

The results presented in Table 1 demonstrate that both the gravitational SHE
and OHE are influenced by the choice of EMT components, suggesting a
dependence of these phenomena on the detection method. While the gravitational
SHE is sensitive to the specific expressions of the EMT, the gravitational OHE
does not exhibit such sensitivity. Furthermore, neither the gravitational SHE
nor the OHE exhibit sensitivity to the details of the vortex wave packets.
These observations imply that when coupled to gravity, the intrinsic OAM of
photons may not be directly equivalent to their spin, highlighting the
complexity of matter-gravity coupling.

Notably, recent findings reveal that a photon with spin $\sigma$, moving from
a finite distance $r_{0}$ to infinity, diverges from the null geodesic plane
by an angle $\theta_{\sigma}$. This angle is approximately described by
$\theta_{\sigma}\simeq\alpha GM\sigma\lambda/2\pi r^{2}_{0}$ Duval  _et al._
(2019); Lian and Chen (2022), a deviation attributed to the gravitational SHE.
Given the gravitational OHE’s similarity to the gravitational SHE near the
gravitational source, analogous behavior is expected for a vortex photon with
a total angular momentum $\sigma+\ell$. As such a vortex photon moves from
$r_{0}$ to infinity, its deviation angle $\theta_{\sigma+\ell}$ from the null
geodesic plane can be approximated as:

$\theta_{\sigma+\ell}\simeq\frac{(\alpha\sigma+\beta\ell)GM\lambda}{2\pi
r^{2}_{0}}.$ (21)

This equation indicates that photons carrying identical total angular momenta,
$\sigma+\ell$, might also separate to each other. Consider three photons, each
with a total angular momentum of $\sigma+\ell=1$. The first photon has a spin
$\sigma=1$ and an intrinsic OAM $\ell=0$; the second photon has $\sigma=0$ and
$\ell=1$; and the third has $\sigma=-1$ and $\ell=2$. The angles
$\theta_{\sigma+\ell}$ corresponding to these three photons are given by:

$\theta_{1,0}\simeq\frac{2GM\lambda}{2\pi
r_{0}^{2}},\quad\theta_{0,1}\simeq-\frac{GM\lambda}{2\pi
r_{0}^{2}},\quad\theta_{-1,2}\simeq-\frac{4GM\lambda}{2\pi r_{0}^{2}},$ (22)

where the center of the vortex photon is defined by using the energy density
of the symmetric EMT. Hence, a photon’s trajectory is influenced not only by
its total angular momentum but also by the specific contributions of its spin
and intrinsic OAM.

###  III.2 SAM VS intrinsic OAM

Photon dynamics are governed by Maxwell’s equations and can be described
through electromagnetic wave packets. For a spin-polarized photon, its SAM is
characterized by the polarization of the electromagnetic field $A^{\mu}$,
specifically represented by a polarization vector $\vec{\epsilon}$.
Importantly, $\vec{\epsilon}$ could be independent of its wave packet’s
spatial distribution, as suggested by left-handed and right-handed circularly
polarized photons. This independence suggests that its SAM remains physically
meaningful when a photon is approximated as a point-like particle.

However, for a vortex photon, its intrinsic OAM is captured by a phase factor
$\psi(\vec{x},t)\propto\exp(i\ell\phi)$ within the electromagnetic wave
packet, where $\phi$ is the azimuthal angle. This angle, and consequently the
phase factor, are dependent on the wave packet’s spatial distribution. Thus,
when the wave packet’s size is theoretically reduced to zero, the phase factor
becomes undefined, suggesting that the intrinsic OAM might lose its physical
meaning in the approximation of a photon as a point-like particle.

The MPD equations are widely employed to analyze the dynamics of particles
with angular momentum, particularly for describing the motion of spin-
polarized photons. In these equations, photon characteristics are simplified
to angular momentum, employing a geodesic-like equation for motion
description. However, by reducing photon dynamics to geodesic-like equations,
this method simplifies the spatial distribution of wave packets into a mere
representation of angular momentum. Consequently, it fails to adequately
capture the intrinsic OAM of vortex photons, which is significantly influenced
by the spatial distribution of their wave packets. This discussion suggests
that the MPD equations may have limitations in accurately describing the
motion of vortex photons, possessing intrinsic OAM. Our investigation into the
gravitational OHE for vortex photons highlights these limitations.

##  IV Discussion

In this study, we investigate the dynamics of a vortex photon modeled as a
vortex Laguerre-Gaussian electromagnetic wave packet within the Lense-Thirring
metric by solving Maxwell’s equations. The trajectory of the photon is
determined by analyzing the center of its EMT. Distinct from the gravitational
SHE, where photons with opposite spins predominantly diverge perpendicular to
the null geodesic plane, vortex photons with different intrinsic OAM also
demonstrate separation within this plane. This separation, driven by intrinsic
OAM, is identified as the gravitational orbital Hall effect (OHE). This
differentiation between the gravitational spin and orbital Hall effects
suggests that, in gravitational fields, the intrinsic OAM of photons may not
align with their spin angular momentum, highlighting a complex coupling
between photon angular momentum and gravity.

For a spin-polarized vortex photon moving from a finite distance $r_{0}$ to
infinity, its deviation angle $\theta_{\sigma+\ell}$ from the null geodesic
plane can be effectively approximated by the expression
$\theta_{\sigma+\ell}\simeq(\alpha\sigma+\beta\ell)GM\lambda/2\pi r_{0}^{2}$.
Here, $r_{0}$ denotes the photon’s initial distance from the gravitational
source. The symbols $\sigma$ and $\ell$ represent the SAM and intrinsic OAM of
the vortex photon, respectively. The coefficients $\alpha$ and $\beta$ depend
on the chosen expression and component of the EMT that are used to define the
center of the photon. Specifically, when the energy density of the symmetric
EMT is applied, these coefficients are established as $\alpha=2$ and
$\beta=-1$. This leads to a specific expression for the deviation angle:
$\theta_{\sigma+\ell}\simeq(2\sigma-\ell)GM\lambda/2\pi r_{0}^{2}$.

The expression $\theta_{\sigma+\ell}\simeq(2\sigma-\ell)GM\lambda/2\pi
r_{0}^{2}$ implies that two photons with equal total angular momentum
$\sigma+\ell=0$, when propagating through a gravitational field, may diverge
from one another. Specifically, a photon with SAM $\sigma=0$ and intrinsic OAM
$\ell=0$ is predicted to follow the null geodesic, resulting in a vanishing
deviation angle $\theta_{0+0}$. Conversely, a photon with SAM $\sigma=1$ and
intrinsic OAM $\ell=-1$ is expected to deviate from the null geodesic,
exhibiting a deviation angle $\theta_{1+(-1)}\simeq 3GM\lambda/2\pi
r_{0}^{2}$. Thus, despite possessing equal total angular momentum
$\sigma+\ell$, the trajectories of these two photons may differ, demonstrating
a non-zero separation angle.

The intrinsic OAM of a vortex photon is characterized as a phase factor
$\psi\propto\exp(i\ell\phi)$ in its electromagnetic wave packet, with
$\ell\phi$ indicating the dependence on the spatial distribution of its wave
packet. Our findings indicate that, when coupled to gravity, the intrinsic OAM
of photons might not align with their SAM. This discrepancy suggests that the
spatial distribution of the vortex photon’s wave packet extends beyond a
simple angular momentum representation, necessitating a more detailed
examination of its wave packet characteristics.

## Acknowledgement

We are grateful to Bei Liu, Li-Li Yang and Qi-Liang Zhao for the helpful
discussions, and we thank Peng-Cheng Zhao for his suggestions during the
numerical calculations. This work was supported by the Fundamental Research
Funds for the Central Universities, Sun Yat-sen University.

##  Appendix A Equations of motion for vortex photons in Lense-Thirring metric

Within the weak field approximation, characterized by $GM/r\ll 1$, the
dynamics of the perturbation term $\widetilde{A}^{\nu}$ for vortex photons are
governed by Eq. (7). By adhering to the first order of $x^{i}/b$ as stipulated
by Eqs. (II.1), Eq. (7) simplifies to

$\displaystyle\left(-\partial_{x}^{2}+\nabla^{2}\right)\tilde{A}^{0}$
$\displaystyle=\frac{2GM}{b^{3}}\left(z\partial_{t}\bar{A}^{3}+(b-2x)\partial_{t}\bar{A}^{1}+y\partial_{t}\bar{A}^{2}\right)+\frac{2GMa}{b^{5}}(-6by(\partial_{x}\bar{A}^{1}-\partial_{y}\bar{A}^{2}))$
$\displaystyle+\frac{2GMa}{b^{5}}(3(b^{2}-3bx)(\partial_{y}\bar{A}^{1}+\partial_{x}\bar{A}^{2})+3bz(\partial_{z}\bar{A}^{2}+\partial_{y}\bar{A}^{3})),$
(23) $\displaystyle\left(-\partial^{2}_{t}+\nabla^{2}\right)\tilde{A}^{1}$
$\displaystyle=\frac{2GM}{b}\left(\partial^{2}_{t}\bar{A}^{1}+\partial_{i}\partial^{i}\bar{A}^{1}\right)+\frac{2GM}{b^{3}}(2(b-2x)\partial_{x}\bar{A}^{1}+y\left(\partial_{y}\bar{A}^{1}+\partial_{x}\bar{A}^{2}\right))$
$\displaystyle+\frac{2GMz}{b^{3}}\left(\partial_{z}\bar{A}^{1}+\partial_{x}\bar{A}^{3}\right)+\frac{4aGM}{b^{3}}\left(-y\partial_{t,x}\bar{A}^{1}+(b-2x)\partial_{t,y}\bar{A}^{1}\right)$
$\displaystyle+\frac{2aGM}{b^{3}}\left(1-\frac{3ax}{b}\right)\partial_{t}\bar{A}^{2},$
(24) $\displaystyle\left(-\partial_{t}^{2}+\nabla^{2}\right)\tilde{A}^{2}$
$\displaystyle=\frac{2GM}{b}\left(\partial_{t}^{2}\bar{A}^{2}+\partial_{i}\partial^{i}\bar{A}^{2}\right)+\frac{2GM}{b^{3}}(2y\partial_{y}\bar{A}^{2}+(b-2x)(\partial_{x}\bar{A}^{2}+\partial_{y}\bar{A}^{1}))$
$\displaystyle+\frac{2GMz}{b^{3}}(\partial_{z}\bar{A}^{2}+\partial_{y}\bar{A}^{3})+\frac{2GMa}{b^{3}}((1-\frac{3x}{b})\partial_{t}\bar{A}^{1}+\frac{3z}{b}\partial_{t}\bar{A}^{3})$
$\displaystyle+\frac{4GMa}{b^{3}}\left(-y\partial_{t,x}\bar{A}^{2}+(b-2x)\partial_{t,y}\bar{A}^{2}\right),$
(25) $\displaystyle(-\partial_{t}^{2}+\nabla^{2})\tilde{A}^{3}$
$\displaystyle=\frac{2GM}{b}(\partial_{t}^{2}\bar{A}^{3}+\partial_{i}\partial^{i}\bar{A}^{3})+\frac{2GM}{b^{3}}(2z\partial_{z}\bar{A}^{3}+(b-2x)(\partial_{y}\bar{A}^{1}+\partial_{x}\bar{A}^{2}))$
$\displaystyle+\frac{2GMy}{b^{3}}(\partial_{z}\bar{A}^{2}+\partial_{y}\bar{A}^{3})+\frac{6GMa}{b^{5}}(y\partial_{t}\bar{A}^{1}-(b-4x)\partial_{t}\bar{A}^{2})$
$\displaystyle+\frac{4GMa}{b^{3}}(-y\partial_{t,x}\bar{A}^{3}+(b-2x)\partial_{t,y}\bar{A}^{3}),$
(26)

where $i=1,2,3$ denotes the spatial indices.

Given the complexity of solving second-order inhomogeneous partial
differential equations, especially due to their challenging nature even with
numerical methods, we opt for a three-dimensional spatial Fourier transform to
facilitate their analysis. Hence, Eq. (A) becomes to

$\displaystyle\left(\partial_{t}^{2}+\omega^{2}\right)\widetilde{A}_{f}^{0}$
$\displaystyle=\frac{2GM}{b^{2}}\left(\left(1-\frac{2i}{b}\frac{\partial}{\partial
k_{x}}\right)\left(i\omega\bar{A}_{f}^{1}\right)+\frac{i}{b}\frac{\partial}{\partial
k_{y}}\left(i\omega\bar{A}_{f}^{2}\right)+\frac{i}{b}\frac{\partial}{\partial
k_{z}}\left(i\omega\bar{A}_{f}^{3}\right)\right)$
$\displaystyle+\frac{6GMa}{b^{3}}\left(\frac{i}{b}\frac{\partial}{\partial
k_{y}}(ik_{x}\bar{A}_{f}^{1}-ik_{y}\bar{A}_{f}^{2})-\left(1-\frac{3i}{b}\frac{\partial}{\partial
k_{x}}\right)(ik_{x}\bar{A}_{f}^{2}+ik_{y}\bar{A}_{f}^{1})\right)$
$\displaystyle-\frac{6GMa}{b^{3}}\left(\frac{i}{b}\frac{\partial}{\partial
k_{z}}\left(ik_{z}\bar{A}_{f}^{2}+ik_{y}\bar{A}_{f}^{3}\right)\right)$ (27)
$\displaystyle\left(\partial_{t}^{2}+\omega^{2}\right)\widetilde{A}_{f}^{1}$
$\displaystyle=\frac{4GM}{b}\left(1-\frac{i}{b}\frac{\partial}{\partial
k_{x}}\right)(\omega^{2}\bar{A}_{f}^{1})+\frac{2GM}{b^{2}}\left(\frac{\partial(k_{z}\bar{A}_{f}^{1}+k_{x}\bar{A}_{f}^{3})}{b\partial
k_{z}}\right)-\frac{4iGM}{b^{2}}$
$\displaystyle\times\left(\left(1-\frac{2i}{b}\frac{\partial}{\partial
k_{x}}\right)(k_{x}\bar{A}_{f}^{1})+\frac{i\partial(k_{x}\bar{A}_{f}^{2}+k_{y}\bar{A}_{f}^{1})}{2b\partial
k_{y}}\right)-\frac{2iGMa}{b^{3}}\left(1-\frac{3i}{b}\frac{\partial}{\partial
k_{x}}\right)$
$\displaystyle\times(\omega\bar{A}_{f}^{2})+\frac{4GMa}{b^{3}}\left(\frac{i\partial}{\partial
k_{y}}(\omega
k_{x}\bar{A}_{f}^{1})\right)-\frac{4GMa}{b^{2}}\left(1-\frac{2i\partial}{b\partial
k_{x}}(\omega k_{y}\bar{A}_{f}^{1})\right)$ (28)
$\displaystyle\left(\partial_{t}^{2}+\omega^{2}\right)\widetilde{A}_{f}^{2}$
$\displaystyle=\frac{4GM}{b}\left(1-\frac{i}{b}\frac{\partial}{\partial
k_{x}}\right)\left(\omega^{2}\bar{A}_{f}^{2}\right)-\frac{2GM}{b^{2}}\left(\left(1-2\frac{i}{b}\frac{\partial}{\partial
k_{x}}\right)(ik_{y}\bar{A}_{f}^{1}+ik_{x}\bar{A}_{f}^{2})\right)$
$\displaystyle-\frac{2GM}{b^{2}}\left(\frac{2\partial(ik_{y}\bar{A}_{f}^{2})}{b\partial
k_{y}}-\frac{\partial(k_{z}\bar{A}_{f}^{2}+k_{y}\bar{A}_{f}^{3})}{b\partial
k_{z}}\right)+\frac{2iGMa}{b^{3}}\left(1-\frac{3i}{b}\frac{\partial}{\partial
k_{x}}\right)(\omega\bar{A}_{f}^{1})$
$\displaystyle-\frac{2GMa}{b^{4}}\frac{\partial(\omega\bar{A}_{f}^{3})}{\partial
k_{z}}-\frac{4GMa}{b^{2}}\left(-\frac{i\partial(\omega
k_{x}\bar{A}_{f}^{2})}{\partial
k_{y}}+\left(1-\frac{2i}{b}\frac{\partial}{\partial k_{x}}\right)(\omega
k_{y}\bar{A}_{f}^{2})\right)$ (29)
$\displaystyle\left(\partial_{t}^{2}+\omega^{2}\right)\widetilde{A}_{f}^{3}$
$\displaystyle=\frac{4GM}{b}\left(1-\frac{i}{b}\frac{\partial}{\partial
k_{x}}\right)(\omega^{2}\bar{A}_{f}^{3})-\frac{2iGM}{b^{2}}\left(1-\frac{2i}{b}\frac{\partial}{\partial
k_{x}}\right)(k_{z}\bar{A}_{f}^{1}+k_{x}\bar{A}_{f}^{3})$
$\displaystyle+\frac{2GM}{b^{3}}\left(\frac{\partial(k_{y}\bar{A}_{f}^{3}+k_{z}\bar{A}_{f}^{2})}{\partial
k_{y}}+\frac{2\partial(k_{z}\bar{A}_{f}^{3})}{\partial
k_{z}}\right)+\frac{6GMa}{b^{4}}\left(\frac{\partial(\omega\bar{A}_{f}^{2})}{\partial
k_{z}}\right)$
$\displaystyle-\frac{4GMa}{b^{2}}\left(\left(1-\frac{2i}{b}\frac{\partial}{\partial
k_{x}}\right)(\omega k_{y}\bar{A}_{f}^{3})-\frac{i}{b}\frac{\partial}{\partial
k_{y}}(\omega k_{x}\bar{A}_{f}^{3})\right),$ (30)

where $\bar{A}^{\rho}_{f}$ and $\widetilde{A}^{\rho}_{f}$ are defined as

$\bar{A}^{\rho}_{f}=\frac{1}{2\pi^{3}}\int\bar{A}^{\rho}\exp(-i\vec{k}\cdot\vec{x})\text{d}^{3}x,\quad\widetilde{A}^{\rho}_{f}=\frac{1}{2\pi^{3}}\int\widetilde{A}^{\rho}\exp(-i\vec{k}\cdot\vec{x})\text{d}^{3}x.\quad$
(31)

## Reference

## References

  * Bliokh (2006) K. Y. Bliokh, Physical Review Letters 97, 43901 (2006).
  * Bliokh  _et al._ (2008) K. Y. Bliokh, A. Niv, V. Kleiner, and E. Hasman, Nature Photonics 2, 748 (2008).
  * Bliokh and Desyatnikov (2009) K. Y. Bliokh and A. S. Desyatnikov, Physical Review A 79, 011807 (2009).
  * Zhou  _et al._ (2013) X. Zhou, J. Zhang, X. Ling, S. Chen, H. Luo, and S. Wen, Physical Review A 88, 053840 (2013).
  * Liu  _et al._ (2015) Y. Liu, X. Ling, X. Yi, X. Zhou, S. Chen, Y. Ke, H. Luo, and S. Wen, Optics letters 40, 756 (2015).
  * Bliokh  _et al._ (2015) K. Y. Bliokh, F. J. Rodríguez-Fortuño, F. Nori, and A. V. Zayats, Nature Photonics 9, 796 (2015).
  * Ling  _et al._ (2023) X. Ling, Z. Zhang, Z. Dai, Z. Wang, H. Luo, and L. Zhou, Laser & Photonics Reviews 17, 2200783 (2023).
  * Fu  _et al._ (2019) S. Fu, C. Guo, G. Liu, Y. Li, H. Yin, Z. Li, and Z. Chen, Physical Review Letters 123, 243904 (2019).
  * Mieling and Oancea (2023) T. B. Mieling and M. A. Oancea, Physical Review Research 5, 023140 (2023).
  * Ling  _et al._ (2017) X. Ling, X. Zhou, K. Huang, Y. Liu, C.-W. Qiu, H. Luo, and S. Wen, Reports on Progress in Physics 80, 066401 (2017).
  * Liu  _et al._ (2022) S. Liu, S. Chen, S. Wen, and H. Luo, Opto-Electronic Science 1, 220007 (2022).
  * Sheng  _et al._ (2023) L. Sheng, Y. Chen, S. Yuan, X. Liu, Z. Zhang, H. Jing, L.-M. Kuang, and X. Zhou, Progress in Quantum Electronics , 100484 (2023).
  * Papapetrou (1951) A. Papapetrou, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 209, 248 (1951).
  * Corinaldesi and Papapetrou (1951) E. Corinaldesi and A. Papapetrou, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 209, 259 (1951).
  * Dixon (1970a) W. G. Dixon, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 314, 499 (1970a).
  * Dixon (1970b) W. G. Dixon, Proc. Roy. Soc. Lond. A 319, 509 (1970b).
  * Dixon (1974) W. G. Dixon, Philosophical Transactions of the Royal Society of London Series A 277, 59 (1974).
  * Dixon (2015) W. Dixon, Equations of Motion in Relativistic Gravity , 1 (2015).
  * Carmeli  _et al._ (1977) M. Carmeli, C. Charach, and M. Kaye, Phys. Rev. D 15, 1501 (1977).
  * Plyatsko  _et al._ (2011) R. Plyatsko, O. Stefanyshyn, and M. Fenyk, Class. Quant. Grav. 28, 195025 (2011), arXiv:1110.1967 [gr-qc] .
  * Hackmann  _et al._ (2014) E. Hackmann, C. Lämmerzahl, Y. N. Obukhov, D. Puetzfeld, and I. Schaffer, Phys. Rev. D 90, 064035 (2014), arXiv:1408.1773 [gr-qc] .
  * Antoniou  _et al._ (2019) I. Antoniou, D. Papadopoulos, and L. Perivolaropoulos, Class. Quant. Grav. 36, 085002 (2019), arXiv:1903.03835 [gr-qc] .
  * Duval  _et al._ (2019) C. Duval, L. Marsot, and T. Schücker, Physical Review D 99, 124037 (2019).
  * Zhang  _et al._ (2020) Y.-P. Zhang, S.-W. Wei, and Y.-X. Liu, Universe 6, 103 (2020), arXiv:2003.10960 [gr-qc] .
  * Costa and Natário (2015) L. F. O. Costa and J. Natário, Fund. Theor. Phys. 179, 215 (2015), arXiv:1410.6443 [gr-qc] .
  * Frolov and Shoom (2012) V. P. Frolov and A. A. Shoom, Phys. Rev. D 86, 024010 (2012).
  * Oancea  _et al._ (2020) M. A. Oancea, J. Joudioux, I. Dodin, D. Ruiz, C. F. Paganini, and L. Andersson, Physical Review D 102, 024075 (2020).
  * Deriglazov (2021) A. A. Deriglazov, Phys. Rev. D 104, 025006 (2021).
  * Harte and Oancea (2022) A. I. Harte and M. A. Oancea, Phys. Rev. D 105, 104061 (2022).
  * Andersson and Oancea (2023) L. Andersson and M. A. Oancea, Class. Quantum Grav. 40, 154002 (2023).
  * Oancea and Harko (2024) M. A. Oancea and T. Harko, Physical Review D 109, 064020 (2024).
  * Lian and Chen (2022) D.-D. Lian and X.-S. Chen, Physical Review D 105, 104008 (2022).
  * Lian and Zhang (2023) D.-D. Lian and P.-M. Zhang, arXiv preprint arXiv:2312.14391 (2023).
  * Obukhov (2001) Y. N. Obukhov, Phys. Rev. Lett. 86, 192 (2001), arXiv:gr-qc/0012102 .
  * Silenko and Teryaev (2005) A. J. Silenko and O. V. Teryaev, Phys. Rev. D 71, 064016 (2005), arXiv:gr-qc/0407015 .
  * Obukhov  _et al._ (2009) Y. N. Obukhov, A. J. Silenko, and O. V. Teryaev, Phys. Rev. D 80, 064044 (2009), arXiv:0907.4367 [gr-qc] .
  * Obukhov  _et al._ (2013) Y. N. Obukhov, A. J. Silenko, and O. V. Teryaev, Phys. Rev. D 88, 084014 (2013), arXiv:1308.4552 [gr-qc] .
  * Wang (2024) Z.-L. Wang, Phys. Rev. D 109, 044060 (2024).
  * Landau and Lifshitz (2013) L. D. Landau and E. M. Lifshitz,  _Quantum mechanics: non-relativistic theory_ , Vol. 3 (Elsevier, 2013).
  * Cohen and Tannoudji (1977) C. Cohen and B. Tannoudji, Diu, and f. laloe, quantum mechanics, vol. i (1977).
  * Beijersbergen  _et al._ (1993) M. W. Beijersbergen, L. Allen, H. Van der Veen, and J. Woerdman, Optics Communications 96, 123 (1993).
  * Bliokh  _et al._ (2017) K. Y. Bliokh, I. P. Ivanov, G. Guzzinati, L. Clark, R. Van Boxem, A. Béché, R. Juchtmans, M. A. Alonso, P. Schattschneider, F. Nori, _et al._ , Physics Reports 690, 1 (2017).
  * Gahagan and Swartzlander (1996) K. Gahagan and G. Swartzlander, Optics Letters 21, 827 (1996).
  * Garcés-Chávez  _et al._ (2003) V. Garcés-Chávez, D. McGloin, M. Padgett, W. Dultz, H. Schmitzer, and K. Dholakia, Physical review letters 91, 093602 (2003).
  * Mair  _et al._ (2001) A. Mair, A. Vaziri, G. Weihs, and A. Zeilinger, Nature 412, 313 (2001).
  * Leach  _et al._ (2010) J. Leach, B. Jack, J. Romero, A. K. Jha, A. M. Yao, S. Franke-Arnold, D. G. Ireland, R. W. Boyd, S. M. Barnett, and M. J. Padgett, Science 329, 662 (2010).
  * Foo  _et al._ (2005) G. Foo, D. M. Palacios, and G. A. Swartzlander, Optics letters 30, 3308 (2005).
  * Harwit (2003) M. Harwit, The Astrophysical Journal 597, 1266 (2003).
  * Kruglov and Vlasov (1985) V. Kruglov and R. Vlasov, Physics Letters A 111, 401 (1985).
  * Swartzlander Jr and Law (1992) G. Swartzlander Jr and C. Law, Physical Review Letters 69, 2503 (1992).
  * Lai  _et al._ (2018) M.-Y. Lai, Y.-L. Wang, G.-H. Liang, F. Wang, and H.-S. Zong, Physical Review A 97, 033843 (2018).
  * Choi  _et al._ (2023) Y.-G. Choi, D. Jo, K.-H. Ko, D. Go, K.-H. Kim, H. G. Park, C. Kim, B.-C. Min, G.-M. Choi, and H.-W. Lee, Nature 619, 52 (2023).
  * Adler (2015) R. J. Adler, Classical and Quantum Gravity 32, 224002 (2015).
  * Wald (2010) R. M. Wald,  _General relativity_ (University of Chicago press, 2010).
  * Liu and Ivanov (2023) B. Liu and I. P. Ivanov, Physical Review A 107, 063110 (2023).