# Challenges of Data-Driven Simulation of Diverse and Consistent Human Driving
Behaviors

Kalle Kujanpää$\ast$$\ddagger$ Daulet Baimukashev$\dagger$ Shibei
Zhu$\ast$$\ddagger$ Shoaib Azam$\dagger$$\ddagger$ Farzeen
Munir$\dagger$$\ddagger$
Gokhan Alcan$\dagger$ Ville Kyrki$\dagger$$\ddagger$

###### Abstract

Building simulation environments for developing and testing autonomous
vehicles necessitates that the simulators accurately model the statistical
realism of the real-world environment, including the interaction with other
vehicles driven by human drivers. To address this requirement, an accurate
human behavior model is essential to incorporate the diversity and consistency
of human driving behavior. We propose a mathematical framework for designing a
data-driven simulation model that simulates human driving behavior more
realistically than the currently used physics-based simulation models.
Experiments conducted using the NGSIM dataset validate our hypothesis
regarding the necessity of considering the complexity, diversity, and
consistency of human driving behavior when aiming to develop realistic
simulators.

## Introduction

Modeling human behavior is pivotal in building intelligent learning systems,
offering profound insights into how humans interact, behave, and respond
within complex environments. This understanding is particularly crucial in the
context of autonomous vehicle (AV) simulation (Feng et al. 2021). Modeling and
simulating human behaviors for building a data-driven simulation environment
for autonomous vehicles is challenging due to human behaviors’ high
variability and diversity (Fuchs, Passarella, and Conti 2022; Chen et al.
2023). As the simulation provides a cost-effective way to evaluate the
performance of autonomous vehicles, achieving statistical realism in such a
simulation environment remains a persistent challenge in the field,
necessitating sophisticated modeling techniques to incorporate human
behaviors.

A key component of creating a high-fidelity driving environment is precisely
modeling human behaviors. There have been notable efforts in the literature to
develop autonomous driving simulators, propelled by rapid advancements in
artificial intelligence, computer graphics and vision, and high-performance
computing devices. Simulators like CARLA (Dosovitskiy et al. 2017), Airsim
(Shah et al. 2018), DriveSim 111https://developer.nvidia.com/drive/simulation,
Wayve Infinity Simulator 222https://wayve.ai/thinking/introducing-wayve-
infinity-simulator/, and SimulationCity
333https://waymo.com/blog/2021/06/SimulationCity.html provide models for
accurate vehicle dynamics and photo-realistic scenarios. However, these
simulators mainly focus on vehicle dynamics and often overlook the driving
environment, particularly the behavior of non-player characters, such as
background road users. The behavior of these background vehicles is typically
either replayed from logs or simulated using heuristic rules, creating a
disparity between the simulation and the real-world driving environment.
Furthermore, the microscopic simulation environments mimic the interaction
between agents using physics-driving models, and hand-crafted rules such as
car-following (Treiber, Hennecke, and Helbing 2000), gap-acceptance
(Mahmassani and Sheffi 1981), and lane-changing models (Kesting, Treiber, and
Helbing 2007; Erdmann 2015) often lack the diversity to model the human
behavior for simulated agents. Therefore, these physics and theory-based
driving models can hardly be generalized and scaled to capture the underlying
diversity of human behavior for driving.

To this end, in this work, we focus our studies on building a high-fidelity
simulation environment statistically representative of the real-world driving
environment by modeling human driving behaviors. We aim to develop a data-
driven simulation model that captures the subtle complexities and the wide
range of variability inherent in human driving behavior in various scenarios
and conditions. In this work, we raise the following research questions open
to discussion:

  1. 1.

How can we accurately model the inherently unpredictable and stochastic nature
of real human driving behavior?

  2. 2.

Considering the diverse nature of human driving behavior, how can this
diversity and variability, especially in different human driving styles, be
effectively incorporated into the simulation environment?

  3. 3.

How to incorporate the temporal consistency of human driving behavior in
simulations?

## Related Work

Human driving behavior is extensively studied, encompassing a broad range of
problem definitions, assumptions in modeling, and diverse methodologies
related to modeling human driving behavior (Brown, Driggs-Campbell, and
Kochenderfer 2020). Driving is a widespread yet intricate dynamic activity
that requires the integration of physiological, mental, and physical aspects
in managing the vehicle to achieve a specific objective. It involves a
repetitive interplay between the driver, the vehicle, and the traffic
environment. The overall task of driving encompasses a wide range of
physiological constructs, behaviors, and action patterns. Modeling the complex
process of how humans execute diverse road maneuvers is challenging. The
variability and diversity in human drivers add to the complexity of developing
Driver Behavior Models (DBMs) (Negash and Yang 2023).

Most classical DBMs such as the Intelligent Driver Model (IDM) (Treiber,
Hennecke, and Helbing 2000), the MOBIL model (Kesting, Treiber, and Helbing
2007), and models by Gipps (Gipps 1981, 1986) are deterministic. These models
often fail to encapsulate the stochastic nature inherent in high-fidelity
human driving behavior. In these models, the parameters are static during
simulations, and the noise is typically disregarded. Despite the calibration
of parameter values with real-world data (Sangster, Rakha, and Du 2013; Li,
Chen, and Zhang 2016; Punzo, Ciuffo, and Montanino 2012; Ma and Qu 2020),
these deterministic models struggle to accurately reflect the variability and
diversity in human driving behavior.

Recently, researchers have shifted their focus towards employing machine
learning-based methods (Xie et al. 2019; Huang, Sun, and Sun 2018; Wang et al.
2017), particularly neural networks, for fitting the DBM. Large naturalistic
driving datasets are utilized to estimate the parameters that reflect observed
trajectories of human drivers. However, they still lack accurate
stochasticity. Introducing the noise into the model can integrate a degree of
stochasticity. Typically, this noise is modeled as Gaussian noise (Kuefler et
al. 2017). However, Gaussian noise may not adequately capture the variability
and diversity inherent in human driving behavior. This is because the
interactions in various driving conditions are highly complex and do not
always conform to the simple Gaussian distribution (Laval, Toth, and Zhou
2014; Treiber and Kesting 2017). In (Yang and Peng 2010), extreme-value and
log-normal distributions are applied. However, the driving behavior exhibited
by humans in different scenarios can vary significantly, thus necessitating a
more sophisticated approach to modeling noise that accurately reflects this
complexity.

In addition to incorporating external noise, model parameters can also be used
to introduce stochasticity. Yang and Peng (2010) accounted for human errors in
car-following by modeling time delay and distraction intervals using
exponential and lognormal distributions, respectively. Game theory was applied
in Talebpour, Mahmassani, and Hamdar (2015) for stochastic lane-changing
behavior modeling. Hamdar, Mahmassani, and Treiber (2015) developed a utility-
based stochastic car-following model with acceleration probabilities following
a continuous logit model. However, these methods often overlook the
distributional accuracy of driving behaviors post-introduction of
stochasticity. Exceptions are found in Wang et al. (2009) and Chen, Li, and
Zhang (2010), where stochastic car-following models were proposed to
accurately capture the distribution of time headway accurately, aligning
simulation environments with real-world data. Nevertheless, these studies
mainly focus on single-lane road scenarios and do not address error
accumulation, which could significantly distort long-term simulations.

Despite various related studies, most existing methodologies fall short in
accurately capturing the distribution of stochastic human driving behaviors.
Consequently, the driving environments generated by these behavior models are
not sufficiently precise in replicating the statistical characteristics of
real-world simulation environments, an aspect crucial for autonomous vehicle
(AV) testing.

## Problem Formulation

Given a set of $n$ demonstrations collected from the real world as
$\mathcal{D}=\\{\tau_{0},\tau_{1},\dots,\tau_{n}\\}$, where each trajectory
$\tau_{i}$ is given as a temporal sequence of state and action pairs
$\tau_{i}=\\{(s_{0},a_{0}),\dots,(s_{T},a_{T})\\}$, the driver behavior
modeling problem can be formulated as follows:

min
$\displaystyle\mathop{\mathbb{E}}_{\tau_{i}\sim\mathcal{D}}\Bigg{(}J(\tau_{i},\hat{\tau}_{i})\Bigg{)}$
(1) subject to
$\displaystyle\hat{\tau}_{i}=\\{(s_{0},\hat{a}_{0}),\dots,(s_{T},\hat{a}_{T})\\},$
$\displaystyle\hat{a}_{t}\sim\pi_{\theta}(s_{t},\psi_{i}),$
$\displaystyle\psi_{i}\sim p(\psi).$

Here, $J(\tau_{i},\hat{\tau}_{i})$ is an objective function measuring the
deviation between an actual trajectory $\tau_{i}$ and a predicted trajectory
$\hat{\tau}_{i}$. $\psi$ represents the parameters of the driver behavior
model, drawn from some distribution $p(\psi)$, and $\theta$ stands for the
general parameters of the policy $\pi$.

Eq.1 represents the general form of modeling the human driving behavior for
the simulated environment. Following this, two possible types of solutions
arise that address the underlying question of how to model the human driving
behavior: theory-based and data-driven models.

In theory-based methods, the policy parameters $\theta$ are generally not
optimized. Rather, they are implicitly determined by the chosen driving model,
such as IDM or its enhancements. Instead, to simulate human driving, the
objective function $J(\tau_{i},\hat{\tau}_{i})$ can be minimized by learning
the distribution $p(\psi)$.

Alternatively, when utilizing generic data-driven policies that are not
structured as theory-based methods, such neural networks, there is an
enhancement in generalization capabilities. In such instances, actions are
often drawn from a policy directly parameterized by $\theta$, as follows:

$\hat{a}_{t}\sim\pi_{\theta}(s(t),\psi)$ (2)

For each time step $t$, the predicted action $\hat{a}_{t}(\theta)$ is
determined by the policy function $\pi_{\theta}$, parameterized by $\theta$,
based on the current state $s(t)$ and a set of parameters $\psi$. These
parameters $\psi$ play a pivotal role in the model, as they can be either
predefined priors or latent variables. This distinction allows the model to
either leverage existing knowledge about driving behaviors (when $\psi$ is a
prior) or to infer and optimize these parameters based on the observed data
(when $\psi$ is treated as a latent variable). We argue that treating the
parameters $\psi$ correctly is essential for accurately modeling the diversity
of actual human behavior, and they should be driver-specific. In the following
case study, we compare a theory-based model, IDM, to a generic data-driven
model in terms of accurately replicating human driving behavior.

## Case Study

Figure 1: We compare the velocity predictions of the two types of IDMs and a
learned XGBoost predictor to the ground truth data in the simulated NGSIM
environments. We find that replicating the human driving behavior is very
difficult, even in this simplistic car-following setting.

b]

Figure 2: We compare the change in velocity predicted by XGBoost to the ground truth change in velocity. The plot indicates that human behavior may be inherently unpredictable, and capturing the behavior seems to be highly challenging. Dataset | Parameters | MSE | SD
---|---|---|---
US 101 | Driver-specific | 215.63 $\pm$ 3.08 | 240.55
| Shared | 289.44 $\pm$ 3.16 | 246.58
I-80 | Driver-specific | 220.16 $\pm$ 5.15 | 387.92
| Shared | 301.19 $\pm$ 5.82 | 438.52
Peachtree | Driver-specific | 58.50 $\pm$ 2.48 | 119.64
| Shared | 107.49 $\pm$ 3.27 | 157.66
Lankershim | Driver-specific | 264.19 $\pm$ 6.08 | 300.14
| Shared | 416.55 $\pm$ 6.14 | 303.25
Table 1: The mean squared error plus-minus standard error of the MSE for IDM
models fitted to the NGSIM data. The results show that fitting driver-specific
models outperforms using shared models, but the performance of both variants
is highly suboptimal. Figure 3: We plot the average acceleration and
deceleration and the preferred minimum time headway, that is, the safety
margin between two vehicles for each driver in each of the four NGSIM
datasets. The plots show that human driving behavior is highly diverse. It can
be fat-tailed and multi-modal and does not simply follow a Gaussian
distribution. Figure 4: We plot the distribution of the IDM parameters for
each of the four NGSIM datasets when each driver is represented by a separate
set of parameters. The black vertical line can be interpreted as the human
mean, whereas the red-shaded area depicts the expected noise in the IDM
fitting process. Figure 5: We analyze the consistency of human drivers on the
four datasets in the NGSIM data by fitting IDM to the longest trajectories in
the dataset. Our results show that in all of the NGSIM datasets, the human
driving style shows statistically significant consistency when we focus on the
human drivers of which most data has been collected. The uncertainty plotted
with error bars is one standard error of the mean.

In this section, we present a series of experiments designed to highlight the
importance and relevancy of the research questions outlined in the
introduction. We discuss the modeling of human driving behavior for its
realistic simulation and argue that human behavior is complex and diverse but
also consistent. Utilizing the Next Generation Simulation (NGSIM) dataset as
our foundation, we focus on the application of the Intelligent Driver Model
(IDM) to predict the diverse driving behaviors observed in the dataset.
Additionally, we scrutinize the limitation of a more robust model, XGBoost, in
capturing realistic driving behavior. The results not only highlight the
inherent diversity and stochastic nature of human driving styles but also
underscore the importance of employing sophisticated methods capable of
robustly addressing both the complexities and consistencies inherent in human
behavior.

### Background

#### Intelligent Driver Model

While modeling complex human driving is challenging, researchers have proposed
various approaches to tackle, for instance, car-following and lane-changing
behaviors. One prominent model in this domain is the Intelligent Driver Model
(IDM), initially introduced by Treiber, Hennecke, and Helbing (2000). Widely
accepted in the literature, IDM is a time-continuous car-following model
focusing on the interaction between vehicles in a traffic flow through
mathematical formulations of vehicle acceleration. Notably, IDM has found
practical applications in traffic management systems, autonomous vehicle
development, and urban planning, and several modifications and improvements
have been proposed over the years (Li et al. 2015; Derbel et al. 2013;
Kesting, Treiber, and Helbing 2010).

IDM describes the dynamics of a vehicle $\alpha$ by updating its location
$x_{\alpha}$ and velocity $v_{\alpha}$ at time $t$ according to the following
two ordinary differential equations:

$\dot{x}_{\alpha}=\frac{dx_{\alpha}}{dt}=v_{\alpha}$ (3)
$\dot{v}_{\alpha}=\frac{dv_{\alpha}}{dt}=a\left[1-\left(\frac{v_{\alpha}}{v_{0}}\right)^{\delta}-\left(\frac{s^{\ast}(v_{\alpha},\Delta
v_{\alpha})}{s_{\alpha}}\right)^{2}\right],$ (4)

where

$\displaystyle s^{\ast}(v_{\alpha},\Delta
v_{\alpha})=s_{0}+v_{\alpha}T+\frac{v_{\alpha}\Delta v_{\alpha}}{2\sqrt{ab}}.$

$\Delta v_{\alpha}$ stands for the approaching rate. $a$ is the desired
acceleration, $v_{0}$ the desired speed, $s^{\ast}$ the spacing between the
ego and leading vehicle, $s_{0}$ the minimum spacing, $b$ the comfortable
braking deceleration, $T$ desired time headway, and $\delta$ a constant.

#### NGSIM Dataset

To study human driving behaviors, we leverage the rich and publicly available
real-world dataset from the Next Generation Simulation (NGSIM) project (U.S.
Department Of Transportation Federal Highway Administration 2017). The dataset
contains data collected from four locations in the US: US 101, I-80,
Lankershim Boulevard, and Peachtree Street. Unlike other real-world driving
datasets, such as those from CARLA, Waymo, and Lyft, which often involve data
collection by a limited number of drivers, the NGSIM dataset provides
extensive data for thousands of drivers, achieving a much more extensive
coverage of human behavior. The vehicle trajectories in the NGSIM dataset have
been recorded with a 100 ms time resolution. Each record within the dataset
includes the vehicle position in both a local and global coordinate system,
lane number, velocity, acceleration, vehicle type, width, and length, and the
vehicle ID of the preceding and following vehicle.

### Complexity of Human Driving Behavior

We hypothesize that real human driving behavior is inherently unpredictable
and stochastic. The limitations of fitting a parametric model, such as the
IDM, to human data are evident in its inability to capture the complexity of
the behavior. Even more performant predictive models, like the tree-boosting
algorithm XGBoost (Chen and Guestrin 2016), face significant challenges in
this task. To evaluate this hypothesis, we applied both the IDM and XGBoost to
the NGSIM dataset. Initially, we individually fitted an IDM to each driver in
the dataset. Subsequently, we fitted a shared IDM to each of the four
scenarios in the dataset.

The fitting process of the IDM to the NGSIM dataset involves several key
steps. First, we organize the dataset into trajectories, with each trajectory
corresponding to an individual driver. Notably, if a driver changes lanes, the
trajectory is split into two distinct trajectories. For each trajectory, we
extract the initial velocity and the gap to the vehicle in front.
Additionally, we include the velocity of the preceding driver as an input to
the model. Subsequently, we simulate the acceleration using the IDM equations
(see Equations 4-3) and update the gap to the leading vehicle based on the new
velocity and the velocity of the other driver. The simulation occurs with a
time step of $dt=0.1$, corresponding to the duration between frames in NGSIM.
The process is akin to solving the IDM ordinary differential equation using
the Euler method, and this is repeated for each frame within the trajectory.
The optimization objective is to minimize the mean squared error between the
predicted velocities and ground truth velocities. The IDM parameters are
optimized using Optuna (Akiba et al. 2019), with 500 trials per optimization.
The parameters subjected to optimization include the desired velocity $v_{0}$,
minimum spacing $s_{0}$, desired time headway $T$, acceleration $a$, and the
comfortable braking deceleration $b$. The exponent $\gamma$ is held constant
at four.

The final mean squared errors from our experiment using the NGSIM dataset are
presented in Table 1. As expected, fitting an IDM for each driver outperforms
the use of shared IDM parameters across the trajectories, with the observed
difference being statistically significant. However, it is crucial to
highlight the remarkably high standard deviation of the mean squared errors,
underscoring the suboptimal nature of the IDM fit.

To qualitatively illustrate the suboptimality of the IDM fit, Figure 1
showcases randomly selected trajectories from each NGSIM dataset. The initial
states were sampled from the initial states of the trajectories in the NGSIM
dataset, and the trajectory was simulated using the learned IDM parameters.
Both the trajectories simulated with driver-specific parameters (in green) and
those with shared IDM parameters (in orange) fail to match the ground truth
data in blue.

Moreover, the intricacy of predicting human driving behavior is further
emphasized by the limitations of the tree-boosting algorithm XGBoost in this
task (Chen and Guestrin 2016). XGBoost was fitted to the NGSIM data by
learning to predict the change in velocity based on the current velocity, the
gap to the preceding vehicle in front in seconds and meters, and the velocity
difference. No driver-specific variables were included in the training data.
Despite training the model for 20,000 tree-boosting iterations, when the
training loss had stabilized, Figure 2 illustrates that the algorithm fails to
achieve a low one-step prediction error, even in the supervised setting.
Consequently, deploying the learned XGBoost driving policy in the NGSIM
environment proves ineffective in predicting realistic driving behaviors, and
this issue is exacerbated by accumulating errors, as depicted by the red
curves in Figure 1.

It is essential to note that the selected methods, IDM and XGBoost, even
struggle to predict the training data, emphasizing the complexity of
replicating human behavior. However, it must be acknowledged that the
difficulty may be influenced by the lower quality of the NGSIM data that is
extracted from raw video files and may contain anomalies. In summary, we
contend that sophisticated methods capable of accommodating the stochasticity
and complexity of human behavior are indispensable for achieving realistic and
accurate human simulations.

### Diversity of Human Driving Behavior

We posit that human drivers exhibit unique driving styles, characterized by
variations in, for instance, aggression, composure, cornering preferences, and
turning behavior. The spectrum of driving styles is extensive. Achieving
accuracy and realism in simulations of human driving behavior requires
considering and modeling this diversity. Furthermore, we assert that simply
sampling, for instance, IDM parameters from a fitted Gaussian may prove
inadequate for capturing the nuanced richness and complexity of human driving
styles.

To rigorously assess our hypothesis regarding the diversity of human driving
styles, we conducted a careful examination of the NGSIM dataset, which
provides valuable insights into the behavior of human drivers. Since the NGSIM
dataset primarily captures car-following behavior, our initial analysis
focuses on three key variables: the average acceleration, the deceleration
during braking, and the preferred headway or safety margin to the vehicle in
front.

In Figure 3, we present the distribution of these variables across the four
NGSIM datasets using histograms. The first row of subplots depicts the
distribution of average acceleration. To derive this, we computed the change
in velocity between two consecutive frames for each driver, including pairs of
frames where the velocity increased by at least two meters per second. Then,
we computed the mean acceleration and plotted the result. Notably, the
Peachtree dataset exhibited a lower percentage (approximately 50 %) of drivers
with such accelerations, in contrast to the US 101, I-80, and Lankershim
datasets, where at least 98 % of the drivers demonstrated these types of
accelerations. This discrepancy should be considered in the interpretation of
our plots. The Peachtree and Lankershim datasets show evidence that the
distribution of mean accelerations resembles a multimodal one, whereas the
histograms for US 101 and I-80 are unimodal.

The second row of Figure 3 displays the distribution of mean decelerations,
following a similar process. Pairs of frames where the speed decreased by at
least two meters per second were filtered, and the average deceleration for
each driver was computed and plotted.

The final row presents plots of the minimum time headway for each driver. This
was determined by dividing the distance to the vehicle in front by the current
velocity of the ego vehicle, with the minimum value taken for each trajectory.
Drivers with a minimum headway of over five seconds were excluded from our
analysis due to insufficient traffic. Notably, anomalies in the dataset may
contribute to apparent peaks at zero seconds, but we believe these are not
representative of actual human behavior.

The modeling of human driving behavior often employs parametric systems such
as IDM. In the second part of our analysis of driving style diversity, we
fitted an IDM for each individual driver and compared the distribution of IDM
parameters across the set of human drivers. Figure 4 illustrates the
distribution of the five IDM parameters for each of the four NGSIM datasets.

In the visual representation, the black vertical line signifies the
parameters’ values when a single model is fitted to every trajectory in the
dataset, essentially representing the human mean. However, it’s also crucial
to account for the inherent noise in fitting the IDM parameters. To estimate
this noise, we conducted multiple fittings of IDM for each diver and computed
the average standard deviation of the parameter values. The red-shaded area in
the figure represents the mean plus-minus two standard deviations. If the
distribution of IDM parameters adhered to a Gaussian distribution, and each
human driver could be adequately represented with the same set of parameter
values, we would expect 95 % of the values to fall within the shaded region.
However, this is not the case, indicating significant diversity in the driving
styles.

It is important to note that, as demonstrated in the previous section, the fit
of IDM to the data is suboptimal, particularly for the maximum acceleration
parameter (a), as illustrated also by Figure 4. This limitation must be
considered when drawing conclusions from the experiment.

In summary, our findings provide compelling evidence of the diverse nature of
human driving styles. When modeling these tyles using parametric models, our
results reveal that the distribution of parameters can exhibit multimodality
or fat-tailed characteristics. This implies that simplistic parameter sampling
from a Gaussian distribution is unrealistic and should be eschewed. Instead,
in realistic simulations, it is advisable to employ more sophisticated
sampling strategies to better capture the nuances of human driving behavior.

### Consistency of Human Driving Behavior

We hypothesize that human driving behavior is consistent, and this consistency
should be accounted for in realistic simulations of human driving. Each human
driver has a specific style and will mostly adhere to that, although within
the constraints of the natural variability of human behavior. This consistency
will also be seen during trajectories; a human will drive similarly at the
beginning of a trajectory as at the end. To evaluate this hypothesis, we
conducted an experiment with the NGSIM dataset and IDM. We select trajectories
from all four datasets in the NGSIM data and split them into half. Then, we
fit an IDM to the first and second halves of the trajectory as described in
the earlier sections of this paper. Then, we compare if the optimization
procedure found the same set of parameters for the two halves. We do this by
calculating the Euclidean distance (L2-norm) between the IDM parameters.

Figure 5 shows the results for the human consistency experiment. The red curve
in the subfigures illustrates the inherent randomness in fitting the IDM
parameters. If the parameters are optimized for each trajectory twice, the
mean Euclidean distance between the found IDM parameters is approximately 15.
The green curve shows the mean Euclidean distance between IDM parameters for
two different drivers, which is, as expected, significantly higher than the
distance between the parameters for the same drivers. The blue curves in the
graphs show that the Euclidean distance between the IDM parameters for the two
halves of the same trajectory is lower than the distance between the IDM
parameters for different drivers, especially in the case of longer
trajectories, which implies that the IDM parameters capture the individual
driving style of each driver and, moreover, that this style is consistent over
time. If we focus on the longest trajectories in the dataset, this difference
in distances is statistically significant. However, this consistency is not
visible for shorter trajectories. We believe that it is due to the
trajectories being too short and the noise of fitting the IDM parameters
starting to dominate the results. Our results support our hypothesis of human
driving behavior being temporally consistent, and therefore, it should be
incorporated in realistic simulations of human driving behavior.

## Conclusion

In this work, we aim to address the importance of human behavior modeling with
specific applications to build high-fidelity data-driven simulators for
autonomous driving. We posited the development of a data-driven simulation
environment as a learning problem with open-ended solutions focused on theory-
based and data-driven models. Our investigation into the complexity,
diversity, variability, and temporal consistency of human driving behavior is
grounded in empirical analysis using the real-world NGSIM dataset.

Our key findings reveal that human behaviors are complex, highlighting that
theory-based models, like the Intelligent Driver Model (IDM), and data-driven
models that fail to explicitly incorporate human driving styles, such as
XGBoost, encounter limitations in fully capturing this complexity. Similarly,
simplistic parameter sampling may not capture the subtle diversity of human
behavior. In addition, our experiments showcased that human driving behavior
exhibits temporal consistency, meaning individual driving styles remain stable
over time. This consistency is evident from the analysis of IDM parameters
across different trajectory halves, emphasizing the need to incorporate this
aspect into simulations.

This research contributes to understanding human driving behavior in the
context of autonomous vehicle simulation and paves the way for future
explorations. It opens up avenues for formalizing data-driven models to
emulate human behavior more effectively and sets the stage for developing
models that capture human driving patterns’ inherent variability and
diversity. These advancements are crucial for creating better human models for
simulation.

## References

  * Akiba et al. (2019) Akiba, T.; Sano, S.; Yanase, T.; Ohta, T.; and Koyama, M. 2019.  Optuna: A next-generation hyperparameter optimization framework.  In _Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery & data mining_, 2623–2631.
  * Brown, Driggs-Campbell, and Kochenderfer (2020) Brown, K.; Driggs-Campbell, K.; and Kochenderfer, M. J. 2020.  Modeling and prediction of human driver behavior: A survey.  _arXiv preprint arXiv:2006.08832_ , 4–45.
  * Chen et al. (2023) Chen, D.; Zhu, M.; Yang, H.; Wang, X.; and Wang, Y. 2023.  Data-driven Traffic Simulation: A Comprehensive Review.  _arXiv preprint arXiv:2310.15975_.
  * Chen and Guestrin (2016) Chen, T.; and Guestrin, C. 2016.  Xgboost: A scalable tree boosting system.  In _Proceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining_ , 785–794.
  * Chen, Li, and Zhang (2010) Chen, X.; Li, L.; and Zhang, Y. 2010.  A Markov model for headway/spacing distribution of road traffic.  _IEEE Transactions on Intelligent Transportation Systems_ , 11(4): 773–785.
  * Derbel et al. (2013) Derbel, O.; Peter, T.; Zebiri, H.; Mourllion, B.; and Basset, M. 2013.  Modified Intelligent Driver Model for driver safety and traffic stability improvement.  _IFAC Proceedings Volumes_ , 46(21): 744–749.  7th IFAC Symposium on Advances in Automotive Control.
  * Dosovitskiy et al. (2017) Dosovitskiy, A.; Ros, G.; Codevilla, F.; Lopez, A.; and Koltun, V. 2017.  CARLA: An open urban driving simulator.  In _Conference on robot learning_ , 1–16. PMLR.
  * Erdmann (2015) Erdmann, J. 2015.  SUMO’s lane-changing model.  In _Modeling Mobility with Open Data: 2nd SUMO Conference 2014 Berlin, Germany, May 15-16, 2014_ , 105–123. Springer.
  * Feng et al. (2021) Feng, S.; Yan, X.; Sun, H.; Feng, Y.; and Liu, H. X. 2021.  Intelligent driving intelligence test for autonomous vehicles with naturalistic and adversarial environment.  _Nature communications_ , 12(1): 748.
  * Fuchs, Passarella, and Conti (2022) Fuchs, A.; Passarella, A.; and Conti, M. 2022.  Modeling Human Behavior Part I–Learning and Belief Approaches.  _arXiv preprint arXiv:2205.06485_.
  * Gipps (1981) Gipps, P. G. 1981.  A behavioural car-following model for computer simulation.  _Transportation research part B: methodological_ , 15(2): 105–111.
  * Gipps (1986) Gipps, P. G. 1986.  A model for the structure of lane-changing decisions.  _Transportation Research Part B: Methodological_ , 20(5): 403–414.
  * Hamdar, Mahmassani, and Treiber (2015) Hamdar, S. H.; Mahmassani, H. S.; and Treiber, M. 2015.  From behavioral psychology to acceleration modeling: Calibration, validation, and exploration of drivers’ cognitive and safety parameters in a risk-taking environment.  _Transportation Research Part B: Methodological_ , 78: 32–53.
  * Huang, Sun, and Sun (2018) Huang, X.; Sun, J.; and Sun, J. 2018.  A car-following model considering asymmetric driving behavior based on long short-term memory neural networks.  _Transportation research part C: emerging technologies_ , 95: 346–362.
  * Kesting, Treiber, and Helbing (2007) Kesting, A.; Treiber, M.; and Helbing, D. 2007.  General lane-changing model MOBIL for car-following models.  _Transportation Research Record_ , 1999(1): 86–94.
  * Kesting, Treiber, and Helbing (2010) Kesting, A.; Treiber, M.; and Helbing, D. 2010.  Enhanced intelligent driver model to access the impact of driving strategies on traffic capacity.  _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 368(1928): 4585–4605.
  * Kuefler et al. (2017) Kuefler, A.; Morton, J.; Wheeler, T.; and Kochenderfer, M. 2017.  Imitating driver behavior with generative adversarial networks.  In _2017 IEEE Intelligent Vehicles Symposium (IV)_ , 204–211. IEEE.
  * Laval, Toth, and Zhou (2014) Laval, J. A.; Toth, C. S.; and Zhou, Y. 2014.  A parsimonious model for the formation of oscillations in car-following models.  _Transportation Research Part B: Methodological_ , 70: 228–238.
  * Li, Chen, and Zhang (2016) Li, L.; Chen, X. M.; and Zhang, L. 2016.  A global optimization algorithm for trajectory data based car-following model calibration.  _Transportation Research Part C: Emerging Technologies_ , 68: 311–332.
  * Li et al. (2015) Li, Z.; Li, W.; Xu, S.; and Qian, Y. 2015.  Stability analysis of an extended intelligent driver model and its simulations under open boundary condition.  _Physica A: Statistical Mechanics and its Applications_ , 419: 526–536.
  * Ma and Qu (2020) Ma, L.; and Qu, S. 2020.  A sequence to sequence learning based car-following model for multi-step predictions considering reaction delay.  _Transportation research part C: emerging technologies_ , 120: 102785\.
  * Mahmassani and Sheffi (1981) Mahmassani, H.; and Sheffi, Y. 1981.  Using gap sequences to estimate gap acceptance functions.  _Transportation Research Part B: Methodological_ , 15(3): 143–148.
  * Negash and Yang (2023) Negash, N. M.; and Yang, J. 2023.  Driver Behavior Modeling Towards Autonomous Vehicles: Comprehensive Review.  _IEEE Access_.
  * Punzo, Ciuffo, and Montanino (2012) Punzo, V.; Ciuffo, B.; and Montanino, M. 2012.  Can results of car-following model calibration based on trajectory data be trusted?  _Transportation research record_ , 2315(1): 11–24.
  * Sangster, Rakha, and Du (2013) Sangster, J.; Rakha, H.; and Du, J. 2013.  Application of naturalistic driving data to modeling of driver car-following behavior.  _Transportation research record_ , 2390(1): 20–33.
  * Shah et al. (2018) Shah, S.; Dey, D.; Lovett, C.; and Kapoor, A. 2018.  Airsim: High-fidelity visual and physical simulation for autonomous vehicles.  In _Field and Service Robotics: Results of the 11th International Conference_ , 621–635. Springer.
  * Talebpour, Mahmassani, and Hamdar (2015) Talebpour, A.; Mahmassani, H. S.; and Hamdar, S. H. 2015.  Modeling lane-changing behavior in a connected environment: A game theory approach.  _Transportation Research Procedia_ , 7: 420–440.
  * Treiber, Hennecke, and Helbing (2000) Treiber, M.; Hennecke, A.; and Helbing, D. 2000.  Congested traffic states in empirical observations and microscopic simulations.  _Physical review E_ , 62(2): 1805.
  * Treiber and Kesting (2017) Treiber, M.; and Kesting, A. 2017.  The intelligent driver model with stochasticity-new insights into traffic flow oscillations.  _Transportation research procedia_ , 23: 174–187.
  * U.S. Department Of Transportation Federal Highway Administration (2017) U.S. Department Of Transportation Federal Highway Administration. 2017.  Next Generation Simulation (NGSIM) Vehicle Trajectories and Supporting Data.
  * Wang et al. (2009) Wang, F.; Li, L.; Hu, J.-M.; Ji, Y.; Ma, R.; and Jiang, R. 2009.  A markov-process inspired ca model of highway traffic.  _International Journal of Modern Physics C_ , 20(01): 117–131.
  * Wang et al. (2017) Wang, X.; Jiang, R.; Li, L.; Lin, Y.; Zheng, X.; and Wang, F.-Y. 2017.  Capturing car-following behaviors by deep learning.  _IEEE Transactions on Intelligent Transportation Systems_ , 19(3): 910–920.
  * Xie et al. (2019) Xie, D.-F.; Fang, Z.-Z.; Jia, B.; and He, Z. 2019.  A data-driven lane-changing model based on deep learning.  _Transportation research part C: emerging technologies_ , 106: 41–60.
  * Yang and Peng (2010) Yang, H.-H.; and Peng, H. 2010.  Development of an errorable car-following driver model.  _Vehicle System Dynamics_ , 48(6): 751–773.


# Wave packets in QFT: leading order width corrections to decay rates and
clock behaviour under Lorentz boosts

Ariel Edery<EMAIL_ADDRESS>
###### Abstract

Decay rates in quantum field theory (QFT) are typically calculated assuming
the particles are represented by momentum eigenstates (i.e. plane waves).
However, strictly speaking, localized free particles should be represented by
wave packets. This yields width corrections to the decay rate and to the clock
behaviour under Lorentz boosts. We calculate the decay rate of a particle of
mass $M$ modeled as a Gaussian wavepacket of width $a$ and centered at zero
momentum. We find the decay rate to be
$\Gamma_{0}\big{[}1-\frac{3a^{2}}{4M^{2}}+\mathcal{O}\big{(}\tfrac{a^{4}}{M^{4}}\big{)}\big{]}$
where $\Gamma_{0}$ is the decay rate of the particle at rest treated as a
plane wave. The leading correction is then of order $\tfrac{a^{2}}{M^{2}}$. We
then perform a Lorentz boost of velocity $v$ on the above Gaussian and find
that its decay rate does not decrease exactly by the Lorentz factor
$\sqrt{1-v^{2}}$. There is a correction of order $\tfrac{a^{2}v^{2}}{M^{2}}$.
Therefore, the decaying wave packet does not act exactly like a typical clock
under Lorentz boosts and we refer to it is a “WP clock” (wave packet clock). A
WP clock does not move with a single velocity relative to an observer but has
a spread in velocities (more specifically, a spread in momenta). Nonetheless,
it is best viewed as a single clock as the wave packet represents a one-
particle state in QFT. WP clocks do not violate Lorentz symmetry and are not
based on new physics: they are a consequence of the combined requirements of
special relativity, quantum mechanics and localized free particles.

##  1 Introduction

Localized particles in Quantum Field Theory (QFT) are represented by wave
packets[1, 2] and localized particles are what take part in decay and
scattering processes111In this work we will focus on decay rates of wave
packets and not scattering cross sections.. In fact, replacing wave packets by
momentum eigenstates (i.e. plane waves) leads to a well-known issue at a
theoretical level: momentum eigenstates cannot be properly normalized in an
infinite volume. This is usually handled using the “finite-volume trick” [2]
where one uses finite (but large) spatial volume $\Omega$ and interactions
over a finite (but large) time interval $T$. One can then calculate the decay
rate using momentum eigenstates because the divergences that naturally occur
due to the squaring of delta functions are now regularized by $\Omega$ and
$T$. In the end, this regularization dependence cancels and does not appear in
quantities of interest. It is these calculations involving momentum
eigenstates that are usually compared to decay experiments. This works in
practice because the wave packets that appear in experiments have an energy
width which is small compared to the particle decay width [2]. The detectors
typically do not have the resolution to measure its effect. As pointed out in
[1], “Real detectors have finite resolution…the measurement of the final-state
momentum is not of such high quality that it can resolve the small variation
of this momentum that results from the momentum spread of the initial wave
packet”. As stated in [2], “the wave packet is classical in the sense that the
resolution of initial position and momentum measurements are much too large to
push the limits of the uncertainty relations”. Therefore, S-matrix
calculations carried out in the plane wave limit have been usually adequate in
practice for comparison with experiment. Nonetheless, we mention at the end of
this introduction that the effects of wave packets may have already appeared
in neutrino oscillation experiments. So the gap between experiment and wave
packet analysis may already be showing signs of narrowing, at least for a
certain type of experiment.

In this work, we obtain analytical results for the leading order correction to
decay rates due to the width of a Gaussian wave packet. There are at least
three reasons to undertake such a calculation. The first reason is better
theoretical understanding: it is beneficial to extract a leading order
analytical result from an otherwise intractable integral that has only a
numerical solution. The second reason is that the resolution of detectors and
the design of experiments keeps improving so that the effects of wave packets
on decay rates may become detectable in next generation particle decay
experiments. The third reason might be the most important: we show that wave
packets behave like a slightly different kind of clock under Lorentz boosts
and we refer to it as a “WP clock” (wave packet clock). This is interesting in
and of itself and may spawn new concepts which we discuss in the conclusion.

We first determine the leading order width correction to the decay rate of a
particle at rest due to the width $a$ of a Gaussian wave packet centered at
zero momentum. The correction is of order $\tfrac{a^{2}}{M^{2}}$ i.e. of order
of the square of the ratio of the Compton wavelength to the inverse width of
the Gaussian. We then perform a Lorentz boost of velocity $v$ on the original
Gaussian wave packet and find that the decay rate does not decrease exactly by
the usual Lorentz factor $\sqrt{1-v^{2}}$. There is a small correction of
order $\tfrac{a^{2}\,v^{2}}{M^{2}}$. This implies that the wave packet does
not act like a typical clock under Lorentz boosts and we refer to it as a “WP
clock” (for wave packet clock). The distinctive feature of a WP clock is that
it does not move with a single velocity relative to an observer but has a
spread in velocities (more specifically, a spread in momenta). It should
nonetheless be viewed as a single clock as the wave packet represents a one-
particle state in QFT. It should be emphasized that there is no breakdown of
Lorentz symmetry here. The WP clock does not stem from new physics but is an
important consequence of combining the requirements of special relativity,
quantum mechanics and localized free particles.

We show that a Gaussian wave packet does not transform into a Gaussian under a
Lorentz boost. We therefore calculate separately the decay rate of a Gaussian
wave packet centered at velocity $v$ (i.e. momentum centered at $M\,\gamma\,v$
where $\gamma=(1-v^{2})^{-1/2}$) and compare it to the decay rate of a
Gaussian centered at zero momentum. Again, we find that the two decay rates
are not related exactly by the Lorentz factor $\sqrt{1-v^{2}}$ except when the
ratio of their respective widths have a particular value.

In this work we do not consider corrections due to processes taking place over
a finite-time interval (instead of an infinite time interval). This has been
well studied in other work [4, 5, 6, 7]. More recently, in [8] a systematic
study was undertaken to derive Fermi’s golden rule in QFT within a Gaussian
wave packet formalism [9, 10, 11, 12, 13]. They succeeded in separating bulk
from boundary time contributions while maintaining the unitarity of the
S-matrix in contrast to earlier work [14]. They showed that the bulk
contribution leads to Fermi’s golden rule while the boundary part can
contribute deviations from it. This provided a solid framework in which to
view the results of the aforementioned work[4, 5, 6, 7].

In [15] a wave packet analysis of quantum correlations in neutrino
oscillations[16, 17, 18] was undertaken to extend previous work [19] that had
been carried out using a plane wave analysis. The claim is that this wave
packet analysis leads to a better agreement with the experimental results,
most notably the MINOS experiments [20, 21] (the corrections due to the wave
packet are more negligible for the Daya Bay experiments[22, 23]). In [24]
neutrino wave packet effects for medium-baseline reactor neutrino oscillation
(MBRO) experiments were also studied. Though it needs to be confirmed, effects
of wave packets may already be showing up in experiments, in particular those
involving neutrino oscillations.

Our paper is organized as follows. We first determine the decay rate of a
Gaussian wave packet centered at zero momentum and determine the correction to
the plane wave result due to the width $a$ of the Gaussian. We then determine
how a wave packet transforms under a Lorentz boost and calculate the decay
rate of the Lorentz-boosted wave packet. This is where the WP clock is
introduced. The Lorentz-boosted wave packet is not a Gaussian and we therefore
calculate separately the decay rate of a Gaussian wave packet centered at
velocity $v$ (i.e. momentum $M\,\gamma\,v$). The conclusion discusses an
important and interesting problem that emerges from this work. The integrals
we encounter in this paper do not have an exact analytical solution but can be
evaluated via a power series expansion up to arbitrary accuracy. The details
are relegated to Appendix A.

##  2 Decay rate of a wave packet

The state $\ket{\phi}$ of a wave packet in quantum field theory can be
expressed as a linear superposition of one-particle states [1]

$\ket{\phi}=\int\dfrac{d^{3}k}{(2\pi)^{3}}\dfrac{1}{\sqrt{2E_{\bf
k}}}\phi({\bf k})\ket{\bf k}$ (1)

where $\ket{\bf k}$ is a one-particle state of momentum ${\bf k}$, $E_{\bf k}$
is the relativistic energy $\sqrt{{\bf k}^{2}+m^{2}}$ and $\phi({\bf k})$ is
the Fourier transform of the spatial wavefunction. We use the Lorentz
invariant normalization [1, 2] $\braket{{\bf p}}{{\bf q}}=2E_{\bf
p}(2\pi)^{3}\delta^{3}({\bf p}-{\bf q})$. We then obtain the conventional
normalization $\braket{\phi}{\phi}=1$ (where the sum of all probabilities is
unity) if

$\int\dfrac{d^{3}k}{(2\pi)^{3}}|\phi({\bf k})|^{2}=1\,.$ (2)

The differential decay rate of an unstable particle of mass $M$ treated as a
plane wave (subscript “0”) with four-momentum $k$ (where $k^{0}=E_{\bf k}$)
decaying into a set of final particles with four-momenta $p_{f}$ is given by
[1]

$\displaystyle d\Gamma_{0}$ $\displaystyle=\dfrac{1}{2E_{\bf
k}}\Big{(}\prod_{f}\dfrac{d^{3}{\bf
p}_{f}}{(2\pi)^{3}}\dfrac{1}{2E_{f}}\Big{)}|\mathcal{M}(k\to\\{p_{f}\\})|^{2}(2\pi)^{4}\delta^{4}(k-\textstyle{\sum}\,p_{f})$
$\displaystyle=\dfrac{1}{2E_{\bf k}}\,d\xi$ (3)

where

$d\xi=\Big{(}\prod_{f}\dfrac{d^{3}{\bf
p}_{f}}{(2\pi)^{3}}\dfrac{1}{2E_{f}}\Big{)}|\mathcal{M}(k\to\\{p_{f}\\})|^{2}(2\pi)^{4}\delta^{4}(k-\textstyle{\sum}\,p_{f})$
(4)

is Lorentz invariant and hence independent of ${\bf k}$ (e.g. if you evaluate
it say at ${\bf k}=0$ then the result is valid for all ${\bf k}$). The
dependence of the decay rate on ${\bf k}$ appears only in the prefactor
$\frac{1}{2E_{\bf k}}$. This factor is not Lorentz invariant and is the part
responsible for time dilation i.e. that the lifetime for a particle increases
with its speed. If we consider the unstable particle as a wave packet, we only
need to integrate over this factor. The differential decay rate for the wave
packet is then given by

$d\Gamma_{\text{packet}}=\Big{(}\int\dfrac{d^{3}k}{(2\pi)^{3}}|\phi({\bf
k})|^{2}\,\dfrac{1}{2E_{\bf k}}\Big{)}\,d\xi\,.$ (5)

###  2.1 Decay rate of Gaussian wave packet centered at zero momentum: width
correction

For a particle of mass $M$ at rest (${\bf k}=0$), we obtain from (3) that

$d\xi=2\,Md\Gamma_{0_{({\bf k}=0)}}$ (6)

where $d\Gamma_{0_{({\bf k}=0)}}$ is the plane wave result for the
differential decay rate of a particle at rest. We will now represent a
“particle at rest” by a Gaussian wavepacket centered at ${\bf k}=0$ with width
$a$ and normalized according to (2). This is given by

$\phi({\bf k})=\Big{(}\dfrac{2}{a}\Big{)}^{3/2}\pi^{3/4}e^{-{\bf
k}^{2}/(2\,a^{2})}\,.$ (7)

Substituting this into (5) and using (6) we obtain that

$\displaystyle d\Gamma_{\text{packet}_{({\bf k}=0)}}$
$\displaystyle=\Big{(}\int_{0}^{\infty}\dfrac{2}{a^{3}\sqrt{\pi}}\,{\bf
k}^{2}e^{-{\bf k}^{2}/a^{2}}\dfrac{1}{\sqrt{{\bf k}^{2}+M^{2}}}d|{\bf
k}|\Big{)}2M\,d\Gamma_{0_{({\bf k}=0)}}$
$\displaystyle=\dfrac{1}{a}\,\text{U}(1/2,0,M^{2}/a^{2})M\,d\Gamma_{0_{({\bf
k}=0)}}$ (8)

where $\text{U}(a,b,z)$ is known as the confluent hypergeometric function of
the second kind[3]. We are interested in obtaining the correction to the plane
wave ($a=0$) result. We therefore perform a series expansion of $\text{U}/a$
about $a=0$ and this yields

$\dfrac{1}{a}\,\text{U}(1/2,0,M^{2}/a^{2})=\frac{1}{M}\Big{(}1-\frac{3\,a^{2}}{4M^{2}}+\mathcal{O}(\tfrac{a^{4}}{M^{4}})\Big{)}\,.$
(9)

Substituting the above into (8) we obtain our first result:

$\Gamma_{\text{packet}_{({\bf k}=0)}}=\Gamma_{0_{({\bf
k}=0)}}\,\Big{(}1-\frac{3a^{2}}{4M^{2}}+\mathcal{O}(\tfrac{a^{4}}{M^{4}})\Big{)}$
(10)

where $\Gamma_{0_{({\bf k}=0)}}$ is the plane wave result for the decay rate
of a particle at rest and $\Gamma_{\text{packet}_{({\bf k}=0)}}$ is the decay
rate of a Gaussian wave packet centered at zero momentum (note that we dropped
the differential symbol because the relation (10) holds for the total decay
rates). We see that the first correction due to the wave packet is of order
$\tfrac{a^{2}}{M^{2}}$ and vanishes in the plane wave limit where $a$ tends to
zero222We work in units where $\hbar=c=1$. Reinstating these two constants the
correction is of order
$\tfrac{\hbar^{2}a^{2}}{M^{2}c^{2}}=\lambda_{c}^{2}\,a^{2}$ where
$\lambda_{c}=\tfrac{\hbar}{Mc}$ is the Compton wavelength of the particle..
Higher-order corrections such as $\tfrac{a^{4}}{M^{4}}$, etc. are expected to
be much smaller and are not included here. The negative sign in (10) makes
physical sense because we expect the lifetime to increase due to the time
dilation produced by non-zero momenta in the wave packet (since there is a
spread of momenta about zero).

###  2.2 The Lorentz boost of a wave packet

We now determine how $\phi({\bf k})$ transforms under a Lorentz boost. The
state of a wavepacket $\ket{\phi}$ is given by (1). After a Lorentz
transformation, we label all quantities by a prime i.e.

$\ket{\phi^{\prime}}=\int\dfrac{d^{3}k^{\prime}}{(2\pi)^{3}}\dfrac{1}{\sqrt{2E_{\bf
k}^{\prime}}}\phi^{\prime}({\bf k}^{\prime})\ket{\bf k^{\prime}}\,.$ (11)

Quantum states transform under a Lorentz transformation $\Lambda$ via a
unitary operator $U(\Lambda)$ [1]. For example, $\ket{\bf k^{\prime}}$ is
given by $U(\Lambda)\ket{\bf k}$. Therefore

$\displaystyle\ket{\phi^{\prime}}$
$\displaystyle=U(\Lambda)\ket{\phi}=\int\dfrac{d^{3}k}{(2\pi)^{3}}\dfrac{1}{\sqrt{2E_{\bf
k}}}\phi({\bf k})\,\,(U(\Lambda)\ket{\bf k})$
$\displaystyle=\int\dfrac{d^{3}k}{(2\pi)^{3}}\dfrac{1}{\sqrt{2E_{\bf
k}}}\phi({\bf k})\,\ket{\bf k^{\prime}}\,.$ (12)

Comparing (12) with (11) we find that

$\dfrac{d^{3}k}{\sqrt{E_{\bf k}}}\,\phi({\bf
k})=\dfrac{d^{3}k^{\prime}}{\sqrt{E_{\bf k}^{\prime}}}\,\phi^{\prime}({\bf
k}^{\prime})\,.$ (13)

In other words, the (continuous) coefficients in the expansion of a wave
packet are Lorentz invariant. Under a Lorentz boost in the 3-direction we have
$k_{3}^{\prime}=\gamma(k_{3}+v\,E)$ and $E^{\prime}=\gamma E(1+v\,k_{3}/E)$ so
that $\frac{d^{3}k^{\prime}}{\sqrt{E_{\bf
k}^{\prime}}}=\frac{d^{3}k}{\sqrt{E_{\bf k}}}[\gamma\,(1+v\,k_{3}/E)]^{1/2}$.
Substituting this into (13) we obtain

$\phi^{\prime}({\bf k}^{\prime})=\phi({\bf
k})\,[\gamma\,(1+v\,k_{3}/E)]^{-1/2}\,.$ (14)

The above describes how $\phi({\bf k})$ transforms under a Lorentz boost. The
wave packet remains normalized, as it should, after the Lorentz transformation
i.e.

$\int\dfrac{d^{3}k^{\prime}}{(2\pi)^{3}}\,|\phi^{\prime}({\bf
k}^{\prime})|^{2}=\int\dfrac{d^{3}k}{(2\pi)^{3}}\,[\gamma\,(1+v\,k_{3}/E)]\,\,\dfrac{|\phi({\bf
k})|^{2}}{[\gamma\,(1+v\,k_{3}/E)]}=\int\dfrac{d^{3}k}{(2\pi)^{3}}\,|\phi({\bf
k})|^{2}=1\,.$ (15)

We now determine $|\phi^{\prime}({\bf k}^{\prime})|^{2}$ in the primed
coordinates when $\phi({\bf k})$ is given by the Gaussian (7). The Lorentz
boost in the 3-direction yields the following relations:
$k_{3}=\gamma\,(k_{3}^{\prime}-vE^{\prime})$ and
$\gamma\,(1+v\,k_{3}/E)=[\gamma\,(1-v\,k_{3}^{\prime}/E^{\prime})]^{-1}$. We
then obtain

$|\phi^{\prime}({\bf
k}^{\prime})|^{2}=\dfrac{8\pi^{3/2}}{a^{3}}\,[\gamma\,(1-v\,k_{3}^{\prime}/E^{\prime})]\,e^{-(k_{1}^{2}+k_{2}^{2}+[\gamma\,(k_{3}^{\prime}-vE^{\prime})]^{2})/a^{2}}\,.$
(16)

This is clearly not a Gaussian. We therefore see that a Gaussian does not
transform into another Gaussian under a Lorentz boost. Though the exponential
peaks at $k_{1}=0,k_{2}=0,k_{3}^{\prime}=M\,\gamma\,v$, (16) does not
represent a Gaussian wave packet centered at momentum $M\,\gamma\,v$. We will
study such a Gaussian in section 2.3. We now turn to determining the decay
rate of the Lorentz-boosted wavepacket.

####  2.2.1 Decay rate of a Lorentz-boosted wave packet and the WP clock

The differential decay rate of a wave packet that is Lorentz boosted (assumed
to be in the 3- direction) is given by (5) with primes added to represent the
Lorentz boosted quantities:

$d\Gamma_{\text{packet}_{boosted}}=\Big{(}\int\dfrac{d^{3}k^{\prime}}{(2\pi)^{3}}|\phi^{\prime}({\bf
k^{\prime}})|^{2}\,\dfrac{1}{2E_{\bf k}^{\prime}}\Big{)}\,d\xi\,.$ (17)

where $d\xi$ is given by the Lorentz-invariant quantity (4). Before boosting,
we start with the Gaussian wavepacket centered at zero momentum given by (7)

$|\phi({\bf
k})|^{2}=\dfrac{8\pi^{3/2}}{a^{3}}\,e^{-(k^{2}+k_{3}^{2})/a^{2}}\,.$ (18)

where $k^{2}=k_{1}^{2}+k_{2}^{2}$ so that $d^{3}k=2\pi\,kdk\,dk_{3}$. From the
calculation carried out in (15) we can see that
$d^{3}k^{\prime}\,|\phi^{\prime}({\bf k^{\prime}})|^{2}=d^{3}k\,|\phi({\bf
k})|^{2}$. Moreover, $E_{\bf k}^{\prime}=\sqrt{k^{2}+k_{3}^{\prime 2}+M^{2}}$
where $k_{3}^{\prime}=\gamma\,(k_{3}+v\,E)$ with
$E=\sqrt{k^{2}+k_{3}^{2}+M^{2}}$. The decay rate of the Lorentz boosted wave
packet is then given by

$d\Gamma_{\text{packet}_{boosted}}=\Big{(}\dfrac{2}{a^{3}\,\sqrt{\pi}}\int_{0}^{\infty}\int_{-\infty}^{\infty}\,dk\,dk_{3}\,\dfrac{k\,e^{-(k^{2}+k_{3}^{2})/a^{2}}}{2\sqrt{k^{2}+[\gamma\,(k_{3}+v\,E)]^{2}+M^{2}}}\Big{)}\,d\xi\,.$
(19)

The integrations converge but do not yield an analytical result in terms of
any well-known function. However, since the exponential peaks at $k=k_{3}=0$
we can perform a series expansion about $k=0$ and $k_{3}=0$ of the function
that multiplies the exponential and then integrate. This leads to an
analytical result as a series in powers of the small parameter $a^{2}/M^{2}$
up to arbitrary accuracy. The result is (see Appendix A)

$\displaystyle\Gamma_{\text{packet}_{boosted}}$
$\displaystyle=\sqrt{1-v^{2}}\,\dfrac{1}{2M}\Big{[}1-\dfrac{3a^{2}}{4M^{2}}\Big{(}1-\dfrac{2\,v^{2}}{3}\Big{)}+\mathcal{O}(\tfrac{a^{4}}{M^{4}})\Big{]}\,d\xi$
$\displaystyle=\sqrt{1-v^{2}}\,\,\Gamma_{0_{({\bf
k}=0)}}\,\Big{[}1-\dfrac{3a^{2}}{4M^{2}}\Big{(}1-\dfrac{2\,v^{2}}{3}\Big{)}+\mathcal{O}(\tfrac{a^{4}}{M^{4}})\Big{]}$
$\displaystyle=\sqrt{1-v^{2}}\,\,\Gamma_{\text{packet}_{({\bf
k}=0)}}\Bigg{[}\dfrac{1-\dfrac{3a^{2}}{4M^{2}}\Big{(}1-\dfrac{2\,v^{2}}{3}\Big{)}+\mathcal{O}(\tfrac{a^{4}}{M^{4}})}{1-\dfrac{3a^{2}}{4M^{2}}+\mathcal{O}(\tfrac{a^{4}}{M^{4}})}\Bigg{]}$
$\displaystyle=\sqrt{1-v^{2}}\,\Gamma_{\text{packet}_{({\bf
k}=0)}}\,\Big{(}1+\dfrac{a^{2}\,v^{2}}{2M^{2}}+\mathcal{O}(\tfrac{a^{4}}{M^{4}})\Big{)}$
(20)

where we used (6) for $d\xi$ and used (10) to express $\Gamma_{0_{({\bf
k}=0)}}$ in terms of $\Gamma_{\text{packet}_{({\bf k}=0)}}$. In the last
equality we performed a binomial expansion and kept terms up to order
$a^{2}/M^{2}$.

The important point to take away from the above result is that the decay rate
of a wave packet after a Lorentz boost is not related to the original decay
rate of the wave packet by exactly the Lorentz factor $\sqrt{1-v^{2}}$. There
is a small correction333Note that compared to the decay rate of a boosted
plane wave the correction is equal to
$-\tfrac{3a^{2}}{4M^{2}}(1-\tfrac{2\,v^{2}}{3})\,.$ equal to
$\tfrac{a^{2}\,v^{2}}{2\,M^{2}}\,$. The correction increases as the boost
velocity increases (reaching a maximum in the limit $v\to 1$).

We therefore see that the decay of a localized particle in QFT does not act
exactly like a typical clock under Lorentz boosts and we refer to it as a “WP
clock” (wave packet clock). In contrast to a typical clock, a WP clock does
not move with a single velocity but has a spread in velocities (or more
specifically, a spread in momenta). It is best to view the WP clock as a
single clock since the wave packet represents a one-particle state in QFT. It
should be emphasized that there is no breakdown of Lorentz symmetry here. The
effect is due to the fact that localized particles in QFT must be represented
by wave packets and not plane waves in accordance with Heisenberg’s
uncertainty principle.

###  2.3 Decay rate of Gaussian wave packet centered at velocity $v$

As we saw, the Lorentz boost by a velocity $v$ of a Gaussian wave packet
centered at zero velocity does not yield a Gaussian wave packet centered at
velocity $v$. It is therefore of interest to determine the decay rate of a
Gaussian wavepacket centered at velocity $v$ and compare it to the decay rate
of a Gaussian at “rest”. The third component $k_{3}$ of the momentum will be
centered at $M\gamma v$ where $\gamma$ is the usual Lorentz factor
$(1-v^{2})^{-1/2}$. The other two components, $k_{1}$ and $k_{2}$, will both
be centered at zero. The width is denoted by $b$ and the Gaussian normalized
again according to (2). This yields

$\phi({\bf k})=\Big{(}\dfrac{2}{b}\Big{)}^{3/2}\pi^{3/4}e^{-[{\bf
p}^{2}+(k_{3}-M\,\gamma v)^{2}]/(2\,b^{2})}$ (21)

where ${\bf p}=(k_{1},k_{2})$ and ${\bf p}^{2}=k_{1}^{2}+k_{2}^{2}$.
Substituting the above into (5) and using
$d^{3}k=dk_{1}\,dk_{2}\,dk_{3}=2\pi|{\bf p}|d|{\bf p}|dk_{3}$ we obtain

$\displaystyle d\Gamma_{\text{packet}_{(|{\bf k}|=M\gamma v)}}$
$\displaystyle=d\xi\,\int_{0}^{\infty}\int_{-\infty}^{\infty}\dfrac{1}{b^{3}\sqrt{\pi}}\,|{\bf
p}|\dfrac{e^{-[{\bf p}^{2}+(k_{3}-M\,\gamma v)^{2}]/b^{2}}}{\sqrt{{\bf
p}^{2}+k_{3}^{2}+M^{2}}}d|{\bf p}|\,dk_{3}\,.$ (22)

The above integral converges but does not yield an analytical result in terms
of any well-known function. However, as in the previous section, we can obtain
an analytical result as a power series expansion in the small parameter
$b^{2}/M^{2}$ up to arbitrary accuracy. The exponential in this case peaks at
$|{\bf p}|=0$ and $k_{3}=M\,\gamma v$. We can therefore perform a series
expansion about $|{\bf p}|=0$ and $k_{3}=M\,\gamma v$ of the function that
multiplies the exponential and then integrate. The result is (see Appendix A
for details):

$\displaystyle d\Gamma_{\text{packet}_{(|{\bf k}|=M\gamma
v)}}=d\xi\,\sqrt{1-v^{2}}\,\,\Big{[}\dfrac{1}{2M}-\dfrac{3b^{2}}{8M^{3}}(1-v^{2})^{2}+\mathcal{O}(\tfrac{b^{4}}{M^{4}})\Big{]}$
(23)

Substituting $d\xi$ from (6) into (23) we obtain

$\displaystyle\Gamma_{\text{packet}_{(|{\bf k}|=M\gamma
v)}}=\sqrt{1-v^{2}}\,\,\Gamma_{0_{({\bf
k}=0)}}\Big{[}1-\dfrac{3b^{2}}{4M^{2}}(1-v^{2})^{2}+\mathcal{O}(\tfrac{b^{4}}{M^{4}})\Big{]}$
(24)

where we dropped the differential symbol as before (since the relation applies
to the total decay rate). Substituting $\Gamma_{0_{({\bf k}=0)}}$ from (10)
into (24) we obtain

$\displaystyle\Gamma_{\text{packet}_{(|{\bf k}|=M\gamma v)}}$
$\displaystyle=\sqrt{1-v^{2}}\,\,\Gamma_{\text{packet}_{({\bf
k}=0)}}\,\Bigg{[}\dfrac{1-\frac{3b^{2}}{4M^{2}}(1-v^{2})^{2}+\mathcal{O}(\tfrac{b^{4}}{M^{4}})}{1-\frac{3a^{2}}{4M^{2}}+\mathcal{O}(\tfrac{a^{4}}{M^{4}})}\Bigg{]}$
$\displaystyle=\sqrt{1-v^{2}}\,\Gamma_{\text{packet}_{({\bf
k}=0)}}\,\Big{[}1+\frac{3}{4M^{2}}\Big{(}a^{2}-b^{2}\,(1-v^{2})^{2}\Big{)}+\mathcal{O}(\tfrac{b^{4}}{M^{4}},\tfrac{a^{4}}{M^{4}},\tfrac{a^{2}\,b^{2}}{M^{4}})\Big{]}\,.$
(25)

We therefore see that $\Gamma_{\text{packet}_{(|{\bf k}|=M\,\gamma\,v)}}$ is
not exactly equal to $\sqrt{1-v^{2}}\,\,\Gamma_{\text{packet}_{({\bf k}=0)}}$.
In other words, the two decay rates are not related exactly by the Lorentz
factor $\sqrt{1-v^{2}}$. For two plane waves (the $a\to 0$ and $b\to 0$ limit)
the relation is exact as expected. For the case where the two widths are equal
($a=b$), the correction to the Lorentz factor is non-zero and equal to
$\tfrac{3\,a^{2}\,v^{2}}{4\,M^{2}}(2-v^{2})$. The correction increases with
$v$ and reaches a maximum in the limit $v\to 1$.

In particular, note that if the width $b$ of the second (moving) Gaussian is
equal to $a\,\gamma^{2}$ where $a$ is the width of the Gaussian at “rest”,
then the leading order correction vanishes and the two decay rates are related
by the Lorentz factor (up to higher order corrections). So the decay rate of a
Gaussian centered at velocity $v$ and one centered at zero velocity are not in
general related exactly by the Lorentz factor except when the ratio of their
widths $b/a$ is equal to $\gamma^{2}$.

##  3 Conclusion

In this paper, we found an analytical result for the leading order width
correction to the decay rate of a Gaussian wave packet centered at zero
momentum. The effect was of order $\tfrac{a^{2}}{M^{2}}$ where $a$ is the
width and $M$ is the mass of the decaying particle. We then showed that the
decay rate of this Gaussian wave packet under a Lorentz boost by velocity $v$
does not decrease exactly by the Lorentz factor $\sqrt{1-v^{2}}$. There is a
correction of order $\tfrac{a^{2}v^{2}}{M^{2}}$. This is important because it
means that a decaying wave packet does not behave like a typical clock under
Lorentz boosts and we refer it to it as a “WP clock”. A WP clock can be
defined as a clock that does not move with a single velocity relative to an
observer but has a spread in velocities (or a spread in momenta). Since the
wave packet is a one-particle state in QFT it is best to view the WP clock as
a single clock. As already mentioned, the WP clock does not imply a breakdown
of Lorentz symmetry: it does not stem from new physics but is an important
consequence of the combined requirements of special relativity, quantum
mechanics and localized free particles.

We found that a Gaussian wave packet which is Lorentz-boosted does not
transform into another Gaussian. This is why we calculate separately the decay
rate of a Gaussian centered at momentum $M\,\gamma\,v$ (with width $b$). We
find that it is not exactly equal to the decay rate of a Gaussian centered at
zero momentum (with width $a$) multiplied by the Lorentz factor
$\sqrt{1-v^{2}}$. There is again a correction that depends on their respective
widths and the velocity $v$. For the case where their respective widths are
equal ($a=b$), the leading order correction to the Lorentz factor is equal to
$\tfrac{3\,a^{2}\,v^{2}}{4\,M^{2}}(2-v^{2})$. However, the leading order
correction vanishes if the ratio of their widths $b/a$ is equal to
$\gamma^{2}$.

There is one thing in particular that would be important to investigate in
future work. In special relativity, there is a symmetry between two inertial
frames. If two frames A and B move with constant velocity relative to each
other, then from A’s perspective, B’s clock runs slower by the usual time
dilation factor and from B’s perspective, A’s clock runs slower by the same
time dilation factor. Consider now a Gaussian wave packet that is centered at
zero momentum relative to an inertial frame A. Let its decay rate in A’s frame
be $\Gamma_{A}$. If we now Lorentz boost the wave packet with velocity $v$ the
decay rate of the Lorentz-boosted wave packet will be $\Gamma_{A}$ multiplied
by the usual Lorentz factor of $\sqrt{1-v^{2}}$ but with a correction due to
the width and speed $v$. In other words, from A’s perspective, the clock now
runs slower but with a small correction to the usual time dilation: this
correction is what defines a WP clock. Now, frame A, by definition, is an
inertial frame with a typical clock. By symmetry, its clock relative to a
frame B attached to the free particle wave packet should also run slower by
the exact same amount as the WP clock. Since A’s clock is a typical clock,
this implies that the correction that appears in a WP clock must now originate
from the nature of frame B. In other words, if we want to preserve the
symmetry inherent to special relativity, frame B attached to the free particle
wave packet must be a slightly different kind of inertial frame (just like a
WP clock is a slightly different kind of clock). One argument in favor of
regarding frame B as an inertial frame is that it is a frame attached to a
free particle. However, in contrast to a free particle in classical mechanics,
a free localized particle in QFT does not have a precise velocity/momentum. So
if frame B can be regarded as an inertial frame, it makes sense that it would
be a slightly different kind of inertial frame. The goal would then be to see
if one can define such an inertial frame so that it is consistent with special
relativity and leads to the desired effects. This is a worthy goal to pursue
with the potential of broadening our scope.

##  Appendix A Integrations via series expansions

###  A.1 Integrals appearing in the decay rate of a boosted wave packet

In section 2.2.1 we encountered the decay rate of a boosted wave packet given
by the double integral (19):

$\displaystyle
d\Gamma_{\text{packet}_{boosted}}=\Big{(}\dfrac{2}{a^{3}\,\sqrt{\pi}}\int_{0}^{\infty}\int_{-\infty}^{\infty}\,dk\,dk_{3}\,\dfrac{k\,e^{-(k^{2}+k_{3}^{2})/a^{2}}}{2\sqrt{k^{2}+[\gamma\,(k_{3}+v\,E)]^{2}+M^{2}}}\Big{)}\,d\xi\,$
$\displaystyle=\Big{(}\dfrac{1}{2\,a^{3}\,\sqrt{\pi}}\int_{0}^{\infty}\int_{-\infty}^{\infty}\,du\,dk_{3}\,\dfrac{e^{-(u+k_{3}^{2})/a^{2}}}{\sqrt{u+[(1-v^{2})^{-1/2}\,(k_{3}+v\sqrt{u+k_{3}^{2}+M^{2}})]^{2}+M^{2}}}\Big{)}\,d\xi$
(A.1)

where we made the substitution $u=k^{2}$, replaced $E$ by its expression
$\sqrt{k^{2}+k_{3}^{2}+M^{2}}=\sqrt{u+k_{3}^{2}+M^{2}}$ and replaced $\gamma$
by $(1-v^{2})^{-1/2}$. The integrals converge but do not yield an analytical
result in terms of any well-known function. However, by expanding the function
that multiplies the exponential about $u=k_{3}=0$ (where the exponential
peaks) as a power series before integrating we can obtain an analytical result
for the integral as a power series in the small parameter $a^{2}/M^{2}$ up to
arbitrary accuracy. The function is (including the normalization pre-factor)

$f(u,k_{3})=\dfrac{1}{2\,a^{3}\,\sqrt{\pi}}\Big{[}u+[(1-v^{2})^{-1/2}\,(k_{3}+v\sqrt{u+k_{3}^{2}+M^{2}})]^{2}+M^{2}\Big{]}^{-1/2}\,.$
(A.2)

Our goal is to obtain the leading correction to the plane wave result for the
decay rate. This means we want the integral (with the normalization pre-
factor) to be accurate up to $a^{2}/M^{2}$ (but do not require it to be
accurate up to $a^{4}/M^{4}$ or higher). This requires that we expand the
function $f$ above only up to second order in $k_{3}$ and first order in $u$.
This yields the function $f_{2}$ given by

$\displaystyle
f_{2}=\dfrac{\sqrt{1-v^{2}}}{\sqrt{\pi}\,a^{3}}\,\dfrac{1}{2\,M}\Big{[}1-\dfrac{u}{2\,M^{2}}+$
$\displaystyle-\dfrac{k_{3}}{M^{3}}(M^{2}-u)v$
$\displaystyle\quad\quad+\dfrac{k_{3}^{2}}{4\,M^{4}}(2M^{2}-3\,u)(-1+2v^{2})\Big{]}\,.$
(A.3)

Multiplying $f_{2}$ by the exponential and then performing the integration
yields

$\int_{0}^{\infty}\int_{-\infty}^{\infty}\,du\,dk_{3}\,f_{2}\,e^{-(u+k_{3}^{2})/a^{2}}=\sqrt{1-v^{2}}\,\dfrac{1}{2M}\Big{(}1-\dfrac{3a^{2}}{4M^{2}}(1-2\,v^{2}/3)+\mathcal{O}(\tfrac{a^{4}}{M^{4}})\Big{)}\,.$
(A.4)

The differential decay rate of a boosted wave packet given by (A.1) is
obtained by multiplying the above result by the factor
$d\xi=2\,Md\Gamma_{0_{({\bf k}=0)}}$ given by (6) where $d\Gamma_{0_{({\bf
k}=0)}}$ is the differential decay rate of a particle at rest treated as a
plane wave. This yields

$\Gamma_{\text{packet}_{boosted}}=\sqrt{1-v^{2}}\,\Gamma_{0_{({\bf
k}=0)}}\Big{(}1-\dfrac{3a^{2}}{4M^{2}}(1-2\,v^{2}/3)+\mathcal{O}(\tfrac{a^{4}}{M^{4}})\Big{)}$
(A.5)

where we have dropped the differential sign since the relation applies also to
the total decay rate. We would like to express the decay rate of the boosted
wave packet in terms of the decay rate of the original wave packet “at rest”
(centered at zero momentum). Using (10) we can express the decay rate of a
plane wave at rest to the decay rate of a wave packet “at rest” i.e.

$\Gamma_{0_{({\bf k}=0)}}=\dfrac{\Gamma_{\text{packet}_{({\bf
k}=0)}}}{1-\frac{3a^{2}}{4M^{2}}+\mathcal{O}(\tfrac{a^{4}}{M^{4}})}\,.$ (A.6)

Substituting the above into (A.5) we obtain

$\displaystyle\Gamma_{\text{packet}_{boosted}}$
$\displaystyle=\sqrt{1-v^{2}}\,\,\Gamma_{\text{packet}_{({\bf
k}=0)}}\Bigg{[}\dfrac{1-\dfrac{3a^{2}}{4M^{2}}(1-2\,v^{2}/3)+\mathcal{O}(\tfrac{a^{4}}{M^{4}})}{1-\dfrac{3a^{2}}{4M^{2}}+\mathcal{O}(\tfrac{a^{4}}{M^{4}})}\Bigg{]}$
$\displaystyle=\sqrt{1-v^{2}}\,\Gamma_{\text{packet}_{({\bf
k}=0)}}\,\Big{(}1+\dfrac{a^{2}\,v^{2}}{2M^{2}}+\mathcal{O}(\tfrac{a^{4}}{M^{4}})\Big{)}\,$
(A.7)

where in the last line we performed a binomial expansion and kept terms only
up to the small parameter $a^{2}/M^{2}$. We therefore see that the Lorentz
boost of a wave packet does not change its decay rate exactly by the Lorentz
factor $\sqrt{1-v^{2}}$. There is a correction of order $a^{2}\,v^{2}/M^{2}$.

###  A.2 Integrals appearing in the decay rate of a Gaussian wave packet
centered at velocity v

In section 2.3 we encountered the decay rate of a Gaussian wave packet
centered at velocity v given by the double integral (22)

$\displaystyle d\Gamma_{\text{packet}_{(|{\bf k}|=M\gamma v)}}$
$\displaystyle=d\xi\,\int_{0}^{\infty}\int_{-\infty}^{\infty}\dfrac{1}{b^{3}\sqrt{\pi}}\,|{\bf
p}|\dfrac{e^{-[{\bf p}^{2}+(k_{3}-M\,\gamma v)^{2}]/b^{2}}}{\sqrt{{\bf
p}^{2}+k_{3}^{2}+M^{2}}}d|{\bf p}|\,dk_{3}$
$\displaystyle=d\xi\,\int_{0}^{\infty}\int_{-\infty}^{\infty}\dfrac{1}{2b^{3}\sqrt{\pi}}\dfrac{e^{-[u+p_{3}^{2}]/b^{2}}}{\sqrt{u+(p_{3}+M\,\gamma
v)^{2}+M^{2}}}du\,dp_{3}$ (A.8)

where we made the substitution $u=|{\bf p}|^{2}$ and $p_{3}=k_{3}-M\,\gamma
v$. The integrals converge but do not yield an analytical result in terms of
any well-known function. Our procedure will be the same as the one used for
the previous integral in the subsection above. The exponential peaks at
$u=p_{3}=0$. We therefore expand the function that multiplies the exponential
as a power series about $u=0$ and $p_{3}=0$. We can then obtain an analytical
result for the integral as a power series in the small parameter $b^{2}/M^{2}$
up to arbitrary accuracy. The function in this case is

$f=\dfrac{1}{2b^{3}\sqrt{\pi}}\Big{[}u+(p_{3}+M\,\gamma
v)^{2}+M^{2}\Big{]}^{-1/2}\,.$ (A.9)

We are interested in obtaining only the leading correction to the decay rate
proportional to $b^{2}/M^{2}$. This requires us to expand $f$ only up to
second order in $p_{3}$ and first order in $u$. We label this $f_{3}$ and it
is given by

$\displaystyle f_{3}=\dfrac{1}{8\,b^{3}\,M^{5}\,\sqrt{\pi}}$
$\displaystyle\Big{(}4\,M^{4}\,\sqrt{1-v^{2}}-4\,M^{3}\,p_{3}\,v\,(1-v^{2})+6\,M\,p_{3}\,u\,v(1-v^{2})^{2}$
$\displaystyle+3\,p_{3}^{2}\,u(1-v^{2})^{5/2}(1-5v^{2})-2\,M^{2}(1-v^{2})^{3/2}(u+p_{3}^{2}(1-3v^{2}))\Big{)}\,.$
(A.10)

The decay rate (A.8) is then given by

$\displaystyle\Gamma_{\text{packet}_{(|{\bf k}|=M\gamma v)}}$
$\displaystyle=d\xi\,\int_{0}^{\infty}\int_{-\infty}^{\infty}\,du\,dp_{3}\,f_{3}\,e^{-(u+p_{3}^{2})/b^{2}}$
$\displaystyle=\sqrt{1-v^{2}}\,\dfrac{d\xi}{2M}\Big{(}1-\dfrac{3b^{2}}{4M^{2}}(1-v^{2})^{2}+\mathcal{O}(\tfrac{b^{4}}{M^{4}})\Big{)}$
(A.11)

where we dropped the differential sign for the decay rate as the result holds
for the total decay rate. Using (6) and (10) we now replace $d\xi$ by

$d\xi=2M\,\Gamma_{0_{({\bf k}=0)}}=2M\,\dfrac{\Gamma_{\text{packet}_{({\bf
k}=0)}}}{1-\frac{3a^{2}}{4M^{2}}+\mathcal{O}(\tfrac{a^{4}}{M^{4}})}\,.$ (A.12)

Substituting the above into the decay rate (A.11) we obtain our final result

$\displaystyle\Gamma_{\text{packet}_{(|{\bf k}|=M\gamma v)}}$
$\displaystyle=\sqrt{1-v^{2}}\,\,\Gamma_{\text{packet}_{({\bf
k}=0)}}\,\Bigg{[}\dfrac{1-\frac{3b^{2}}{4M^{2}}(1-v^{2})^{2}+\mathcal{O}(\tfrac{b^{4}}{M^{4}})}{1-\frac{3a^{2}}{4M^{2}}+\mathcal{O}(\tfrac{a^{4}}{M^{4}})}\Bigg{]}$
$\displaystyle=\sqrt{1-v^{2}}\,\Gamma_{\text{packet}_{({\bf
k}=0)}}\,\Big{[}1+\frac{3}{4M^{2}}\Big{(}a^{2}-b^{2}\,(1-v^{2})^{2}\Big{)}+\mathcal{O}(\tfrac{b^{4}}{M^{4}},\tfrac{a^{4}}{M^{4}})\Big{]}\,.$
(A.13)

## Acknowledgments

The author acknowledges support from a discovery grant of the National Science
and Engineering Research Council of Canada (NSERC). The author also thanks
Michael E. Peskin for valuable communications.

## References

  * [1] M. E. Peskin and D. V. Schroeder, An Introduction to Quantum Field Theory, (Westview Press, USA, 1995).
  * [2] C. Burgess and G. Moore, The Standard Model: A Primer,(Cambridge University Press, UK, 2007).
  * [3] A. B. Olde Daalhuis, Confluent Hypergeometric Functions, NIST Handbook of Mathematical Functions(Cambridge University Press, Cambridge, UK, 2010).
  * [4] K. Ishikawa and Y. Tobita, Matter-enhanced transition probabilities in quantum field theory, Annals Phys. 344 (2014) 118 [arXiv:1206.2593].
  * [5] K. Ishikawa and Y. Tobita, Finite-size corrections to Fermi’s golden rule: I. Decay rates, PTEP 2013 (2013) 073B02 [arXiv:1303.4568].
  * [6] K. Ishikawa and Y. Tobita, Finite-size corrections to Fermi’s Golden rule II: Quasi-stationary composite states, [arXiv:1607.08522].
  * [7] K. Ishikawa, T. Tajima, and Y. Tobita, Anomalous radiative transitions, PTEP 2015 (2015) 013B02 [arXiv:1409.4339].
  * [8] K. Ishikawa and K.Y. Oda, Particle decay in Gaussian wave-packet formalism revisted, PTEP 2018, (2018) 123B01 [arXiv:1809.04285].
  * [9] K. Ishikawa and T. Shimomura, Generalized S-matrix in mixed representations, Prog. Theor. Phys. 114 (2006) 1201 [arXiv:hep-ph/0508303].
  * [10] H. Araki, Y. Munakata, M. Kawaguchi, and T. Gotô, Quantum field theory of unstable particles, Prog. Theor. Phys. 17 (1957) 419.
  * [11] H. Araki, Y. Munakata, M. Kawaguchi, and T. Gotô, Quantum field theory of unstable particles, Prog. Theor. Phys. 18 (1957) 101.
  * [12] K. Yamamoto, Theory of unstable particles in the wave-packet-formalism, Prog. Theor. Phys. 20 (1958) 857.
  * [13] G.C. Cho, H. Kasari, and Y. Yamaguchi, The Time evolution of unstable particles, Prog. Theor. Phys. 90 (1993) 803.
  * [14] E. C. G. Stueckelberg, Relativistic Quantum Theory for Finite Time Intervals, Phys. Rev. 81 (1951) 130.
  * [15] M. Blasone, S. De Siena and C. Matrella, Wave packet approach to quantum correlations in neutrino oscillations, [arXiv:2104.03166].
  * [16] C. Giunti, Coherence and Wave Packets in Neutrino Oscillations, Found. Phys. Lett. 17 (2004) 103 [arXiv:hep-ph/0302026].
  * [17] C. Giunti, Neutrino wave packets in quantum field theory, JHEP 11 (2002) 017 [arXiv:hep-ph/0205014].
  * [18] C. Giunti and C. W. Kim, Coherence of neutrino oscillations in the wave packet approach, Phys. Rev. D 58 (1998) 017301 [arXiv:hep-ph/9711363]
  * [19] F. Ming, X-K. Song, J. Ling, L. Ye and D. Wang, Quantification of quantumness in neutrino oscillations, Eur. Phys. J. C 80 (2020) 275.
  * [20] P. Adamson et al. (MINOS Collaboration), Measurement of Neutrino Oscillations with the MINOS Detectors in the NuMI Beam, Phys. Rev. Lett. 101 (2008) 131802 [arXiv:0806.2237].
  * [21] A. B. Sousa (MINOS and MINOS+Collaborations), First MINOS+ data and new results from MINOS, AIP Conf. Proc. 1666 (2015) 110004.
  * [22] F. N. An et al. (Daya Bay Collaboration),New Measurement of Antineutrino Oscillation with the Full Detector Configuration at Daya Bay, Phys. Rev. Lett. 115 (2015) 111802 [arXiv:1505.03456]. J. Cao and K-B. Luk, An overview of the Daya Bay reactor neutrino experiment, Nucl. Phys. B 908 (2016) 62 [arXiv:1605.01502].
  * [23] Daya Bay Collaboration, Study of the wave packet treatment of neutrino oscillation at Daya Bay, Eur. Phys. J. C 77 (2017) 606 [arXiv:1608.01661].
  * [24] Z. Cheng, W. Wang, C. F.Wong and J. Zhang, Studying the neutrino wave-packet effects at medium-baseline reactor neutrino oscillation experiments and the potential benefits of an extra detector, Nucl. Phys. B 964, (2021) 115304 [arXiv:2009.06450].