11institutetext: Mohamed bin Zayed University of Artificial Intelligence,
Abu Dhabi, United Arab Emirates 22institutetext: Abu Dhabi Health Services
Company (SEHA),
United Arab Emirates
22email<EMAIL_ADDRESS>
# Multi-Task Learning Approach for Unified Biometric Estimation from Fetal
Ultrasound Anomaly Scans

Mohammad Areeb Qazi 11 Mohammed Talha Alam 11 Ibrahim Almakky 11 Werner
Gerhard Diehl 22 Leanne Bricker 22 Mohammad Yaqub 11

###### Abstract

Precise estimation of fetal biometry parameters from ultrasound images is
vital for evaluating fetal growth, monitoring health, and identifying
potential complications reliably. However, the automated computerized
segmentation of the fetal head, abdomen, and femur from ultrasound images,
along with the subsequent measurement of fetal biometrics, remains
challenging. In this work, we propose a multi-task learning approach to
classify the region into head, abdomen and femur as well as estimate the
associated parameters. We were able to achieve a mean absolute error (MAE) of
1.08 mm on head circumference, 1.44 mm on abdomen circumference and 1.10 mm on
femur length with a classification accuracy of 99.91% on a dataset of fetal
Ultrasound images. To achieve this, we leverage a weighted joint
classification and segmentation loss function to train a U-Net architecture
with an added classification head. The code can be accessed through Github.

###### Keywords:

Ultrasound Images Segmentation Classification Fetal Biometry.

##  1 Introduction

Ultrasound is a widely used non-invasive imaging modality in obstetrics for
monitoring pregnancy progression and assessing fetal development. The term
“gestational age” (GA), measured in weeks, is commonly used during pregnancy
to indicate the stage of pregnancy. Accurately estimating gestational age and
fetal biometry parameters through ultrasound images plays a critical role in
evaluating fetal growth progress and identifying any concerns during
development [1]. However, manual measurement of biometry parameters from
ultrasound images is a time-consuming and operator-dependent task, leading to
inconsistent and inaccurate measurements that can affect the interpretation of
fetal growth and development [2].

Deep learning methods have demonstrated impressive performance in various
medical imaging tasks, such as segmentation, registration, and classification
[3]. Despite these advancements, there remains a need to leverage deep
learning techniques for the accurate estimation of gestational age and fetal
biometry parameters from ultrasound images. By automating the process, deep
learning models can provide consistent and reliable measurements, overcoming
the limitations associated with manual estimation. Furthermore, the
integration of deep learning algorithms holds promise in enabling early
detection of fetal growth abnormalities, thereby improving prenatal care and
the management of potential complications [4]. The automation of biometry
estimation in ultrasound imaging has emerged as a compelling area of research
with the potential to enhance the efficiency of prenatal care. Clinicians
recognize the significance of accurately estimating biometric parameters,
including head circumference (HC), femur length (FL), and abdominal
circumference (AC), in assessing fetal growth and development. As a result,
there has been a growing interest in leveraging machine learning methods to
automate biometry estimation, primarily focusing on HC estimation.

The current body of literature on biometry parameter estimation primarily
focuses on HC estimation, with limited attention given to other essential
biometric measurements such as FL and AC. Consequently, there exists a
research gap regarding the comprehensive assessment of multiple biometric
parameters in ultrasound analysis [5]. In this paper, we aim to fill this gap
by developing an end-to-end network that incorporates the estimation of
multiple biometric parameters. Instead of solely emphasizing HC estimation,
our framework will incorporate the calculation of FL and AC measurements. This
integrated approach, encompassing anatomical region identification and
subsequent biometric parameter estimation, has the potential to be a valuable
tool for clinicians in their assessment of fetal development.

In this work, our main contributions are:

  1. 1.

We introduce a comprehensive multi-task framework that combines classification
and accurate estimation of biometric parameters, namely, HC, AC and FL.

  2. 2.

We conduct an investigation into determining the optimal weight allocations
for the losses within our multi-task setup. This analysis sheds light on the
crucial interplay between losses, enhancing our understanding of the system’s
overall performance.

##  2 Methods

Figure 1: (a) Overview of the proposed multi-task learning network. B: Brain,
A: Abdomen, F: Femur.

This research introduces a multitask learning method for 2D ultrasound images.
The proposed model proposes an end-to-end Convolutional Neural Network (CNN)
architecture that does classification as well as segmentation. The model
predicts the class and the segmentation mask which are forwarded for post-
processing, producing the bio-parameter associated with that class. The
proposed framework is shown in Fig. 1. We use U-Net [6] as the backbone
network for the task due to its excellent performance in 2D medical image
segmentation.

The U-Net design can be separated into three parts: an encoding path, a
decoding path, and skip connections between them. U-Net incorporates stacked
convolution layers in each level, as shown in Fig. 1. In our approach, we
employed five downsampling steps in the encoder to capture sufficient
features. Correspondingly, the decoder used five upsampling steps to restore
the features to the original size, leveraging the encoded information to
predict the target’s segmentation mask. Skip connections linked feature maps
from the encoder to the decoder, allowing spatial information to flow and
improving segmentation results. All convolution layers used 3x3 kernels,
followed by batch normalization and ReLU activation. Additionally, we utilized
2x2 pooling to downsample feature maps effectively.

###  2.1 Multi-task Learning Network

Taking inspiration from the established practice in image classification,
where prominent CNN models like VGG [7] and ResNet [8] incorporate high-level
feature maps, our approach capitalizes on this insight. We extend the idea to
our proposed multi-task learning network, enabling the extraction of shared
features from U-Net for segmentation and classification. Within this network
architecture, we introduce a classification branch at the base of the U-Net,
as depicted in Fig. 1. This branch utilizes feature maps from the bottleneck
layer, which are flattened and directed through a fully connected (FC) layer,
followed by a softmax layer. This ensemble empowers the network to predict the
input image class into brain, abdomen, or femur regions, which also helps in
the segmentation task.

For the classification task, we employed the cross-entropy loss function. This
loss function quantifies the dissimilarity between predicted and actual class
probabilities, aiming to minimize the discrepancy during training, and it is
formulated as follows:

$\mathcal{L}_{\mathrm{cls}}(x,y)=-\frac{1}{N}\sum_{n=1}^{N}\log\left(\frac{\exp\left(x_{n,y_{n}}\right)}{\sum_{c=1}^{C}\exp\left(x_{n,c}\right)}\right)$
(1)

where $x$ represents the input, $y$ signifies the ground truth class, $C$ is
the total number of classes, and $N$ is the batch size.

For the segmentation task, we employed the Dice loss function, which is
formulated as follows:

$\mathcal{L}_{\mathrm{seg}}(p,q)=1-\frac{2\sum_{i}p(i)\cdot
q(i)}{\sum_{i}p(i)^{2}+\sum_{i}q(i)^{2}}$ (2)

where the variable $p$ signifies the ground truth binary mask, whereas $q$
stands for the predicted segmentation mask. The index $i$ pertains to the
individual elements within the mask. Within our methodology, we integrate the
classification loss and the segmentation loss through a linear combination,
regulated by a hyperparameter denoted as $\lambda$. This combined loss, known
as the multi-task loss, is precisely defined as:

$\mathcal{L}_{\mathrm{joint}}=\lambda\cdot\mathcal{L}_{\mathrm{cls}}+(1-\lambda)\cdot\mathcal{L}_{\mathrm{seg}}$
(3)

The significance of $\lambda$ lies in its capacity to assign weights to
individual losses. This strategic weight allocation hinges on the concept that
tasks characterized by higher training losses merit increased attention during
the optimization process. By systematically varying $\lambda$ and conducting
an ablation study, we gain insights into its impact and determine the optimal
value that aligns with our objectives.

###  2.2 Estimating the Bio-Parameters

To calculate the bio parameter, we post-process the segmentation masks. We
observed that the estimation process of HC and AC differs from the FL
estimation process. This distinction arises due to the distinct shape
characteristics of these anatomical regions. The head and abdomen regions can
be assumed to have an elliptical shape, allowing us to apply ellipse-fitting
techniques to extract the border of the head/abdomen from the segmented
regions. Then the circumference can be approximated using the following
equation:

$\displaystyle
circumference=\pi\frac{a+b+3-(a-b)^{2}}{10a+10b+\sqrt{a^{2}+14ab+b^{2}}}$ (4)

where $a,b,x,y$ is the semi-minor and semi-major axis lengths and the centre
coordinate of the ellipse.

In contrast to the elliptical shape of the brain and abdomen regions, we found
the femur region to have a rectangular shape. To estimate the dimensions of
the femur, we apply a rectangle fitting technique, allowing us to extract the
length and breadth of the rectangular region. Subsequently, we calculate the
perimeter of the rectangle using the obtained length and breadth measurements.

##  3 Data and Experiments

###  3.1 Dataset

Our dataset consists of images collected from a hospital over one year. The
dataset comprises a total of 5,872 images from 2206 subjects, with 2,023
images depicting the brain, 2,128 images representing the abdomen, and the
remaining 1,721 images focusing on the femur. To ensure comprehensive
evaluation, we divided the dataset on the basis of subjects into an 80%
training set and a 20% testing set. Additionally, we allocated 10% of the
training set as a separate validation set for fine-tuning and performance
monitoring.

The dataset included images and some associated information. Specifically, for
the brain and abdomen, the dataset contained coordinates for the center,
length of the the major and minor axes with an angle with the x-axis. Using
these values, elliptical masks were created to represent these regions in the
images. For the femur, the data had endpoints of it. These were connected to
create a line with a fixed width, essentially forming a mask for the femur
class. To make the data consistent, both the images and the masks were resized
to a common size of 256 x 256 before using them for the model.

###  3.2 Experiments

To evaluate the effectiveness of the proposed method, we conducted ablation
experiments. The results for each of these experiments are shown in Tab. 1. We
started by training a separate U-Net model for each of the three classes which
serves as the baseline and later helped us understand the benefit of training
a single model for all three classes. Then, we added the classification head
in the model and trained only the classification head by keeping the weight of
the dice loss 0. Similarly, we also trained the model by keeping the
classification loss weight as 0. Then we trained different models by tweaking
the values of $\lambda$.

All models were implemented using Pytorch [9] and were trained and tested on
Quadro RTX 6000 GPUs (25G) using the ADAMax optimizer with a decaying learning
rate initialized at 5e-4. The batch size and the number of epochs were set to
16 and 100 for each experiment respectively.

In this context, considering the relatively manageable nature of the
classification task and the balanced dataset, we opted to utilize only
accuracy as the metric for evaluating the classification component of the
model. For the segmentation task, we used the biometric parameter error metric
i.e. Mean Average Error (MAE). We apply the post-processing steps on the mask
and then calculate MAE. MAE gauges the average magnitude of discrepancies
between predicted and true values. The MAE values served as objective metrics
for quantifying the results of our models in capturing these measurements. The
formula for MAE is:

$MAE=\frac{1}{n}\sum_{i=1}^{n}|y_{\text{true}}^{(i)}-y_{\text{pred}}^{(i)}|$
(5)

Here, $n$ stands for the number of samples, $y_{\text{true}}^{(i)}$ represents
the true biometric parameter value for the ith sample, and
$y_{\text{pred}}^{(i)}$ is the corresponding predicted value.

Table 1: The effects of $\lambda$ on Multi-Task Learning: The table presents
an analysis of varying $\lambda$ values on the Multi-Task Learning framework.
The performance metrics encompass accuracy percentages and MAE for distinct
organ classes. As $\lambda$ is adjusted, notable trends emerge in both
accuracy and estimation errors across the brain, abdomen, and femur classes.

$\lambda$ | Accuracy (%) | Brain (MAE in mm) | Abdomen (MAE in mm) | Femur (MAE in mm)
---|---|---|---|---
1 | 100 | - | - | -
0.8 | 100 | 1.50 ±1.6156 | 2.33 ±4.9390 | 1.50 ±2.1424
0.6 | 100 | 1.88 ±1.7502 | 2.78 ±2.9089 | 1.65 ±2.6943
0.4 | 99.91 | 1.39 ±1.2199 | 2.08 ±3.9290 | 1.53 ±6.3734
0.2 | 100 | 1.41 ±1.2420 | 1.83 ±2.0092 | 1.13 ±1.6356
0.1 | 100 | 1.19 ±1.0825 | 1.55 ±1.7488 | 1.14 ±1.5085
0.05 | 99.91 | 1.46 ±1.2163 | 2.11 ±3.7542 | 1.14 ±1.7677
0.025 | 100 | 1.24 ±1.1018 | 1.73 ±2.7895 | 1.31 ±1.9683
0.01 | 100 | 1.22 ±1.1179 | 1.61 ±1.8859 | 1.19 ±1.4767
0.001 | 99.91 | 1.08 ±0.9702 | 1.44 ±1.9269 | 1.13 ±1.4498
0.00001 | 99.33 | 1.18 ±1.0052 | 1.54 ±1.7283 | 1.10 ±1.1298
0 | - | 11.24 ±23.1575 | 14.94 ±26.8828 | 65.12 ±56.0734
Separate Models | 1.71 ±1.35 | 1.55 ±1.88 | 1.08 ±1.18

##  4 Results and Discussion

Tab. 1 provides an overview of our ablation experiment outcomes. Initially,
when we used separate models for each class, we obtained MAE values of 1.71 mm
for the brain, 1.55 mm for the abdomen, and 1.08 mm for the femur. We then
trained only the encoder for classification across all classes (see Tab. 1,
$\lambda$ = 1), successfully classifying all test images accurately. Moreover,
we also trained the U-Net without the classification head (see Tab. 1,
$\lambda$ = 0) to examine its effect on segmentation when classification is
excluded. The results were notably worse, yielding 11.24 mm MAE for the brain,
14.95 mm MAE for the abdomen, and 65.12 mm MAE for the femur. These outcomes
are considered the baseline for the ablation study, guiding us in finding the
optimal balance between the weights of the two losses.

The value of $\lambda$ in Eq. 3 affects the trade-off between classification
and segmentation performance in multi-task learning. We experimented with
various $\lambda$ values and observed interesting trends in Tab. 1. Larger
$\lambda$ values gave excellent classification accuracy but had some impact on
segmentation. Lowering $\lambda$ improved segmentation without significantly
compromising classification. Surprisingly, very low $\lambda$ values produced
good segmentation. Among all tested values, $\lambda$ = 0.001 struck the best
balance. Although classification performance was consistently strong across
values near 1 and smaller, the MAE saw the most improvement from 1 to 0.001.
Further, lowering $\lambda$ further led to diminishing returns. Fig. 2
illustrates this relationship graphically, showing how reducing $\lambda$
enhanced MAEs.

Figure 2: Effect of $\lambda$ on the Bio-Parameters: The image illustrates the
impact of varying $\lambda$ values on the bio-parameters. Noticeably, the
model’s performance is better with lower $\lambda$ values.

Fig. 3 visually presents the impact of varying $\lambda$ on the convergence of
different losses. In Fig. 3(a), we observe that a lower $\lambda$ value led to
quicker convergence of the dice loss. This can be attributed to the model’s
increased emphasis on the dice loss, which in turn accelerated the learning
process. Conversely,higher $\lambda$ values caused greater fluctuations in the
cross-entropy loss. This behavior suggests that the model found learning the
classes relatively easier when the $\lambda$ was low. Consequently, the
overall loss, as depicted in Fig. 3(b), exhibited smoother convergence when
$\lambda$ was set to a smaller value. The optimal $\lambda$ value is 0.001 at
which learning was notably seamless for both tasks, and the results reached
their peak performance.

The pronounced improvement observed with lower values of $\lambda$ concludes
that the classification task is relatively more manageable. This is attributed
to the distinct structural characteristics of the classes, allowing the model
to effectively grasp features and distinguish them. The multi-task learning
framework plays a pivotal role in enhancing performance across individual
classes, leveraging the increased volume of training data to capture a broader
spectrum of variations.

(a) Dice & Cross-Entropy Loss w.r.t $\lambda$

(b) Total Loss w.r.t $\lambda$

Figure 3: Training Loss Graphs with various values of $\lambda$.

The integration of a classification head yielded substantial advantages,
particularly evident in improved segmentation outcomes and more precise
estimations of the bio-parameter. This enhancement was prominently visible in
the error reductions for specific classes: the MAE for the brain class
diminished from 1.71 mm to 1.08 mm, and for the abdomen class, it decreased
from 1.55 mm to 1.44 mm. Interestingly, the femur class exhibited relatively
consistent performance, similar to training a standalone model. This could be
attributed to the structural similarity between the brain and abdomen classes,
allowing the model to leverage knowledge gained from one class to benefit the
learning process for the other.

Conversely, the distinct characteristics of the femur class prevented similar
cross-class learning benefits. The output generated by our model when
processing images from various classes is depicted in Fig. 4. In cases of
success, as shown in Fig. 4(a), the model’s predictions closely align with the
ground truth, resulting in favorable Mean Average Error (MAE) values. However,
challenges are evident in Fig. 4(b), which showcases instances of the model
struggling with accurate segmentation, subsequently affecting the estimation.
Notably, for the brain class, the segmentation mask appears sound, yet the MAE
reaches around 6mm. Conversely, for the femur, complications arise due to
mixed outputs with other classes, adversely impacting the predictions.

(a) Success cases

(b) Failure cases

Figure 4: Model Predictions: Each prediction incorporates printed estimations
from the model. The prediction showcases the estimated biometric parameter
(left) alongside the corresponding ground truth value (right). Further, the
predicted class (left) and the actual class (right) are presented below these
values.

##  5 Conclusion

In this work, an integrated approach is proposed to simultaneously classify,
segment, and estimate parameters of various body organs in a fetus. The method
utilizes a multi-task learning technique. By merging the cross-entropy loss
with the dice loss, our model is designed to effectively estimate the
parameters associated with each class. Unlike previous studies, our approach
involves the calculation of multiple bio-parameters and capitalizes on
insights from other classes to enhance overall performance. Additionally, we
conduct experiments to evaluate the impact of varying loss weights within the
multi-task learning framework. To evaluate the effectiveness of our approach,
we compared it with individually trained segmentation models for each organ
class. The results demonstrated noticeable improvements in terms of MAE and
classification accuracy. The integrated model outperformed the individual
models across evaluated organ classes. One of the key advantages of our method
is its versatility, as it can be applied to multiple body organs, rather than
focusing on a single organ. This offers clinicians a comprehensive tool to
assess various bio-parameters simultaneously. Additionally, the computational
cost of our model remains manageable, as it covers multiple organ parameters
within a single framework. Future work can be done by extending this framework
to more organ estimations.

## References

  * [1] Jamie Perin, Amy Mulick, Diana Yeung, Francisco Villavicencio, Gerard Lopez, Kathleen L Strong, David Prieto-Merino, Simon Cousens, Robert E Black, and Li Liu.  Global, regional, and national causes of under-5 mortality in 2000–19: an updated systematic analysis with implications for the sustainable development goals.  The Lancet Child & Adolescent Health, 6(2):106–115, 2022.
  * [2] Leo Joskowicz, D Cohen, N Caplan, and Jacob Sosna.  Inter-observer variability of manual contour delineation of structures in ct.  European radiology, 29:1391–1399, 2019.
  * [3] Andrés Anaya-Isaza, Leonel Mera-Jiménez, and Martha Zequera-Diaz.  An overview of deep learning in medical imaging.  Informatics in medicine unlocked, 26:100723, 2021.
  * [4] Mehmet Murat Seval and Bulut Varlı.  Current developments in artificial intelligence from obstetrics and gynecology to urogynecology.  Frontiers in Medicine, 10, 2023.
  * [5] Maria Chiara Fiorentino, Francesca Pia Villani, Mariachiara Di Cosmo, Emanuele Frontoni, and Sara Moccia.  A review on deep-learning algorithms for fetal ultrasound-image analysis.  Medical Image Analysis, 83:102629, jan 2023.
  * [6] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In Nassir Navab, Joachim Hornegger, William M. Wells, and Alejandro F. Frangi, editors, Medical Image Computing and Computer-Assisted Intervention – MICCAI 2015, pages 234–241, Cham, 2015. Springer International Publishing.
  * [7] Karen Simonyan and Andrew Zisserman.  Very deep convolutional networks for large-scale image recognition.  arXiv preprint arXiv:1409.1556, 2014.
  * [8] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
  * [9] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al.  Pytorch: An imperative style, high-performance deep learning library.  Advances in neural information processing systems, 32, 2019.


[a,b]Rainer Sommer

# Log-enhanced discretization errors in integrated correlation functions

Leonardo Chimirri  Nikolai Husung

###### Abstract

Integrated time-slice correlation functions $G(t)$ with weights $K(t)$ appear,
e.g., in the moments method to determine $\alpha_{s}$ from heavy quark
correlators, in the muon g-2 determination or in the determination of smoothed
spectral functions.

For the (leading-order-)normalised moment $R_{4}$ of the pseudo-scalar
correlator we have non-perturbative results down to $a=10^{-2}$ fm and for
masses, $m$, of the order of the charm mass in the quenched approximation. A
significant bending of $R_{4}$ as a function of $a^{2}$ is observed at small
lattice spacings.

Starting from the Symanzik expansion of the integrand we derive the asymptotic
convergence of the integral at small lattice spacing in the free theory and
prove that the short distance part of the integral leads to $\log(a)$-enhanced
discretisation errors when $G(t)K(t)\raisebox{-0.43057pt}{
$\stackrel{{\scriptstyle\small{t\to 0}}}{{\sim}}$}\,t$ for small $t$. In the
interacting theory an unknown, function $K(a\Lambda)$ appears.

For the $R_{4}$-case, we modify the observable to improve the short distance
behavior and demonstrate that it results in a very smooth continuum limit. The
strong coupling and the $\Lambda$-parameter can then be extracted. In general,
and in particular for $g-2$, the short distance part of the integral should be
determined by perturbation theory. The (dominating) rest can then be obtained
by the controlled continuum limit of the lattice computation.

##  1 Introduction

We consider a $\bf p=0$ (spatial momentum zero) correlator

$G(t,M,a)=a^{3}\sum_{\bf x}\langle
P^{\mathrm{RGI}}(x)\overline{P}^{\mathrm{RGI}}(0)\rangle=G(t,M,0)+\Delta
G(t,M,a)\,,$ (1)

of a renormalization group invariant (RGI) local field $P^{\mathrm{RGI}}$ of
dimension three with non-trivial quantum numbers such that the vacuum does not
contribute as intermediate state. The RGI mass of the theory (or the set of
masses) is denoted by $M$ and $\Delta O$ denotes the lattice artefact of an
observable $O$. Weighted integrals $\int G(t,M,0)\,K(t)\,\mathrm{d}t$ , such
as moments, need a weight $K(t)\raisebox{-0.43057pt}{
$\stackrel{{\scriptstyle\small{t\to 0}}}{{\sim}}$}t^{n}\,,n>2$ to ensure
convergence at small $t$.111Fields of other dimensions or integrals of the
type $\int\langle
P^{\mathrm{RGI}}(x)\overline{P}^{\mathrm{RGI}}(0)\rangle\,\tilde{K}(x)\,\mathrm{d}^{4}x$
lead to trivial changes of our discussion.  Specializing to moments with
$n>2$, one can then also consider

$\mathcal{M}_{n}(M,a)=a\sum_{t}t^{n}\,G(t,M,a)=\mathcal{M}_{n}(M,0)+\Delta\mathcal{M}_{n}(M,a)\,,$
(2)

with a finite continuum limit $\mathcal{M}_{n}(M,0)$. The case $n=4$ will be
discussed in detail since it is of particular interest for computing
$\alpha_{s}$, when $P$ is a heavy-quark bilinear [1] and furthermore the
hadronic vacuum polarization contribution to $g-2$ of the muon has the form
above in the time-momentum representation [2] with a
$K(t)\raisebox{-0.43057pt}{ $\stackrel{{\scriptstyle\small{t\to
0}}}{{\sim}}$}t^{4}$. We will comment on other moments as we go along. In the
following we assume mass-degenerate quarks to simplify the notation.

Note that in the heavy quarks moments method for determining $\alpha_{s}$ one
typically considers the dimensionless

$\overline{\mathcal{M}}_{4}(M,a)=M^{2}\mathcal{M}_{4}(M,a)\,,$ (3)

with $M$ the RGI-mass, such that also $\overline{\mathcal{M}}_{4}$ is scale
invariant. Specifically one chooses
$P^{\mathrm{RGI}}=Z^{\mathrm{RGI}}P^{\mathrm{bare}},\;P^{\mathrm{bare}}=\bar{c}\gamma_{5}c^{\prime}$
and for discretizations with enough chiral symmetry the renormalization factor
$Z^{\mathrm{RGI}}$ is not needed due to
$MP^{\mathrm{RGI}}=m_{\mathrm{bare}}P^{\mathrm{bare}}$. The correlator $G$,
eq. 1, is even under time-refections, $G(t,M,0)=G(-t,M,0)$. Thus moments for
odd $n$ vanish and only moments with $n\geq 4$ are finite.

In an $\mathrm{O}(a)$-improved theory, the Symanzik effective theory
prediction (SymEFT) [3, 4, 5] is

$\Delta G\raisebox{-0.43057pt}{ $\stackrel{{\scriptstyle\small{a\to
0}}}{{\sim}}$}a^{2}[\alpha_{s}(1/a)]^{\hat{\gamma}_{\mathrm{lead}}}\,.$ (4)

Naively one may expect that this also leads to
$\Delta\mathcal{M}_{n}\raisebox{-0.43057pt}{
$\stackrel{{\scriptstyle\small{a\to
0}}}{{\sim}}$}a^{2}[\alpha_{s}(1/a)]^{\hat{\gamma}_{\mathrm{lead}}}\,.$ Here
we discuss that this is not the case and show that a safe continuum limit
cannot even be taken with lattice spacings down to $a=10^{-2}$fm (section 2).
We derive that already in the free theory an $a^{2}\log(aM)$ term is present
(section 3) and sketch what changes in the SymEFT prediction in the
interacting theory (section 4). Since the general conclusion is that integrals
such as the one defining $\mathcal{M}_{4}$ cannot be computed reliably on the
lattice, we then propose a modification for $\mathcal{M}_{4}$ (section 6.1)
and demonstrate that it works very well. Finally we also make a simple and
practical proposal which solves the issue for the HVP contribution to the muon
$g-2$ (section 6.3).

##  2 Demonstration of the deviations from simple $a^{2}$ scaling

We computed $\overline{\mathcal{M}}_{4}$ (and other moments [6]) in the
quenched approximation on ensembles sft7 - sft4 [7] and q_b649 - q_b616 [6]
with lattice spacings $a={0.01\,\mathrm{fm}}\times 2^{n/2},n=0\ldots 6$, i.e.
$0.01\,\mathrm{fm}\leq a\leq 0.08\,\mathrm{fm}$. The property
$MP^{\mathrm{RGI}}=m_{\mathrm{bare}}P^{\mathrm{bare}}$ is guaranteed by using
the twisted mass formulation at maximal twist and double insertions of the
Pauli term in SymEFT are avoided by including the Sheikholeslami-Wohlert term
[8] with non-perturbative improvement coefficient [9]. Further details are
given in [10].

In fig. 1 we show the lattice spacing dependence of

$R_{4}=\frac{\overline{\mathcal{M}}_{4}}{[\overline{\mathcal{M}}_{4}]_{g=0,a>0}}\,.$
(5)

The normalization by the lattice leading perturbative order (finite $a>0$) is
crucial as seen by the points with continuum norm,
$\overline{\mathcal{M}}_{4}/[\overline{\mathcal{M}}_{4}]_{g=0,a=0}$. Again we
refer to [6, 10] for more details. However, despite the strong reduction of
discretisation errors by the lattice norm, a continuum extrapolation with data
in the range $a\in[0.02,0.04]$ fm (linear fit in fig. 1) where $R_{4}$
seemingly scales with $a^{2}$ corrections, clearly leads to a wrong result.
This is seen by the $a=0.01$ fm data point and corroborated by our method
sketched in section 6.1. Such a behavior is the nightmare of numerical
analysis. Note that the mass $M\approx M_{\mathrm{charm}}$ is not that high.

Figure 1: Lattice dependence of
$R_{4}=\overline{\mathcal{M}}_{4}/[\overline{\mathcal{M}}_{4}]_{g=0}$, where
the normalization is performed with $[\overline{\mathcal{M}}_{4}]_{g=0}$ at
finite $a$ (“lattice norm”) or at $a=0$ (“continuum norm”). The quark mass is
around the charm quark mass.

##  3 Derivation of the $a^{2}\log(aM)$ term in the free theory

In this and the following section we study the small $t$ behavior where mass-
effects can be neglected and we first consider the contribution to
$\mathcal{M}_{4}$ from a range $t_{1}\leq t\leq t_{2}\ll 1/M$,

$\displaystyle\Delta I(t_{1},t_{2})$ $\displaystyle=$ $\displaystyle
2a\sum_{t=t_{1}}^{t_{2}}w_{\mathrm{T}}(t)\,t^{4}\,G(t,M,a)-I_{\mathrm{cont}}(t_{1},t_{2})\,,\quad
t_{1}M\ll 1,\;t_{2}M\ll 1\,,$ (7) $\displaystyle
I_{\mathrm{cont}}(t_{1},t_{2})=2\int_{t_{1}}^{t_{2}}\mathrm{d}t\;t^{4}\,G(t,M,0)\,,\quad
t_{1}M\ll 1,\;t_{2}M\ll 1\,.$

The weight $w_{\mathrm{T}}(t)$ implements the trapezoidal rule: it is $1/2$ at
the boundaries and $1$ otherwise.

In order to gain understanding, we start with the free theory, $g=0$. This
case is illuminating and at the same time we can get the relevant result by
dimensional reasoning alone.

We split

$\displaystyle\Delta I(0,t)$ $\displaystyle=$ $\displaystyle\Delta
I(0,t_{1})+\Delta I(t_{1},t)\,,$ (8)

discuss the second term and then add the first one. The SymEFT prediction for
the cutoff effects of $G$ are

$\displaystyle\Delta G$ $\displaystyle=$ $\displaystyle
k_{\mathrm{L}}\,\frac{a^{2}}{t^{5}}+\mathrm{O}(a^{4})+\mathrm{O}(M^{2}t^{2})\,,$
(9)

with a constant $k_{\mathrm{L}}$ which depends on the fermion
discretization.222This form is simply due to dimensional counting. $G(t)$ has
mass dimension $-3$ and therefore behaves like $\sim t^{-3}$ for small $t$ in
the free theory. In the interacting theory there are log-corrections to that
functional form due to anomalous dimensions of $P$ and the SymEFT operators.
Relative cutoff effects are $\sim a^{2}/t^{2}$, again because for $tM\ll 1$
the only dimensionful parameter apart from $a$ is $t$. Performing an explicit
leading order computation, expanded in $a/t$ in the Wilson regularization we
find $k_{\mathrm{L}}=1$. Since mass-effects are irrelevant, $k_{\mathrm{L}}=1$
holds irrespective of whether we choose a twisted mass term or a standard one.
Not indicating the higher order corrections in $a$ and $M$ any further we get

$\displaystyle\Delta I(t_{1},t)$ $\displaystyle\overset{a\ll t_{1}}{\sim}$
$\displaystyle\;k_{\mathrm{L}}\,a^{2}\
\,\int_{t_{1}}^{t}\,\mathrm{d}s\,s^{-1}+\Delta I_{\mathrm{T}}(t_{1},t)\,$ (10)
$\displaystyle=$ $\displaystyle k_{\mathrm{L}}a^{2}\,\,\log(t/t_{1})+\Delta
I_{\mathrm{T}}(t_{1},t)=k_{\mathrm{L}}a^{2}\,[\log(t/a)-\log(t_{1}/a)]+\Delta
I_{\mathrm{T}}(t_{1},t)\,.$ (11)

Here, $\Delta I_{\mathrm{T}}\sim a^{2}$ is the error in using the trapezoidal
rule for the integral. We drop it because it does not play a role in the
following; it is regular as $t\to 0$ and does not introduce a log. We then
obtain

$\displaystyle\Delta I(0,t)$ $\displaystyle=$ $\displaystyle\underbrace{\Delta
I(0,t_{1})-k_{\mathrm{L}}a^{2}\log(t_{1}/a)}_{=k\,a^{2}}+k_{\mathrm{L}}\,a^{2}\,\log(t/a)=a^{2}\,[k+k_{\mathrm{L}}\log(t/a)]\,,$
(12)

with another dimensionless constant $k$ depending on the regularization. The
first term, $ka^{2}$, has this form because it neither depends on $t_{1}$ nor
on $t$ and $a$ is the only dimensionful parameter.

For the full moment, $t$ gets replaced by the only physics scale of the
integral, namely $1/M$. We thus arrive at

$\displaystyle\frac{\Delta\mathcal{M}_{4}(M,a)}{\mathcal{M}_{4}(M,0)}$
$\displaystyle=$ $\displaystyle
a^{2}M^{2}\,[k^{\prime}-k_{\mathrm{L}}\log(M\,a)]+\mathrm{O}(M^{4}a^{4})\,\,.$
(13)

We note that [11] have argued for the presence of a $\log(t/a)$ term in the
same discretised integral (in the context of HVP). In contrast to their
argumentation, we never work with divergent integrals or with the Symanzik
expansion for $a/t=\mathrm{O}(1)$.

It is instructive to add higher order terms,
$k_{d}\,(a^{2}/t^{5})\,(a/t)^{d-2}$ with $d>2$ terms in the SymEFT for $\Delta
G$. They yield

$\displaystyle\Delta I(t_{1},t)$ $\displaystyle\overset{a\ll t_{1}}{\sim}$
$\displaystyle\;k_{\mathrm{L}}a^{2}\,[\log(t/a)-\log(t_{1}/a)]+a^{2}\sum_{d>2}\frac{k_{d}}{d-2}\,[\,(a/t_{1})^{d-2}-(a/t)^{d-2}\,]\,.$
(14)

and

$\displaystyle\Delta I(0,t)$ $\displaystyle\sim$
$\displaystyle\;a^{2}k^{\prime\prime}\,+\,k_{\mathrm{L}}a^{2}\,\log(t/a)-a^{2}\sum_{d>2}\frac{k_{d}}{d-2}(a/t)^{d-2}\,,$
(15)

where $k^{\prime\prime}$ now receives contributions also from the $d>2$ terms
in $\Delta G$. Note that the reasoning for the term $a^{2}k^{\prime}$ is
unchanged. It is simply the dimension of $\Delta I$ inherited from the one of
$I$ and the independence on $t_{1}$. This means that Symanzik improvement does
not hold for the integral: we could improve $\Delta G$ such that all $\sim
a^{2}$ terms are removed, but $\Delta\mathcal{M}_{4}$ would remain of order
$a^{2}$ due to the $d>2$ terms in (14). "Only" the log-term at order $a^{2}$
disappears by improvement of the integrand.

Consider for a moment the moment

$N_{3}=a\sum_{t\geq 0}t^{3}\,G(t,M,a)\,.$ (16)

In this case, we obtain $\mathrm{O}(a)$ effects, irrespective of how the
theory was improved. It is relevant to investigate whether such terms appear
in some (sub-)integrals in representations of light-by-light scattering
evaluated on the lattice [12].

##  4 SymEFT analysis beyond the free theory

It is not difficult to follow the above steps for the interacting theory. One
has to write $\Delta G$ as in eq. 4 and also the short distance behavior
changes due to anomalous dimension effects. These modifications introduce
powers of $\alpha_{s}(1/a)$ and $\alpha_{s}(1/t)$, respectively, but are not
of prime relevance. More important is that the step analogous to eq. 12 is
modified to

$ka^{2}\to K(a\Lambda)\,a^{2}\,,$ (17)

with a function $K(a\Lambda)$ which is not restricted by simple arguments.
Without knowing the behavior of $K$ at the origin, nothing can be concluded
about $M$-independent $a$-effects of the integral. The structure of external
scale dependent cutoff effects will be discussed in a publication [13]. The
basic reason for the difficulty is of course that the interacting theory has a
dynamical scale, $\Lambda$, which makes the dimensional analysis much less
restrictive.

##  5 Higher moments $\mathcal{M}_{n},\;n>4$

With $n>4$, the $a^{2}\log(aM)$ term is absent in the free theory. Still,
$\log(aM)$ dependences are present, but they are pushed to a higher order in
$a$,

$\mathcal{M}_{n}=\ldots+\mathrm{const.}\,\times
a^{n-2}\log(aM)+\ldots\,.\qquad$ (18)

##  6 Solutions

Our discussion shows that integrals of the considered type cannot be computed
on the lattice in the straight-forward way. The best solution to this problem
is to avoid integrands which have a behavior $\sim t^{k}\,,k<2$. First we
describe a specific solution for $\overline{\mathcal{M}}_{4}$ for which we
have a complete numerical demonstration. Then we propose a general solution,
which in particular will be useful for HVP.

###  6.1 A practical solution for $\overline{\mathcal{M}}_{4}$

Our simple solution for the moment $\overline{\mathcal{M}}_{4}$ uses two
different masses in the form (dropping the $a$-dependence)

$\displaystyle\rho(M_{1},M_{2})$
$\displaystyle=\frac{2\pi^{2}}{3}\,(1-r^{2})^{-1}\,[\overline{\mathcal{M}}_{4}(M_{1})-r^{2}\overline{\mathcal{M}}_{4}(M_{2})]$
(19)
$\displaystyle=\frac{2\pi^{2}}{3}\,(1-r^{2})^{-1}M_{1}^{2}[\mathcal{M}_{4}(M_{1})-\mathcal{M}_{4}(M_{2})]\,,\quad
r=M_{1}/M_{2}>1.$ (20)

The second line shows that the small $t$ asymptotics of the integrand is
improved via,

$\displaystyle t^{4}\,[G(t,M_{1})-G(t,M_{2})]\sim
t^{4}\,(t^{2}M_{2}^{2}-t^{2}M_{1}^{2})+\mathrm{O}(t^{8})\,.$ (21)

There are log-corrections to this equation in the interacting theory, which
are not relevant here. Due to the extra two powers of $t$, which come with the
mass-effects, the quantity $\rho(M_{1},M_{2})$ has no log-enhanced $a^{2}$
effects (they will appear only at the level $a^{4}$).

For the purpose of extracting $\alpha_{s}$ it is now relevant to choose
$M_{1}$ and $M_{2}$ not too different. Then the perturbative expansion, which
is given in terms of the one of $\overline{\mathcal{M}}_{4}$, does not contain
large logs of $M_{2}/M_{1}$. We write the perturbative expansion in terms
333We implicitly define
$m_{\star}=\overline{m}_{\rm\overline{MS}}(m_{\star})$. In practice, to
evaluate $\overline{\mathcal{M}}_{4}(M_{2})$ we choose
$\alpha_{\rm\overline{MS}}(m_{2\star})$ as expansion variable and use the
5-loop running of the coupling and quark mass to relate $m_{2\star}$ to
$M_{2}$. One could also obtain expansion coefficients which depend on $r$.  of
$\alpha_{s}(m_{2\star})$ with $M_{1}>M_{2}$. The
$\frac{2\pi^{2}}{3}\,(1-r^{2})^{-1}$ normalization in eq. 19 ensures

$\rho(M_{1},M_{2})=1+c_{1}\alpha_{\rm\overline{MS}}(m_{2\star})+\ldots\,,$
(22)

where $c_{1}=0.74272...$ is the same expansion coefficient as the one of
$R_{4}=\frac{2\pi^{2}}{3}\overline{\mathcal{M}}_{4}$ and higher order ones are
easily obtained. We expand in $\alpha_{\rm\overline{MS}}(m_{2\star})$ because
the difference is dominated somewhat more by long distances and $M_{2}$ is the
smaller of the masses.

Figure 2: Continuum limit extrapolations of $\rho$ and its TL improved
version, $\rho(M_{1},M_{2})^{\mathrm{Latnorm}}$. Masses, specified in units
$z_{i}=M_{i}\sqrt{8t_{0}}$, are $z_{1}=4.5,\;z_{2}=3$ (left) and
$z_{1}=13.5,\;z_{2}=9$ (right).

We show continuum limit extrapolations in fig. 2. They are almost straight in
$a^{2}$ at small $a$ which makes them quite easy to do. They can be further
improved by dividing $\rho$ by the same function evaluated at leading order,
i.e. $g=0$. There is a choice which masses to insert into the leading order
formula. A good choice is again $m_{\star}$. Precisely we define

$R_{4}^{\mathrm{TL}}(a\mu)=R_{4}|_{g=0}$ (23)

with $\mu$ the twisted mass and then

$\displaystyle\rho^{\mathrm{Latnorm}}(M_{1},M_{2})$ $\displaystyle=$
$\displaystyle\frac{3}{2\pi^{2}}\,(1-r_{\star}^{2})\,\frac{\rho(M_{1},M_{2})}{R_{4}^{\mathrm{TL}}(am_{\star
1})-r_{\star}^{2}R_{4}^{\mathrm{TL}}(am_{\star 2})}$ (24)

with

$r_{\star}=\frac{m_{\star 1}}{m_{\star 2}}\,.$ (25)

In principle it is important that $r_{\star}$ is given by the ratio of the
masses that appear in $R_{4}^{\mathrm{TL}}$ for the log-term to cancel. But
numerically, replacing $r_{\star}\to r$ makes only a small difference.
Examples for how the discretization errors are reduced can be seen in fig. 2.
For all our values of $M_{1},M_{2}$, the leading order improved
$\rho(M_{1},M_{2})^{\mathrm{Latnorm}}$ has a rather convincing continuum
extrapolation.

After the continuum extrapolation, one straight-forwardly extracts the
effective $\Lambda$-parameter and arrives at the red circles in fig. 3. These
values are computed from three-loop perturbation theory (i.e. including
$\alpha^{3}$ in $R_{4}$) at finite $\alpha(m_{\star})$. They then have a
residual dependence

$\Lambda_{\mathrm{eff}}=\Lambda+\mathrm{O}(\alpha^{2}(m_{\star}))\,,$ (26)

on $m_{\star}$ and we call them “effective”. The comparison to the Dalla Brida
and Ramos value [14], extracted at $\alpha^{2}<0.01$ with the help of a finite
size step scaling method, shows that $\Lambda$ computed from $\rho$ has at
most small (on the scale of our uncertainties) corrections at the largest
mass. That mass is given by $z=9$ or $m_{\star}\approx 2.7\,\mathrm{GeV}$.

Figure 3: $\Lambda_{\rm\overline{MS}}$ computed from
$\alpha_{\rm\overline{MS}}(m_{2\star})$, where the latter is obtained from the
non-perturbative $\rho$. The dotted line is a fit to all points including the
Dalla Brida / Ramos one [14]. The reconstructed data points are described in
the text.

###  6.2 Reconstruction of
$R_{4}=\frac{2\pi^{2}}{3}\overline{\mathcal{M}}_{4}(M)$ from $\rho$.

From the definition eq. 19 of $\rho$ it is clear that given
$\rho(M_{1},M_{2})$ and $R_{4}(M_{2})$ one can determine $R_{4}(M_{1})$. This
can be exploited by using $\rho$ to go from
$R_{4}(M_{\mathrm{ref}}\gg\Lambda)$, where perturbative uncertainties are
suppressed the most, to smaller masses.444In the opposite direction all
uncertainties in $\rho$ get enhanced, quickly leading to uncontrolled results.
We insert the known [14] $\Lambda$-parameter into the three-loop (i.e.
including $\alpha^{3}$) perturbative expression for $R_{4}$ at our highest
mass, $z_{\mathrm{ref}}=13.5$ and obtain

$\displaystyle R_{4}^{\mathrm{reconstructed}}(M)$ $\displaystyle=$
$\displaystyle(1-r^{-2})\,\rho(M,M_{\mathrm{ref}})+r^{-2}\,R_{4}^{\mathrm{3-loop}}(M_{\mathrm{ref}})\,,\;$
(28) $\displaystyle
r=M/M_{\mathrm{ref}},\;z_{\mathrm{ref}}=\sqrt{8t_{0}}M_{\mathrm{ref}}=13.5\,.$

Note that perturbative errors are small in $R_{4}(M_{\mathrm{ref}})$ as seen
in the analysis of $\rho$. They get further suppressed by a factor
$r^{-2}\approx 1/20$ when we go to $z=3$. This means that we obtain the non-
perturbative dependence of $\Lambda_{\mathrm{eff}}$ (as of now computed from
$R_{4}$ and therefore with somewhat different $\mathrm{O}(\alpha^{2})$ terms)
on $\alpha$. We remind the reader that a direct computation of $R_{4}$ was
impossible due to the $a^{2}\log(aM)$ effects.

###  6.3 Proposal for the HVP contribution to the muon $g-2$

The discussion in the previous section is easily transferred to the case of
the muon $g-2$, working with differences of the HVP integral for different
(artificial) muon masses. Additionally, we would like to advocate a very
simple solution for this and similar cases, where the short distance
contribution to the integral is subdominant. In contrast to the
$\mathcal{M}_{4}$-case the goal is not to determine $\alpha_{s}$ or other
short-distance parameters.

It is then advisable to split the integral into a short-distance part
evaluated by continuum perturbation theory and a long-distance one to be
computed on the lattice:

$\int_{0}^{\infty}\mathrm{d}tF(t)=\underbrace{\int_{0}^{\infty}\mathrm{d}t\,[1-\chi(t)]\,F(t)}_{\text{continuum
PT}}\,+\,\underbrace{a\sum_{t=0}^{\infty}\,\chi(t)\,\,F(t)}_{\text{continuum
limit of lattice
results}}\,,\quad\chi(t)\,\sim\,\begin{cases}\mathrm{O}(t^{2})&t\Lambda_{\mathrm{\overline{MS}}}\ll
1\\\ 1&t\Lambda_{\mathrm{\overline{MS}}}\gg 1\end{cases}\,.$ (29)

For example the function $\chi$ can be taken as

$\chi(t)=\frac{(M_{\mathrm{cut}}t)^{k}}{(M_{\mathrm{cut}}t)^{k}+1}\,,\;M_{\mathrm{cut}}\gg\Lambda_{\mathrm{\overline{MS}}}$
(30)

or also as a step-function, $\chi(t)=\theta(tM_{\mathrm{cut}}-1)$. The smooth
version seems advantageous for perturbation theory as well as for the lattice
discretization of the integral. The use of perturbation theory for the small
$t$-part of the integral has already been anticipated in [2]. Our discussion
adds further motivation and understanding. It suggests a smooth function
$\chi$ such as eq. 30.
Acknowledgements. We thank the organizers for a very pleasant and successful
conference. Discussions with H. Meyer, S. Kuberski and A. Risch on HVP are
gratefully acknowledged. We also thank C. Lehner for an email exchange on that
subject. LC and RS acknowledge funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Skłodowska-Curie grant
agreement No. 813942.

## References

  * [1] HPQCD collaboration, _High-Precision Charm-Quark Mass and QCD coupling from Current-Current Correlators in Lattice and Continuum QCD_ , _Phys. Rev. D_ 78 (2008) 054513 [0805.2999].
  * [2] D. Bernecker and H.B. Meyer, _Vector Correlators in Lattice QCD: Methods and applications_ , _Eur. Phys. J. A_ 47 (2011) 148 [1107.4388].
  * [3] K. Symanzik, _Continuum limit and improved action in lattice theories:(i). principles and $\varphi$4 theory_, _Nuclear Physics B_ 226 (1983) 187.
  * [4] N. Husung, P. Marquard and R. Sommer, _Asymptotic behavior of cutoff effects in Yang–Mills theory and in Wilson’s lattice QCD_ , _Eur. Phys. J. C_ 80 (2020) 200 [1912.08498].
  * [5] N. Husung, _Logarithmic corrections to O( $a$) and O($a^{2}$) effects in lattice QCD with Wilson or Ginsparg-Wilson quarks_, 2206.03536.
  * [6] L. Chimirri and R. Sommer, _Investigation of the Perturbative Expansion of Moments of Heavy Quark Correlators for $N_{f}=0$_, _PoS_ LATTICE2021 (2022) 354 [2203.07936].
  * [7] N. Husung, M. Koren, P. Krah and R. Sommer, _SU(3) Yang Mills theory at small distances and fine lattices_ , _EPJ Web Conf._ 175 (2018) 14024 [1711.01860].
  * [8] B. Sheikholeslami and R. Wohlert, _Improved Continuum Limit Lattice Action for QCD with Wilson Fermions_ , _Nucl. Phys. B_ 259 (1985) 572.
  * [9] M. Lüscher, S. Sint, R. Sommer, P. Weisz and U. Wolff, _Nonperturbative O(a) improvement of lattice QCD_ , _Nucl. Phys. B_ 491 (1997) 323 [hep-lat/9609035].
  * [10] L. Chimirri and R. Sommer, _A quenched exploration of heavy quark moments and their perturbative expansion_ , 2023.
  * [11] T. Harris, M. Cè, H.B. Meyer, A. Toniato and C. Török, _Vacuum correlators at short distances from lattice QCD_ , _PoS_ LATTICE2021 (2021) 572 [2111.07948].
  * [12] E.-H. Chao, R.J. Hudspith, A. Gérardin, J.R. Green, H.B. Meyer and K. Ottnad, _Hadronic light-by-light contribution to $(g-2)_{\mu}$ from lattice QCD: a complete calculation_, _Eur. Phys. J. C_ 81 (2021) 651 [2104.02632].
  * [13] L. Chimirri, N. Husung and R. Sommer, _in preparation_ , 2023.
  * [14] M. Dalla Brida and A. Ramos, _The gradient flow coupling at high-energy and the scale of SU(3) Yang–Mills theory_ , _Eur. Phys. J. C_ 79 (2019) 720 [1905.05147].


# The Resolution of Ambiguities in
Light Perturbation by Gravitational Waves

Thomas B. Mieling orcid:0000-0002-6905-0183, email:
<EMAIL_ADDRESS>University of Vienna, Faculty of Physics, Vienna,
Austria  TURIS Research Platform, Vienna, Austria  Stefan Palenta
orcid:0000-0002-6541-9537 University of Vienna, Faculty of Physics, Vienna,
Austria

###### Abstract

Some previously published expressions for the perturbation of light by
gravitational waves exhibit pathological behaviour in the limit of parallel
propagation. We show that this is caused by similarly pathological initial or
boundary data and can thus be remedied by implementing better-behaved initial
conditions.

##  1 Introduction

The perturbation of light by weak gravitational waves (GW’s) has been
discussed extensively in the literature [1, 2, 3, 4, 5, 6]. However, some
previously published results yield ill-behaved expressions in the parallel
limit, where the electromagnetic and gravitational waves propagate in the same
direction. Specifically, this problem arises due to the occurrence of
expressions of the form

$\varPsi:=\frac{1}{2}\frac{\t A{{}^{i}{}^{j}}\t m{{}_{i}}\t m{{}_{j}}}{1-\t
n{{}^{k}}\t m{{}_{k}}}\,,$ (1.1)

where $m$ and $n$ are unit three-vectors (in the flat background geometry)
along the light ray and the GW, respectively, and $\t A{{}^{i}{}^{j}}$ is the
GW polarisation tensor (which is transverse to $n$ and traceless). Indeed,
setting

$\displaystyle(\t m{{}_{i}})$ $\displaystyle=\begin{pmatrix}\cos\theta\\\
\sin\theta\cos\phi\\\ \sin\theta\sin\phi\end{pmatrix}\,,$ $\displaystyle(\t
n{{}_{i}})$ $\displaystyle=\begin{pmatrix}1\\\ 0\\\ 0\end{pmatrix}\,,$
$\displaystyle(\t A{{}^{i}{}^{j}})$ $\displaystyle=\begin{pmatrix}0&0&0\\\
0&A_{+}&A_{\times}\\\ 0&A_{\times}&-A_{+}\end{pmatrix}\,,$ (1.2)

one obtains

$\lim_{\theta\to 0}\varPsi=\lim_{\theta\to
0}\cos^{2}(\theta/2)[A_{+}\cos(2\phi)+A_{\times}\sin(2\phi)]=A_{+}\cos(2\phi)+A_{\times}\sin(2\phi)\,.$
(1.3)

But in the limit $\theta\to 0$, the angle $\phi$ becomes meaningless so that
this limit depends on the direction from which $m$ tends towards $n$. Hence,
$\varPsi$ _cannot_ be extended continuously to $m=n$.

Such ill-behaved expressions were obtained in various perturbative
calculations concerned with null geodesics [Eq. (3.15)]Finn2009[Sect.
III]Angelil2015, the electromagnetic field tensor [Eq. (2.20)]Lobo1992[Eq.
(3.5)]Montanari1998[Sect. 3]Calura1999, and the electromagnetic gauge
potential [Eq. (2.19)]Calura1999[Eq. (3.11)]Codina1980[Eq. (26)]Park2021. (In
Ref. [11] the singularity was found to be removable by a similarly singular
gauge transformation.) Here, we show that such pathological behaviour can be
traced back to ill-behaved initial data of the geodesic equation, or boundary
data for the eikonal equation, and can be remedied by implementing well-
behaved data instead.

Specifically, we consider the perturbation of null geodesics in Section 2 and
perturbations of the eikonal function in Section 3. In both cases, we show
explicitly how these ill-behaved expressions arise and give corrected
solutions (constructed from well-behaved data) which are free of this
pathological behaviour.

##  2 Null Geodesics

Within the framework of linearised gravity where the metric $\t
g{{}_{\mu}{}_{\nu}}$ is close to the Minkowski metric
$(\t\eta{{}_{\mu}{}_{\nu}})=\operatorname{diag}(-1,1,1,1)$, plane
gravitational waves can be described in the standard transverse-traceless (TT)
gauge as

$\t g{{}_{\mu}{}_{\nu}}=\t\eta{{}_{\mu}{}_{\nu}}+\varepsilon\t
h{{}_{\mu}{}_{\nu}}(\t\kappa{{}_{\sigma}}\t
x{{}^{\sigma}})+O(\varepsilon^{2})\,,$ (2.1)

where $\varepsilon\ll 1$ is the GW amplitude, $\t\kappa{{}_{\mu}}$ its wave
vector, and $\t h{{}_{\mu}{}_{\nu}}$ satisfies $\t h{{}_{\mu}{}_{0}}=0$, $\t
h{{}_{\mu}{}_{\nu}}\t\eta{{}^{\nu}{}^{\rho}}\t\kappa{{}_{\rho}}=0$ as well as
$\t h{{}_{\mu}{}_{\nu}}\t\eta{{}^{\mu}{}^{\nu}}=0$. Denoting by $\omega_{g}$
the GW frequency and setting

$\displaystyle\t\kappa{{}_{\mu}}$ $\displaystyle=\omega_{g}(-1,+\t
n{{}_{i}})\,,$ $\displaystyle\t{\bar{\kappa}}{{}_{\mu}}$
$\displaystyle=\omega_{g}(-1,-\t n{{}_{i}})\,,$ (2.2)

whose indices are raised and lowered with $\t\eta{{}_{\mu}{}_{\nu}}$ ($\t
n{{}_{i}}$ is the direction of GW propagation), we define the null coordinates
$U$, $V$, and transverse spatial coordinates $(y^{a})=(y^{1},y^{2})$:

$\displaystyle\omega_{g}U$ $\displaystyle=-\t\kappa{{}_{\mu}}\t
x{{}^{\mu}}\,,$ $\displaystyle\omega_{g}V$
$\displaystyle=-\t{\bar{\kappa}}{{}_{\mu}}\t x{{}^{\mu}}\,,$
$\displaystyle(y^{a})$ $\displaystyle=\t*b{{}^{(a)}_{i}}\t x{{}^{i}}\,,$ (2.3)

where $(n,\t b{{}^{(1)}},\t b{{}^{(2)}})$ is any (fixed) orthonormal basis of
Euclidean space (the time components of $n$ and $\t b{{}^{(a)}}$ are set to
zero). The metric then takes the form

$g=-\mathrm{d}U\,\mathrm{d}V+(\t\delta{{}_{a}{}_{b}}+\varepsilon\t
h{{}_{a}{}_{b}}(-\omega_{g}U))\,\mathrm{d}\t y{{}^{a}}\,\mathrm{d}\t
y{{}^{b}}+O(\varepsilon^{2})\,,$ (2.4)

which can be regarded as an approximate version of a Rosen metric

$g=-\mathrm{d}U\,\mathrm{d}V+\t g{{}_{a}{}_{b}}(U)\,\mathrm{d}\t
y{{}^{a}}\,\mathrm{d}\t y{{}^{b}}\,,$ (2.5)

which also allows for solutions to the full nonlinear Einstein vacuum
equations [§ 109]Landau1994[Sect. 24.5]Stephani2003Cropp2010. Since the
computation of light rays in the metric (2.5) is no more difficult than in
(2.4), we consider the exact case first and specialise to the perturbative
setting afterwards.

###  2.1 Non-Perturbative Solution

Using the obvious Killing vector fields of (2.5), one obtains the constant
momenta

$\displaystyle\t p{{}_{V}}$ $\displaystyle=-\tfrac{1}{2}\dot{U}\,,$
$\displaystyle\t p{{}_{a}}$ $\displaystyle=\t
g{{}_{a}{}_{b}}(U)\t{\dot{y}}{{}^{a}}\,.$ (2.6)

Integrating the equations $\t p{{}_{V}}=\text{const.}$ and $\t
p{{}_{a}}=\text{const.}$ determines $U$ and $\t y{{}^{a}}$ as a function of
the affine parameter $s$. Finally, $V(s)$ is obtained from the equation $4\t
p{{}_{U}}\t p{{}_{V}}=\t g{{}^{a}{}^{b}}\t p{{}_{a}}\t p{{}_{b}}$, leading to
the result

$\displaystyle U(s)$ $\displaystyle=U(0)-2\t p{{}_{V}}s\,,$ (2.7)
$\displaystyle V(s)$ $\displaystyle=V(0)+\frac{1}{2\t p{{}_{V}}}\int_{0}^{s}\t
g{{}^{a}{}^{b}}(U(s^{\prime}))\mathrm{d}s^{\prime}\,\t p{{}_{a}}\t
p{{}_{b}}\,\,,$ (2.8) $\displaystyle\t y{{}^{a}}(s)$ $\displaystyle=\t
y{{}^{a}}(0)+\int_{0}^{s}\t
g{{}^{a}{}^{b}}(U(s^{\prime}))\,\mathrm{d}s^{\prime}\,\t p{{}_{b}}\,,$ (2.9)

cf. Ref. [17]. (We assume $\t p{{}_{V}}\neq 0$ and consider the limiting case
afterwards.) Defining

$\t G{{}^{\mu}{}^{\nu}}(s)=\int_{0}^{s}\t
g{{}^{\mu}{}^{\nu}}(U(s))\,\mathrm{d}s^{\prime}\,,$ (2.10)

the solution can be written concisely as

$\t x{{}^{\mu}}(s)=\t x{{}^{\mu}}(0)+\left(\t
G{{}^{\mu}{}^{\nu}}(s)\t\kappa{{}^{\sigma}}-\tfrac{1}{2}\t\kappa{{}^{\mu}}\t
G{{}^{\nu}{}^{\sigma}}(s)\right)\frac{\t p{{}_{\nu}}(0)\t
p{{}_{\sigma}}(0)}{\t\kappa{{}^{\lambda}}\t p{{}_{\lambda}}(0)}\,,$ (2.11)

where $\t p{{}_{\mu}}(0)$ is the initial momentum. This equation is readily
verified by noting that eq. 2.7 is obtained by contracting with
$\t\kappa{{}_{\mu}}$, eq. 2.8 by contraction with $\t{\bar{\kappa}}{{}_{\mu}}$
and eq. 2.9 by contraction with $\t b{{}^{(a)}_{\mu}}$. Note that this
equation is covariant under linear coordinate transformations.

###  2.2 Perturbative Solution

Consider now the perturbative setting in TT coordinates:

$\t g{{}^{\mu}{}^{\nu}}=\t\eta{{}^{\mu}{}^{\nu}}-\varepsilon\t
h{{}^{\mu}{}^{\nu}}(U)+O(\varepsilon^{2})\,.$ (2.12)

Since the coordinate transformation given in (2.3) is linear, eq. 2.11 also
applies to the TT coordinates. Inserting (2.12) into (2.11) and expanding the
initial momentum as $\t
p{{}_{\mu}}(0)=\t*p{{}^{{(0)}}_{\mu}}+\varepsilon\t*p{{}^{{(1)}}_{\mu}}+O(\varepsilon^{2})$,
one obtains the perturbed light rays

$\displaystyle\t x{{}^{\mu}}(s)$ $\displaystyle=\t
x{{}^{\mu}}(0)+s\,\t\eta{{}^{\mu}{}^{\nu}}\t*p{{}^{{(0)}}_{\nu}}+\varepsilon\t
x{{}^{{(1)}}{}^{\mu}}(s)+O(\varepsilon^{2})\,,$ (2.13) where $\displaystyle\t
x{{}^{{(1)}}{}^{\mu}}(s)$
$\displaystyle=s\,\t\eta{{}^{\mu}{}^{\nu}}\t*p{{}^{{(1)}}_{\nu}}-\bigg{(}\tfrac{1}{2}s\t\kappa{{}^{\mu}}\t
h{{}^{\nu}{}^{\sigma}}(u_{0})+\t
H{{}^{\mu}{}^{\nu}}(s)\t\kappa{{}^{\sigma}}-\tfrac{1}{2}\t\kappa{{}^{\mu}}\t
H{{}^{\nu}{}^{\sigma}}(s)\bigg{)}\frac{\t*p{{}^{{(0)}}_{\nu}}\t*p{{}^{{(0)}}_{\sigma}}}{\t\kappa{{}^{\lambda}}\t*p{{}^{{(0)}}_{\lambda}}}\,,$
(2.14) $\displaystyle\t H{{}^{\mu}{}^{\nu}}(s)$ $\displaystyle=\int_{0}^{s}\t
h{{}^{\mu}{}^{\nu}}(U^{(0)}(s))\,\mathrm{d}s^{\prime}\,,$ (2.15)

where $u_{0}=-\t\kappa{{}_{\rho}}\t x{{}^{\rho}}(0)/\omega$. Here, we have
used $(\t\eta{{}^{\mu}{}^{\nu}}-\varepsilon\t h{{}^{\mu}{}^{\nu}}(u_{0}))\t
p{{}_{\mu}}(0)\t p{{}_{\nu}}(0)=O(\varepsilon^{2})$.

While the last two terms in eq. 2.14 are uniformly bounded in $s$ in typical
applications (for monochromatic waves, $\t H{{}^{\mu}{}^{\nu}}$ is oscillatory
with zero average), the first two terms grow linearly with the affine
parameter $s$, a behaviour which was interpreted as an indication of the
breakdown of the perturbative expansion in Ref. [7]. These terms can be
eliminated by an appropriate choice of initial momentum:

$\t*p{{}^{*}_{\mu}}=\t*p{{}^{{(0)}}_{\mu}}+\tfrac{1}{2}\varepsilon\t\kappa{{}_{\mu}}\t
h{{}^{\nu}{}^{\sigma}}(u_{0})\frac{\t*p{{}^{{(0)}}_{\nu}}\t*p{{}^{{(0)}}_{\sigma}}}{\t\kappa{{}^{\lambda}}\t*p{{}^{{(0)}}_{\lambda}}}+O(\varepsilon^{2})\,,$
(2.16)

which leads to the following _tentative_ first-order correction of the
trajectory:

$\t*x{{}_{*}^{{(1)}}{}^{\mu}}=\bigg{(}\tfrac{1}{2}\t\kappa{{}^{\mu}}\t
H{{}^{\nu}{}^{\sigma}}(s)-\t
H{{}^{\mu}{}^{\nu}}(s)\t\kappa{{}^{\sigma}}\bigg{)}\frac{\t*p{{}^{{(0)}}_{\nu}}\t*p{{}^{{(0)}}_{\sigma}}}{\t\kappa{{}^{\lambda}}\t*p{{}^{{(0)}}_{\lambda}}}\,,$
(2.17)

which was also obtained in Ref. [7] using renormalisation techniques. However,
this expression has a pathological parallel limit (recall that we have
excluded the case $\t p{{}_{V}}=0$). To see this, consider a polarised
monochromatic GW with wave vector $\t\kappa{{}_{\mu}}=\omega_{g}(-1,\t
n{{}_{i}})$ and waveform ${\t h{{}^{a}{}^{b}}=\t
A{{}^{a}{}^{b}}\cos(\omega_{g}U)}$. Choosing the initial momentum
$\t*p{{}^{{(0)}}_{\mu}}=(-1,\t m{{}_{i}})$, and using the parametrisation
(1.2), one finds

$\t H{{}^{a}{}^{b}}(s)=\t
A{{}^{a}{}^{b}}\frac{\sin(\omega_{g}s(1-\cos\theta))}{\omega_{g}(1-\cos\theta)}\,,$
(2.18)

which tends to $s\t A{{}^{a}{}^{b}}$ as $\theta\to 0$. In this limit one
obtains

$\lim_{\theta\to
0}\t*x{{}_{*}^{{(1)}}{}^{\mu}}=-s\,\t\kappa{{}^{\mu}}\omega_{g}^{-1}[A_{+}\cos(2\phi)+A_{\times}\sin(2\phi)]\,,$
(2.19)

which depends on the geometrically undefined angle $\phi$. This is the
pathological behaviour described in the introduction.

In fact, this can be traced back entirely to the similarly pathological
behaviour of the initial momentum $p^{*}$: inserting the above parametrisation
into eq. 2.16 one finds

$\lim_{\theta\to
0}\t*p{{}^{*}_{\mu}}=\t\kappa{{}_{\mu}}\omega_{g}^{-1}\left(1-\varepsilon[h_{+}(u_{0})\cos(2\phi)+h_{\times}(u_{0})\sin(2\phi)]\right)+O(\varepsilon^{2})\,,$
(2.20)

which also depends on the geometrically undefined angle $\phi$. It is thus
clear that the problematic behaviour of eq. 2.17 can be cured only by choosing
a different initial momentum.

The fact that the parallel limit yields the correct momentum only up to an
ill-behaved multiplicative factor suggests that the problem can be remedied by
changing the parametrisation of the curve. This is indeed the case: adding a
multiple of $\t*p{{}^{{(0)}}_{\nu}}$ to $\t*p{{}^{*}_{\nu}}$ such that the
overall momentum satisfies $\t p{{}_{0}}(0)=-1$, say, one is led to the
_improved_ initial momentum

$\t*p{{}^{\dagger}_{\mu}}=\t*p{{}^{{(0)}}_{\mu}}+\tfrac{1}{2}\varepsilon(\t\kappa{{}_{\mu}}-\omega_{g}\t*p{{}^{{(0)}}_{\mu}})\t
h{{}^{\nu}{}^{\sigma}}(0)\frac{\t*p{{}^{{(0)}}_{\nu}}\t*p{{}^{{(0)}}_{\sigma}}}{\t\kappa{{}^{\lambda}}\t*p{{}^{{(0)}}_{\lambda}}}+O(\varepsilon^{2})\,,$
(2.21)

which _is_ well-behaved in the collinear limit: $\lim_{\theta\to
0}\t*p{{}^{\dagger}_{\mu}}=\t*p{{}^{{(0)}}_{\mu}}$. With this choice, the
first-order correction of the light ray now takes the form

$\t x{{}^{{(1)}}{}^{\mu}}=\bigg{(}\tfrac{1}{2}\t\kappa{{}^{\mu}}\t
H{{}^{\nu}{}^{\sigma}}(s)-\t
H{{}^{\mu}{}^{\nu}}(s)\t\kappa{{}^{\sigma}}-\tfrac{1}{2}\omega_{g}s\,\t
p{{}^{{(0)}}{}^{\mu}}\t
h{{}^{\nu}{}^{\sigma}}(0)\bigg{)}\frac{\t*p{{}^{{(0)}}_{\nu}}\t*p{{}^{{(0)}}_{\sigma}}}{\t\kappa{{}^{\lambda}}\t*p{{}^{{(0)}}_{\lambda}}}\,.$
(2.22)

One readily verifies that this expression vanishes for $\theta\to 0$: as one
would expect, light rays are not perturbed by co-propagating gravitational
waves. Because the linearly growing term here is proportional to the
unperturbed tangent vector, it _does not_ describe a physical separation from
the unperturbed ray: the physical first-order perturbation of the curve always
remains bounded whenever $\t H{{}^{\mu}{}^{\nu}}$ does.

While the construction of $\t*p{{}^{\dagger}_{\mu}}$ given here might appear
ad-hoc, it resolves the issues explained in the introduction. It is also
characterised geometrically by the following three conditions: (i)
$p^{\dagger}$ is lightlike; (ii) $\t*p{{}^{\dagger}_{0}}=-1$ so that the
initial energy is normalised; (iii) there is a constant $\lambda$ such that
for all three Killing vectors $K=\partial_{v},\partial_{y},\partial_{z}$ one
has $\t*p{{}^{\dagger}_{\mu}}\t K{{}^{\mu}}=\lambda\t*p{{}^{{(0)}}_{\mu}}\t
K{{}^{\mu}}$.

(a) Pathological light cones computed from the curves (2.17).

(b) Regular light cones computed from the curves (2.22).

Figure 1: Spatial light cones of decreasing opening angle. The left panel (a)
demonstrates the direction-dependence of the limit $\theta\to 0$ of eq. 2.17:
coming from above, one covers a larger coordinate distance than coming from
the side. Contrarily, the cones constructed from eq. 2.22 are free of this
ambiguity (b). For visualisation, we have set $\omega_{g}=1$,
$\varepsilon=0.1$ and have drawn the light rays up to an affine parameter of
$s=2\pi$.

Figure 1 shows _spatial_ light cones of decreasing opening angle obtained by
varying $\phi$ over the range $(0,2\pi)$ and varying the affine parameter $s$
up to one GW wavelength. Figure 1(a) shows the cones obtained from eq. 2.17
which does not extend continuously to the parallel case: when the opening
angle is shrunk to zero, the covered coordinate distance depends on the
direction from which the limit is taken. Conversely, eq. 2.22 is free of this
issue, as shown in Figure 1(b).

(a) Pathological light sphere computed from the curves (2.17).

(b) Regular light sphere computed from the curves (2.22).

Figure 2: “Light Spheres” obtained by following light rays from a point along
all possible directions for an affine parameter of $s=2\pi$ (with
$\omega_{g}=1$ and $\varepsilon=0.1$ for visualisation). The discontinuity of
eq. 2.17 produces a “fold” in the direction of GW propagation (a), while the
sphere constructed from eq. 2.22 is continuous there (b).

This discontinuity can be visualised also by constructing “spheres” by varying
the emission angles $\theta,\phi$ while keeping the affine parameter, $s$,
fixed. The resulting (spatial) wave fonts are shown in Figure 2. The
discontinuity of eq. 2.17 results in a “notch”, see Figure 2(a), while eq.
2.22 gives rise to the regular surface shown in Figure 2(b).

As the two curves corresponding to the initial momenta $p^{*}$ and
$p^{\dagger}$ are equivalent in the sense that one can be obtained from the
other by re-parametrisation, their predictions for the light trajectories
coincide. However, care must be taken when the affine parameter is interpreted
as an optical phase (as is the case in geometrical optics). Indeed, similar
pathologies arise in the eikonal equation, as discussed in the next section.

##  3 Eikonal Equation

One of the fundamental quantities of geometrical optics is the eikonal
function $\psi$, satisfying the eikonal equation

$\t g{{}^{\mu}{}^{\nu}}(\t\nabla{{}_{\mu}}\psi)(\t\nabla{{}_{\nu}}\psi)=0\,.$
(3.1)

In this description, light rays are modelled by integral curves of
$\t\nabla{{}^{\mu}}\psi$, which are null geodesics along which the eikonal is
constant. To describe the perturbation of plane light waves, consider
perturbative solutions of the form

$\psi=\t*k{{}^{{(0)}}_{\mu}}\t
x{{}^{\mu}}+\varepsilon\psi^{(1)}+O(\varepsilon^{2})\,,$ (3.2)

where $\t*k{{}^{{(0)}}_{\mu}}$ is a constant null vector in the flat
background geometry. Inserting this into eq. 3.1 with the metric given in eq.
2.12 and neglecting terms of order $\varepsilon^{2}$, one obtains

$\t*k{{}^{{(0)}}{}^{\mu}}\t\partial{{}_{\mu}}\psi^{(1)}=\tfrac{1}{2}\t
h{{}^{\mu}{}^{\nu}}(\t\kappa{{}_{\sigma}}\t
x{{}^{\sigma}})\t*k{{}^{{(0)}}_{\mu}}\t*k{{}^{{(0)}}_{\mu}}\,,$ (3.3)

where
$\t*k{{}^{{(0)}}{}^{\mu}}=\t\eta{{}^{\mu}{}^{\nu}}\t*k{{}^{{(0)}}_{\nu}}$. For
a unique solution one must prescribe boundary values on a suitable
hypersurface $\Sigma$, for which we introduce some notation. Given any point
$\t x{{}^{\mu}}$, denote by $\t*x{{}_{\text{ret.}}^{\mu}}$ the point obtained
by following the unperturbed light rays back in time (along $-\t
k{{}^{{(0)}}{}^{\mu}}$) until one reaches $\Sigma$. This can be written as

$\t*x{{}_{\text{ret.}}^{\mu}}=\t x{{}^{\mu}}-s(x)\t k{{}^{{(0)}}{}^{\mu}}\,,$
(3.4)

where the “parametric distance” $s$ satisfies

$\displaystyle\t k{{}^{{(0)}}{}^{\mu}}\t\partial{{}_{\mu}}s$
$\displaystyle=1\,,$ $\displaystyle s|_{\Sigma}=0\,.$ (3.5)

The solution to eq. 3.3 with prescribed boundary data
$\psi^{(1)}|_{\Sigma}=\delta\psi$ can then be written as

$\psi^{(1)}=\frac{\t k{{}^{{(0)}}_{\mu}}\t
k{{}^{{(0)}}_{\mu}}}{2\t\kappa{{}^{\lambda}}\t
k{{}^{{(0)}}_{\lambda}}}\int_{\t\kappa{{}_{\sigma}}\t*x{{}_{\text{ret.}}^{\sigma}}}^{\t\kappa{{}_{\sigma}}\t
x{{}^{\sigma}}}\t
h{{}^{\mu}{}^{\nu}}(u)\,\mathrm{d}u+\delta\psi(x_{\text{ret.}})\,.$ (3.6)

This is the general expression for the eikonal perturbation for arbitrary
boundary values and GW waveforms. To get some insight into this equation,
consider a monochromatic polarised GW with $\t h{{}^{\mu}{}^{\nu}}=\t
A{{}^{\mu}{}^{\nu}}\cos(U)$. In this case, one obtains

$\psi^{(1)}=\frac{\t A{{}^{\mu}{}^{\nu}}\t k{{}^{{(0)}}_{\mu}}\t
k{{}^{{(0)}}_{\mu}}}{2\t\kappa{{}^{\lambda}}\t
k{{}^{{(0)}}_{\lambda}}}[\sin(\t\kappa{{}_{\sigma}}\t
x{{}^{\sigma}})-\sin(\t\kappa{{}_{\sigma}}\t
x{{}^{\sigma}}+s(x)\t\kappa{{}_{\sigma}}\t
k{{}^{{(0)}}{}^{\sigma}})]+\delta\psi(x-sk^{(0)})\,.$ (3.7)

Further, setting $\t*k{{}^{{(0)}}_{\mu}}=\omega(-1,\t m{{}_{i}})$ for some
constant vector $\t m{{}_{i}}$ (normalised in the flat background geometry),
and choosing $\Sigma=\\{\t m{{}_{i}}\t x{{}^{i}}=0\\}$ (on which the phase is
spatially constant in the unperturbed case), one has $s=\t m{{}_{i}}\t
x{{}^{i}}/\omega$, and the eikonal perturbation takes the form

$\psi^{(1)}=-\frac{\omega\t A{{}^{i}{}^{j}}\t m{{}_{i}}\t
m{{}_{j}}}{2\omega_{g}(1-\t n{{}^{k}}\t
m{{}_{k}})}[\sin(\t\kappa{{}_{\sigma}}\t
x{{}^{\sigma}})-\sin(\t\kappa{{}_{\sigma}}\t x{{}^{\sigma}}-\omega_{g}\t
m{{}_{l}}\t x{{}^{l}}(1-\t n{{}^{k}}\t m{{}_{k}}))]+\delta\psi(x-sk^{(0)})\,.$
(3.8)

The solutions given in Refs. [Eq. (40)]Angelil2015[Eq. (26)]Park2021 and [9,
Eq. (2.19)] (using a further unidimensional approximation) correspond to the
following choice of boundary values:

$\delta\psi_{*}=-\frac{\omega\t A{{}^{i}{}^{j}}\t m{{}_{i}}\t
m{{}_{j}}}{2\omega_{g}(1-\t n{{}^{k}}\t
m{{}_{k}})}\sin(\t\kappa{{}_{\sigma}}\t x{{}^{\sigma}})\bigg{|}_{\Sigma}\,,$
(3.9)

so that the last two terms in eq. 3.8 cancel. This leads to the expression

$\displaystyle\psi^{(1)}_{*}=-\frac{\omega\t A{{}^{i}{}^{j}}\t m{{}_{i}}\t
m{{}_{j}}}{2\omega_{g}(1-\t n{{}^{k}}\t
m{{}_{k}})}\sin(\t\kappa{{}_{\sigma}}\t x{{}^{\sigma}})\,,$ (3.10)

which has precisely the pathological behaviour described in the introduction.
As for single light rays, this is caused by the pathological behaviour of the
boundary data (3.9): it has the same form as (3.10) and thus suffers from the
same problems.

This issue is resolved by choosing well-behaved boundary data instead. As
explained in detail in Refs. [18, 19], normal emission of light from $\Sigma$
with constant frequency $\omega$ corresponds to $\delta\psi\equiv 0$, which
yields instead

$\psi^{(1)}=-\frac{\omega\t A{{}^{i}{}^{j}}\t m{{}_{i}}\t
m{{}_{j}}}{2\omega_{g}(1-\t n{{}^{k}}\t
m{{}_{k}})}[\sin(\t\kappa{{}_{\sigma}}\t
x{{}^{\sigma}})-\sin(\t\kappa{{}_{\sigma}}\t x{{}^{\sigma}}-\omega_{g}\t
m{{}_{l}}\t x{{}^{l}}(1-\t n{{}^{k}}\t m{{}_{k}}))]\,.$ (3.11)

Inserting the parametrisation (1.2), it is readily verified that the parallel
limit of this expression is well-behaved: $\lim_{\theta\to 0}\psi^{(1)}=0$,
which is in agreement with the expectation light rays are not perturbed by co-
propagating gravitational waves.

###  3.1 Inadequacy of Separation of Variables

Given the simple solution of the geodesic equation in Section 2.1, one might
hope for a similarly simple solution of the eikonal equation

$-4(\t\partial{{}_{u}}\psi)(\t\partial{{}_{v}}\psi)+\t
g{{}^{a}{}^{b}}(u)(\t\partial{{}_{a}}\psi)(\t\partial{{}_{b}}\psi)=0\,,$
(3.12)

without approximations. Indeed, the standard method of separation of variables
(as commonly used in the Hamilton-Jacobi formulation of mechanics) directly
leads to the complete integral

$\psi=\t p{{}_{V}}V+\t p{{}_{a}}\t y{{}^{a}}+\frac{\t p{{}_{a}}\t
p{{}_{b}}}{4\t p{{}_{V}}}\int\t g{{}^{a}{}^{b}}(u)\,\mathrm{d}u\,,$ (3.13)

with constants $\t p{{}_{V}}$ and $\t p{{}_{a}}$, as also found in Ref. [17,
Eq. (7)]. Expanding the metric tensor $\t g{{}^{a}{}^{b}}$ to first order in
$\varepsilon$, one obtains

$\displaystyle S$ $\displaystyle=S^{(0)}+\varepsilon
S^{(1)}+O(\varepsilon^{2})\,,$ (3.14) where $\displaystyle S^{(0)}$
$\displaystyle=\t p{{}_{V}}V+\t p{{}_{a}}\t
y{{}^{a}}+\tfrac{1}{4}u\t{\delta}{{}^{a}{}^{b}}(u)\t p{{}_{a}}\t p{{}_{b}}/\t
p{{}_{V}}\,,$ (3.15) $\displaystyle S^{(1)}$ $\displaystyle=-\frac{\t
p{{}_{a}}\t p{{}_{b}}}{4\t p{{}_{V}}}\int\t
h{{}^{a}{}^{b}}(u)\,\mathrm{d}u\,.$ (3.16)

Evidently, $S^{(0)}$ satisfies the Eikonal equation in flat space and has a
well-behaved parallel limit. However, the parallel limit of $S^{(1)}$ is ill-
behaved, as it is of the same form as discussed in the previous section. Thus,
the method of separation of variables is inadequate to compute the eikonal for
arbitrary GW incidence angles.

##  4 Conclusion

We have noticed that some previously published expressions for the
perturbations of light rays and phase functions, due to gravitational waves,
exhibit pathological behaviour in the limit where light is emitted in parallel
to the gravitational wave. Contrary to the interpretation that this issue
signals a break-down of perturbative expansions [7] or even a limit in the
range of validity of geometrical optics [13], our analysis shows that this is
due the use of ad-hoc particular solutions of the equations without checking
their physical plausibility.

While we have restricted the discussion to single light rays and the eikonal
equation for simplicity of exposition, the argument carries over to more
accurate models based on wave equations (specifically those implied by
Maxwell’s equations). As shown in Ref. [19], a careful choice of boundary data
for Maxwell’s equations is necessary to obtain similarly well-behaved
solutions to the full electromagnetic field emitted from a radiating surface.

## Acknowledgements

We thank Piotr Chruściel for helpful discussions. Research supported in part
by the Austrian Science Fund (FWF), Project P34274 and Grant TAI 483-N, as
well as the Vienna University Research Platform TURIS. TBM is a recipient of a
DOC Fellowship of the Austrian Academy of Sciences at the Faculty of Physics
of the University of Vienna and acknowledges partial support from the Vienna
Doctoral School in Physics (VDSP).

## References

  * [1] Jerzy Plebański “Electromagnetic Waves in Gravitational Fields” In _Physical Review_ 118.5, 1960, pp. 1396–1408 DOI: 10.1103/PhysRev.118.1396
  * [2] F.. Estabrook and H.. Wahlquist “Response of Doppler spacecraft tracking to gravitational radiation.” In _General Relativity and Gravitation_ 6.5, 1975, pp. 439–447 DOI: 10.1007/BF00762449
  * [3] David Christie “Electromagnetic waves in the Kerr Schild metric” In _General Relativity and Gravitation_ 40.1, 2008, pp. 57–79 DOI: 10.1007/s10714-007-0515-2
  * [4] M. Rakhmanov “On the round-trip time for a photon propagating in the field of a plane gravitational wave” In _Classical and Quantum Gravity_ 26.15, 2009, pp. 155010 DOI: 10.1088/0264-9381/26/15/155010
  * [5] Michael J. Koop and Lee Samuel Finn “Physical response of light-time gravitational wave detectors” In _Phys. Rev. D_ 90.6, 2014, pp. 062002 DOI: 10.1103/PhysRevD.90.062002
  * [6] Arkadiusz Błaut “Gauge independent response of a laser interferometer to gravitational waves” In _Classical and Quantum Gravity_ 36.5, 2019, pp. 055004 DOI: 10.1088/1361-6382/ab01ad
  * [7] Lee Samuel Finn “Response of interferometric gravitational wave detectors” In _Phys. Rev. D_ 79.2, 2009, pp. 022002 DOI: 10.1103/PhysRevD.79.022002
  * [8] Raymond Angélil and Prasenjit Saha “Geometrical versus wave optics under gravitational waves” In _Phys. Rev. D_ 91.12, 2015, pp. 124007 DOI: 10.1103/PhysRevD.91.124007
  * [9] J.. Lobo “Effect of a weak plane GW on a light beam” In _Classical and Quantum Gravity_ 9.5, 1992, pp. 1385–1394 DOI: 10.1088/0264-9381/9/5/019
  * [10] Enrico Montanari “On the propagation of electromagnetic radiation in the field of a plane gravitational wave” In _Classical and Quantum Gravity_ 15.8, 1998, pp. 2493–2507 DOI: 10.1088/0264-9381/15/8/024
  * [11] Mirco Calura and Enrico Montanari “Exact solution to the homogeneous Maxwell equations in the field of a gravitational wave in linearized theory” In _Classical and Quantum Gravity_ 16.2, 1999, pp. 643–652 DOI: 10.1088/0264-9381/16/2/025
  * [12] J.. Codina, J. Graells and C. Martín “Electromagnetic phenomena induced by weak gravitational fields. Foundations for a possible gravitational wave detector” In _Phys. Rev. D_ 21.10, 1980, pp. 2731–2735 DOI: 10.1103/PhysRevD.21.2731
  * [13] Chan Park and Dong-Hoon Kim “Observation of gravitational waves by light polarization” In _European Physical Journal C_ 81.1, 2021, pp. 95 DOI: 10.1140/epjc/s10052-021-08893-4
  * [14] Lev Davidovich Landau and E.. Lifshitz “The Classical Theory of Fields” 2, Course of theoretical physics Pergamon Press, 1975
  * [15] Hans Stephani et al. “Exact solutions of Einstein’s field equations” Cambridge University Press, 2003 DOI: 10.1017/cbo9780511535185
  * [16] Bethan Cropp and Matt Visser “General polarization modes for the Rosen gravitational wave” In _Classical and Quantum Gravity_ 27.16, 2010, pp. 165022 DOI: 10.1088/0264-9381/27/16/165022
  * [17] Nikodem J. Popławski “A Michelson interferometer in the field of a plane gravitational wave” In _Journal of Mathematical Physics_ 47.7, 2006, pp. 072501 DOI: 10.1063/1.2212670
  * [18] Thomas B. Mieling “The response of laser interferometric gravitational wave detectors beyond the eikonal equation” In _Classical and Quantum Gravity_ 38.17, 2021, pp. 175007 DOI: 10.1088/1361-6382/ac15db
  * [19] Thomas B. Mieling, Piotr T. Chruściel and Stefan Palenta “The electromagnetic field in gravitational wave interferometers” In _Classical and Quantum Gravity_ 38.21, 2021, pp. 215004 DOI: 10.1088/1361-6382/ac2270