[table]position=top

[type=editor, orcid=0009-0003-5583-0431] [1]

1]organization=India Meteorological Department, addressline=Ministry of Earth
Sciences, state=New Delhi, country=India

2]organization=Indian Institute of Information Technology, Allahabad,
city=Prayagraj, state=Uttar Pradesh, country=India

[cor1]Corresponding author

# Estimation of the Area and Precipitation Associated with a Tropical Cyclone
‘Biparjoy’ by using Image Processing

Shikha Verma<EMAIL_ADDRESS>Kuldeep Srivastava
<EMAIL_ADDRESS>Akhilesh Tiwari<EMAIL_ADDRESS>Shekhar
Verma<EMAIL_ADDRESS>[  [

###### Abstract

The rainfall associated with Topical Cyclone(TC) contributes a major amount to
the annual rainfall in India. Due to the limited research on the quantitative
precipitation associated with Tropical Cyclones (TC), the prediction of the
amount of precipitation and area that it may cover remains a challenge. This
paper proposes an approach to estimate the accumulated precipitation and
impact on affected area using Remote Sensing data. For this study, an instance
of Extremely Severe Cyclonic Storm, ‘Biparjoy’ that formed over the Arabian
Sea and hit India in 2023 is considered in which we have used the satellite
images of Integrated Multi-satellitE Retrievals for GPM (IMERG)-Late Run of
Global Precipitation Measurement (GPM). Image processing techniques were
employed to identify and extract precipitation clusters linked to the cyclone.
The results indicate that Biparjoy contributed a daily average rainfall of
53.14 mm/day across India and the Arabian Sea, with the Indian boundary
receiving 11.59 mm/day, covering an extensive 411.76 thousand square
kilometers. The localized intensity and variability observed in states like
Gujarat, Rajasthan, Madhya Pradesh, and Uttar Pradesh highlight the need for
tailored response measures, emphasizing the importance of further research to
enhance predictive models and disaster readiness, crucial for building
resilience against the diverse impacts of tropical cyclones.

###### keywords:

Tropical CycloneGlobal Precipitation MeasurementIMERG Image Processing

##  1 Introduction

Tropical Cyclones (TCs) are one of the most destructive natural disasters that
affect coastal regions worldwide (Ghosh and Mistri, 2023). The 5400 km
coastline of India is affected by TC of varying intensities. The North Indian
Ocean (NIO) is a basin located between Africa and the Indian subcontinent and
is one of the regions that experience TCs. The North Indian Ocean (NIO) basin
is divided into two sub-basins, the Bay of Bengal ( Bay of Bengal (BoB)) and
the Arabian Sea. The Bay of Bengal sub-basin is more active than the Arabian
Sea, with the majority of the TCs forming in this region. The NIO basin
experiences cyclones from April to December, with the peak season occurring
from September to November. The cyclones in this region are classified by to
the India Meteorological Department ( India Meteorological Department (IMD))
scale, which takes into account the maximum sustained wind speed near the
center of the cyclone.

The impacts of TCs on the NIO region are significant, including damage to
infrastructure, loss of life, and economic losses. The frequency and intensity
of TCs in the NIO have been observed to be increasing in recent decades. This
increase in TC activity has been attributed to the rising sea surface
temperatures, which provide favorable conditions for the formation and
intensification of TCs (Ali, 1996; Michaels et al., 2006). The NIO region is
also vulnerable to storm surges, which are caused by the combination of high
winds and low atmospheric pressure during a TC. The impact of storm surges is
often more severe than the impact of the strong winds associated with TCs
(Jakobsen and Madsen, 2004; Hubbert et al., 1991; Needham et al., 2015).
Though the disaster caused by a TC cannot be avoided, its impact can be
greatly mitigated through a variety of management approaches, such as
response, recovery, prevention/reduction, and preparedness (Hoque et al.,
2017). However, this makes it essential to understand the patterns and
behavior of TCs in the NIO basin to improve early warning systems.

Rainfall associated with TCs has a significant impact on the environment and
human activities in many regions of the world. These intense weather systems
often cause severe flooding, landslides, and other hazards that can cause
widespread damage to infrastructure, agriculture, and human settlements.
Therefore, the study of rainfall patterns associated with TC is crucial for
understanding the characteristics and behavior of these storms and their
potential impact on the affected areas. This may also allow meteorologists to
prepare precipitation climatology by determining the average downpour over the
globe. The long-term precipitation averages may be compared against short-term
precipitation events to calculate deviation. These deviations or departures
from climatology are an important aspect of interpreting the changing weather
patterns. This will help to assess the severity of weather phenomena that are
accompanied by severe precipitation occurrences such as TC.

It was observed that in 50 years (1970-2019), floods and tropical cyclones
have been two major disasters causing mortality (Ray et al., 2021). One of the
most significant impacts of tropical cyclones is the abundant rainfall
produced by them (Marchok et al., 2007; Rogers et al., 2009). While forecasts
of TC track and intensity have improved significantly, less attention has been
paid to improve the estimate and forecast of the quantitative precipitation
from TCs (Marchok et al., 2007; BISWAS et al., 2013; Kidder et al., 2005). It
has been observed that the precipitation associated with TCs accounts for
$\sim 5-6\%$ over the globe, $\sim 6-9\%$ of total precipitation over the
Tropics, and about $5\%$ over the NIO which includes both the Bay of Bengal
(BoB) and the Arabian Sea (AS) (Bhatla et al., 2020; Jiang and Zipser, 2010;
Lin et al., 2015).

Because of the dense population in peninsular India, the effects of TCs
rainfall are more severe over the NIO (Singh et al., 2019). The spatial
distribution of rainfall (Osuri et al., 2012) and the time of occurrence of
rainfall (Mohapatra et al., 2013) are sensitive to cyclone tracks. The IMD and
Joint Typhoon Warning Centre (JTWC), USA monitors the TCs over NIO to estimate
the cyclone’s path, intensity, landfall location, and time of landfall with
the help of an optimal observational network, including satellite, radar,
surface, and upper air observations (Mohapatra et al., 2012).

The prediction of rainfall associated with a TC is challenging because of TC’s
dynamic characteristics. Rainfall estimates based on manual observation,
radar, Numerical Weather Prediction (NWP), and other remote sensors present
intriguing opportunities for improvement. However, there has been minimal
prior research confirming the model QPF performance in TCs. One of the key
impediments to improving QPF in TCs is the lack of a description of the
distribution of rain in space and time (Tang and Matyas, 2018).

Kidder et al. (2000) estimated the rain rate in TC with track forecasts using
Tropical Rainfall Potential (TRaP) and calculated precipitation potential from
Advanced Microwave Sounding Unit ( Advanced Microwave Sounding Unit (AMSU))
data that appeared to be useful in forecasting the heavy precipitation
associated with a landfalling TC. Lonfat et al. (2004) studied the
precipitation distributions in TCs from National Aeronautics and Space
Administration (NASA)’s TRMM to understand the precipitation structure of
tropical cyclones and the dynamics governing the rainfall structure. In an
attempt by Luitel et al. (2018), the performance of NWP was evaluated in
forecasting rainfall associated with 15 U.S. landfalling TCs during the
2007-2012 period. It was observed that there was a substantial error in the
precipitation prediction associated with a TC. Hence, prior to the prediction,
it is important to have a detailed study on quantitative analysis of rainfall
observed during the lifetime of a TC. Understanding the quantitative aspects
of rainfall, such as its intensity, volume, and distribution, meteorologists
can better predict and assess the potential risks associated with these
cyclones, enabling them to issue timely warnings and take appropriate
mitigation measures. Additionally, it is crucial for climate studies and
research. Climate change is known to influence the behavior of cyclones,
including their rainfall characteristics. By analyzing and quantifying
rainfall patterns over time, scientists can study the impact of climate change
on tropical cyclones and assess any potential changes in their behavior, such
as increased rainfall rates or shifts in the distribution of rainfall. By
analyzing historical rainfall data and patterns, scientists can develop models
that provide insights into the potential impact of rainfall associated with
tropical cyclones, enabling authorities to issue timely warnings, plan
evacuation strategies, and allocate necessary resources. Furthermore, it helps
in assessing the hydrological risks associated with tropical cyclones, aiding
in the design and implementation of appropriate infrastructure, drainage
systems, and flood mitigation measures.

One of the basic requirements of rainfall study associated with any weather
phenomenon is to identify the area of proneness on the basis of vulnerability
and hazard. Since areas affected by cyclonic rainfall can be influenced by a
number of factors, such as environmental humidity, low-level vorticity,
vertical wind shear, TC latitude, and TC intensity (Lin et al., 2015), it is
important for any country to identify the prone locations in order to lower
down the adverse impact of the cyclone. In this context, the government
organizations of India such as the Building Materials and Technology Promotion
Council (BMTPC) and National Disaster Management Authority (NDMA) have listed
vulnerable districts of different states of India based on cyclone-proneness.
The IMD has identified districts for which cyclone warnings must be issued.
Taking it to further, Mohapatra et al. (2012) have attempted to segregate the
cyclone-prone districts of India into very high, high, moderate, and low
categories of proneness. The categorization has been done on the basis of
Probable Maximum Winds (PMW), Probable maximum storm surge (PMSS), and
Probable maximum precipitation (PMP) for a day over any station in the
district considering the other parameters such as enhanced rainfall in
association with monsoon systems.

The manual observation of rainfall in the basin is not possible. Moreover, the
weather radars installed at coasts have range limitations. Due to limited
observations, it is difficult to accurately compute the actual rainfall that
occurred over the basin during a cyclonic event. Fortunately, an earth-
observing satellite Global Precipitation Measurement (GPM), provides frequent
estimates of precipitation at a global scale including rainfall occurred in
the basin during a cyclonic period. The Integrated Multi-satellitE Retrievals
for GPM is designed to compensate for the limited sampling available from
single Low-Earth-Orbit (LEO) -satellites by using as many LEO-satellites as
possible, and then augmenting with Geosynchronous-Earth-Orbit (
Geosynchronous-Earth-Orbit (GEO)) Infrared (IR) estimates. The IMERG system
runs twice in near-real time providing “Late” multi-satellite product 14 hr
after observation time and is appropriate for daily and longer applications.

In this context, an attempt has been made to provide quantitative
precipitation estimation using image processing techniques over IMERG products
and infer the complex relation between cyclones and associated precipitation.
An instance of the Extremely Severe Cyclonic Storm, ‘Biparjoy’ is taken to
explain the approach used in this study. The Biparjoy Cyclone occurred over
the Arabian Sea and hit India in 2023. The daily average precipitation
associated with Biparjoy over the Arabian Sea and India has been estimated
along with the statistics of states affected by it.

The document begins by explaining the identified problem in Section 2,
shedding light on it’s importance. Moving on to Section 3.1, it delves into
the specific geographic focus of the study, while Section 3.2 provides details
about the data used. In Section 4, the step-by-step process employed to
address the problem is detailed. The results are presented in Section 5,
followed by an in-depth discussion in Section 6. Finally, the conclusion is
presented in the Section 7.

##  2 Problem Description

Currently, there is a significant knowledge gap regarding TC-induced rainfall
patterns, hindering precise predictions and effective mitigation strategies.
The lack of specific TC-related rainfall data poses challenges in assessing
the intensity and spatial distribution of precipitation within these cyclonic
systems. To address this, a novel approach utilizing image processing
techniques has been proposed to estimate TC-associated rainfall and its
spatial distribution. The segmentation of rainfall areas within the TC system
enables a quantitative assessment of rainfall extent and intensity. The
significance of this research lies in improving our understanding of TC-
induced rainfall and enhancing the accuracy of rainfall predictions, aiding
meteorologists, disaster management agencies, and policymakers in making
informed decisions. Furthermore, this research focuses specifically on the
Indian region, where tropical cyclones frequently make landfall and have the
potential to cause significant damage. By developing a model that can
accurately estimate the rainfall patterns associated with TCs over India, this
study contributes to the existing knowledge base and improves the overall
understanding of cyclone dynamics in this region. This research can also serve
as a foundation for further studies and advancements in TC rainfall
estimation, ultimately leading to more accurate predictions and better
preparedness for cyclone-related rainfall events. However, the reliance on
satellite imagery introduces limitations such as potential cloud cover, image
resolution, and real-time data availability, which can affect the accuracy and
timeliness of rainfall estimates. Addressing these limitations is crucial for
robust and accurate TC rainfall estimations.

##  3 Study Area and Dataset

Figure 1: Study area

###  3.1 Study Area

The study area chosen for this research is the Arabian Sea (AS) of the NIO
basin as shown in shaded region of Figure 1. The AS is divided into 6 zones
namely NWA, NEA, WCA, ECA, SWA, and SEA. The region is known to be prone to
tropical cyclones, particularly during the pre-monsoon and post-monsoon
seasons, and therefore, it is important to understand the behavior of these
cyclones to better prepare and mitigate their impacts. The cyclogenesis of
tropical cyclone Biparjoy of Extremely Severe Cyclonic Storm ( Extremely
Severe Cyclonic Storm (ESCS)) category was observed in AS on 06-06-2023, which
was the first Cyclonic Storm (CS) over the Arabian Sea in the year 2023. As
per the report of IMD, it crossed Saurashtra and Kutch coast close to Jakhau
Port near latitude 23.28°N and longitude 68.56°E with Maximum Sustained Wind
(MSW) of 65 knots gusting to 75 knots on 15-06-2023. Finally, it weakened into
a depression ( Depression (D)) over South Rajasthan and adjoining North
Gujarat on 18-06-2023 and into a WML over central parts of Northeast Rajasthan
and neighbourhood on 19-06-2023 at 03 UTC.

###  3.2 Dataset

The GPM_3IMERGHHL (GPM IMERG Late Precipitation L3 Half Hourly) dataset has
been selected as the primary data set for this study. This choice is motivated
by the dataset’s capability to provide comprehensive and detailed information
on precipitation dynamics. Utilizing advanced algorithms such as GPROF2017,
the dataset offers half-hourly precipitation estimates with a spatial
resolution of 0.1° x 0.1°. Notably, the inclusion of "late" observations,
occurring around 14 hours after the initial observation time, adds temporal
depth to the dataset. The meticulous processing, including gridding,
intercalibration to the GPM Combined Ku Radar-Radiometer Algorithm (CORRA)
product, and bias correction, enhances the dataset’s accuracy and reliability,
making it an ideal choice for the study’s comprehensive analysis of
precipitation patterns (NASA GES DISC, 2019).

In adherence to the WMO standards, this study considers the Date of
Acquisition (DOA) for 24-hour precipitation from 07-06-2023 to 20-06-2023
recorded at 03 UTC to ensure a consistent and uniform observation protocol.
The detail of the dataset along with Date_id (D_ID) is listed in Table 1 and
will be used as a nomenclature in the further sections.

Table 1: Details of the dataset used SN |  Date of Acquisition (DOA) | D_ID
---|---|---
1 | 07-06-2023 | D1
2 | 08-06-2023 | D2
3 | 09-06-2023 | D3
4 | 10-06-2023 | D4
5 | 11-06-2023 | D5
6 | 12-06-2023 | D6
7 | 13-06-2023 | D7
8 | 14-06-2023 | D8
9 | 15-06-2023 | D9
10 | 16-06-2023 | D10
11 | 17-06-2023 | D11
12 | 18-06-2023 | D12
13 | 19-06-2023 | D13
14 | 20-06-2023 | D14

##  4 Methodology

The proposed method uses satellite imagery to analyze rainfall and determine
the extent of the Biparjoy cyclone’s coverage. By processing the satellite
data using image processing techniques, we can identify specific rainfall
patterns and understand the impact of the cyclone more accurately. This
approach overcomes the limitations of traditional methods and provides a
clearer picture of the cyclone’s effects, even in areas with limited ground-
based observations. Hence, it offers a more comprehensive understanding of the
cyclone’s impact and enables improved decision-making and response strategies.

The step-wise workflow of the proposed approach is shown in Figure 2. It shows
that the data is downloaded which is pre-processed for further processing.
Then, the mask of the precipitation blob has been generated from which the
major blob was extracted based on the connected component. The precipitation
from the extracted cluster has been computed for the study site. Based on the
spatial resolution of IMERG data of rainfall, the total area affected by
rainfall associated with cyclone ‘Biparjoy’ in India has been computed. A
comprehensive explanation of each step is given below.

Figure 2: Workflow of the proposed approach

###  4.1 Data download

The IMERG product GPM_3IMERHHL v06 is downloaded from the data server of NASA.
The details of the dataset are given in section 3.2. IMERG (Integrated Multi-
satellite Retrievals for GPM) is a unified satellite precipitation product
that provides estimated surface precipitation over most of the globe. It is a
product of Global Precipitation Measurement (GPM), which is an international
satellite mission jointly launched by NASA and Japan Aerospace Exploration
Agency (JAXA) to provide unified precipitation measurements of rain and snow
every three hours.

###  4.2 Pre-processing

In this section, the pre-processing steps for the downloaded data are
explained. The pre-processing steps ensure that the data is in the appropriate
format and ready for further analysis. These steps are important in any
analysis involving geospatial data, as they can significantly impact the
accuracy and reliability of the results.

The first step in pre-processing is re-projection, which involves warping the
images from their original projection space to a selected projection space.
This process requires the calculation of a geometric transformation model,
which accounts for the relative geometry of the two projections and local
scale and resolution changes in the images. The amount of distortion incurred
during re-projection depends on these factors. In this study, the projection
of downloaded raster data is converted from sinusoidal to geographic latitude-
longitude for further processing (Steinwand et al., 1995).

The next step in pre-processing is spatial subsetting, which involves
extracting the region of interest (ROI). This process returns a new object
with only features that are spatially related to the Region of Interest (ROI)
(Ayoobi and Tangestani, 2017). This step is important because it reduces the
data size and processing time, making it easier to analyze the data. Figure 3
shows the images for D1 to D14 after re-projection and spatial subsetting.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

Figure 3: Pre-processed images for the D_ID (a)-(n): D1-D14

###  4.3 Masking

After the pre-processing step, a two-step thresholding approach is employed to
filter out irrelevant information and highlight the crucial features in the
satellite imagery. Initially, pixels with values exceeding 0.9 mm are masked
out, focusing on regions where precipitation is significant. Subsequently, a
finer threshold of 0.3 mm is employed, emphasizing the areas with a more
pronounced cyclonic structure. This sequential approach ensures that the
masking strategy evolves from a broader emphasis on significant precipitation
to a more nuanced focus, particularly when the cyclone’s eye becomes
prominent. The outcome of this process is a binary mask image, which contains
pixels with either a value of 0 or 1, representing the areas that meet the
specified condition and those that do not. To visually represent the
effectiveness of this preprocessing step, the masked images for the D1 to D14
are shown in Figure 4.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

Figure 4: Cloud mask for D_ID (a)-(n): D1-D14

###  4.4 Identification and Extraction of Cyclonic Precipitation Cluster

In this step of the study, the focus is on labeling the results obtained from
the previous step based on the area they cover. The objective is to identify
the precipitation clusters associated with the Biparjoy cyclone, which can
then be further processed. The employed methodology involves leveraging
Connected Component Analysis ( Connected Component Analysis (CCA)). This
technique is utilized to effectively identify and label connected regions
within a binary image, enabling a detailed examination of the precipitation
patterns associated with the cyclone. It involves finding the objects in a
binary image by creating a labeled image where pixels with identical integer
values belong to the same class of object. The algorithm works by scanning the
binary image pixel by pixel and checking whether each pixel is connected to
any of its neighbors. A pixel is considered connected to its neighbors if they
have the same value, and they are adjacent to each other either horizontally
or vertically. When a group of connected pixels is identified, they are
assigned a unique label, and the process is repeated until all pixels in the
image have been labeled. The result of the CCA algorithm is a labeled image,
where each cluster of connected pixels is assigned a unique label. This
labeled image can then be used for further analysis and processing.
Mathematically, the CCA algorithm can be represented by the following
equation:

$L=bwlabel(BW)$

where L is the labeled image, BW is the binary image, and bwlabel is the
function that performs the CCA operation. The function bwlabel assigns a
unique integer label to each connected cluster of pixels in the binary image
BW.

In the present study, the binary image represents the precipitation clusters
associated with the cyclone. The False value in the binary image represents
background pixels, while the True value represents object pixels.

Moreover, to ascertain the precise association of identified precipitation
clusters with the Biparjoy cyclone and to mitigate any potential ambiguity
arising from other weather events in the region, we are incorporating the best
track information of 03 UTC provided by the IMD, as illustrated in Figure 5.
This meticulous consideration of the Biparjoy cyclone’s best track is
essential for the accurate labeling of the identified precipitation clusters,
laying a crucial foundation for the subsequent in-depth analysis of the
associated rainfall patterns.

The extracted images of cyclonic precipitation clusters so obtained are shown
in Figure 6. These images serve as pivotal visual aids for further analysis
and interpretation of the cyclonic activity attributed to the Biparjoy. Moving
forward, Figure 7, presents a vector image of the study area projected onto
the extracted cyclonic precipitation clusters, which effectively illustrates
the distribution of precipitation in various Indian states.

Figure 5: Best Track of Biprajoy Cyclone

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

Figure 6: Extracted Cyclonic Precipitation Clusters for the D_ID (a)-(n):
D1-D14

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

Figure 7: Visualization of Extracted Cyclonic precipitation cluster over the
boundary of India for the D_ID (a)-(n): D1-D14

Having identified and labeled the cyclonic precipitation clusters, the
subsequent step involves computing pertinent statistics associated with these
clusters. Specifically, the total precipitation associated with the Biparjoy
cyclone is computed by summing up the precipitation values within the
identified clusters. Additionally, the total area affected by the clusters and
the total area of the affected states is also computed for each day. These
computed metrics provide valuable insights into the cyclone’s impact on the
affected regions, offering essential information for disaster response and
management efforts.

##  5 Result

In this section, we estimate the total rainfall associated with cyclone
Biparjoy on each day along with the area covered by it. Moreover, the centroid
of extracted precipitation cluster and Maximum Rainfall Point are computed.
Thereafter, the result is analyzed and discussed in detail.

###  5.1 Pixel Count

The pixel counts in the processing of IMERG product containing the number of
pixels in the mask and the number of pixels in extracted precipitation
clusters after applying image processing techniques from D1 to D14 are shown
in Figure 8.

Figure 8: Pixel Count of Mask and Extracted Cyclonic Precipitation Cluster

The analysis of the pixel data in conjunction with the distinctive
characteristics of Cyclone Biparjoy reveals intricate relationships between
meteorological phenomena and the observed variations in precipitation
clusters. On analyzing 8, the findings reveal significant insights into the
spatial distribution of precipitation clusters associated with the Biparjoy
cyclone. Notably, the count of pixels within the binary mask consistently
exceeds the count of pixels in the extracted precipitation cluster. This
discrepancy is indicative of areas within the mask that do not contribute to
the identified precipitation clusters. Further investigation indicates that a
subset of pixels in the mask is associated with convection patterns related to
the cross Equatorial southwest monsoon flow over the southeast Arabian Sea.
Additionally, the rainfall observed over the Kerala and Lakshadweep area is
attributed to the onset of the southwest monsoon, which commenced over most
parts of Kerala and the entire Lakshadweep on June 8th.

The temporal evolution of the pixel counts provides valuable information about
the dynamic nature of precipitation clusters throughout the observed period.
On D9, we observe a notable increase in the number of pixels within the binary
mask, suggesting an intensification of weather patterns. This escalation is
likely attributed to a weather phenomenon occurring over North-Eastern part of
India. Concurrently, the count of pixels in the extracted precipitation
cluster remains relatively consistent, signifying a concentration of
precipitation in specific regions. Furthermore, the diurnal fluctuations in
convective clouds, further emphasize the dynamic nature of the precipitation
clusters, reflected in the pixel counts.

Following landfall, the sustained pixel counts within the binary mask and
extracted precipitation cluster illustrate the cyclone’s persistence even
after making landfall. These pixel trends reflect the post-landfall intensity
and gradual weakening till midnight of June 16th, providing a visual
representation of the cyclone’s behavior.

Consistently, the pixel data reveals variations, underscoring the complexity
of Cyclone Biparjoy and the challenges in predicting its path. The
interconnectedness of pixel variations and meteorological characteristics
highlights the potential for advanced methodologies in pixel segmentation and
the integration of additional data sources for a more comprehensive
understanding of cyclonic events. This holistic approach contributes to
improved forecasting and enhanced disaster management capabilities.

###  5.2 Estimation of Rainfall

Since the IMERG product acquires data from sparse locations, the rainfall data
of cyclones over the ocean is available. After pre-processing, the total rain
due to the Biparjoy cyclone has been computed. Knowing the total rainfall over
India after landfall is important to prepare the rainfall climatology.
Therefore, the day-wise total rainfall and the corresponding area are computed
before and after landfall over India by using Simple Mean Method. Though the
Simple Mean method is easy to use and provides a quick estimate of the average
value, it has limitations in situations where the data is highly skewed or
contains outliers.

On the basis of extracted cyclonic precipitation cluster, the average daily
rainfall is calculated over India by using Simple Mean. In this calculation,
rainfall measurements greater than 0.1mm are considered significant, while
measurements less than or equal to 0.1mm are considered as trace. The result
is shown in Figure 9.

Figure 9: Average Daily Rainfall due to Biparjoy Cyclone over India Figure 10:
Rainfall received by different states of India

Figure 9 provides an insight into the variation of average daily rainfall
during the Biparjoy cyclone event and highlights the importance of considering
both the extracted cyclonic precipitation cluster and the area over which the
rainfall is averaged when analyzing rainfall patterns.

The daily average rainfall due to Cyclone Biparjoy, computed from the
extracted precipitation clusters, stands at 53.14 mm/day, encompassing both
India and the Arabian Sea region, whereas the Indian boundary received an
average of 11.59 mm/day. This aggregated metric serves as a key indicator of
the overall intensity and extent of precipitation influenced by Cyclone
Biparjoy during the analyzed timeframe. Upon closer examination of the broader
geographical scope, the impact on India as a whole is characterized by
significant variations in rainfall patterns across different states.

During the influence of Cyclone Biparjoy, four states were significantly
affected namely Gujarat, Rajasthan, Madhya Pradesh, and Uttar Pradesh
experiencing distinct patterns of precipitation as depicted in Figure 10.
Cyclone Biparjoy made landfall on D6 along the coast of Gujarat, marking a
significant meteorological event. The impact of the cyclone on Gujarat was the
most pronounced impact across 6 days, marked by a substantial 31.25 mm/day
downpour on D7 and an average daily rainfall of 11.03 mm/day. This rainfall
persisted in Gujarat on the subsequent day, showcasing the sustained influence
of Cyclone Biparjoy on the region. Intriguingly, on D9, a temporary respite
from rainfall was observed over the entire Indian region. However, this
reprieve was short-lived as the cyclone resumed affecting Gujarat and extended
its impact to Rajasthan from the following day onwards. The meteorological
data indicates a noteworthy shift in the cyclone’s trajectory and influence.
The absence of recorded rainfall over India on D9 signifies a brief interlude
in the cyclone’s activity, offering a temporary relief from precipitation
across the country. Nevertheless, the subsequent resurgence of rainfall,
specifically impacting Gujarat and Rajasthan, underscores the dynamic nature
of Cyclone Biparjoy’s trajectory and its ability to affect diverse regions
over time. Rajasthan, spanning 5 days, experienced peak precipitation on D13,
recording a notable 42.89 mm/day rainfall, with an average daily rainfall of
8.21 mm. Madhya Pradesh, impacted over 2 days, observed its highest
precipitation of 10.05 mm on D13, featuring an average daily rainfall of 2.62
mm/day. Similarly, Uttar Pradesh, affected for 1 day, witnessed peak rainfall
of 19.15 mm/day on D14, displaying a comparatively lower average daily
rainfall of 1.37 mm/day. Notably, Rajasthan received the highest single-day
rainfall among the four states, recording 42.89 mm/day on D13. This
information emphasizes the localized intensity and variability of
precipitation during the cyclone’s impact on different regions. The state-
specific statistics provide a comprehensive and detailed assessment of
Biparjoy’s influence on regional precipitation patterns. These detailed
insights play a pivotal role in formulating well-informed disaster management
and response strategies tailored to the specific impact experienced by each
state, enhancing overall preparedness and resilience.

###  5.3 Estimation of Area Covered

In this step, we conducted a thorough analysis of the spatial distribution of
rainfall during the Biparjoy cyclone, aiming to understand its impact on
different regions. To achieve this, we calculated the total area covered by
the precipitation cluster over India. The results are visually presented in
Figure 11, providing valuable insights into the geographical extent of
rainfall associated with the cyclone.

The comprehensive area covered by the Biparjoy cyclone was determined by
unifying all the extracted raster images. The union operation allowed us to
merge these individual clusters, providing a consolidated view of the spatial
extent of precipitation associated with the cyclone. This approach ensures
that the cumulative impact area is accurately captured, considering the
evolving nature of the cyclone over time. This methodology contributes to a
comprehensive analysis of the cyclone’s influence on various regions, laying
the groundwork for a thorough understanding of its spatial and temporal
characteristics.

Figure 11: Area Covered by Biparjoy Cyclone over Arabian Sea and India Figure
12: Total Area Covered by Extracted Precipitation Cluster of Biparjoy Cyclone
and Over India

In the course of our analysis of Cyclone Biparjoy’s impact on different
regions, we meticulously examined the daily coverage area across Gujarat,
Rajasthan, Madhya Pradesh, and Uttar Pradesh. The visual representation is
shown in 12.

Our findings reveal that the Biparjoy covered a substantial 411.76 thousand
square kilometers in total, with Gujarat experiencing an impact over 154.75
thousand square kilometers, followed by Rajasthan with 194.80 thousand square
kilometers. Madhya Pradesh and Uttar Pradesh were also significantly affected,
with impact areas of 40.78 and 21.43 thousand square kilometers, respectively.

Throughout the progression of Cyclone Biparjoy, distinct patterns of
precipitation coverage were observed. On D7 and D8, a notable area of Gujarat,
totaling 39.1 thousand sq. km, experienced the impact of the cyclone.
Subsequently, the influence intensified on D10 and D11, extending to cover
substantial regions in Gujarat, Rajasthan, and Madhya Pradesh. The zenith of
the cyclone’s impact was observed on D11, with widespread precipitation
exceeding 150.7 thousand sq. km in Gujarat and 57.3 thousand sq. km in
Rajasthan. This impactful trend persisted on D12 and D13, revealing
significant coverage in Gujarat, Rajasthan, and Madhya Pradesh.

Madhya Pradesh, in particular, witnessed substantial coverage from D12 to D14,
recording precipitation areas of 207.0 thousand sq. km, 138.9 thousand sq. km,
and 77.9 thousand sq. km, respectively. Even on the concluding day (D14), as
the cyclone’s influence gradually diminished, a considerable area remained
covered in Rajasthan, Madhya Pradesh, and Uttar Pradesh, highlighting the
lasting effect of Cyclone Biparjoy across these regions. Notably, Rajasthan
sustained coverage of 40.8 thousand sq. km, while Uttar Pradesh experienced
precipitation over an area of 21.4 thousand sq. km. This detailed analysis
provides valuable insights into the temporal progression of the cyclone and
its varying impact on different states, contributing to our understanding for
further research and disaster management strategies.

##  6 Discussion

In this study, we estimated the quantitative precipitation associated with
cyclone Biparjoy which provided insights into the spatial distribution of the
rainfall and its impact on the affected areas. We used image processing
techniques on IMERG raster data to identify and extract precipitation clusters
associated with the Biparjoy cyclone. Specifically, we computed the total
rainfall for each day and also determined the area covered by the cyclone. We
then proceeded to conduct a detailed analysis of the results obtained.

The findings of this study provide a significant contribution to the
understanding of Cyclone Biparjoy’s impact on precipitation patterns across
India. The daily average rainfall of 53.14 mm/day reflects the cyclone’s
substantial influence, with Gujarat experiencing a particularly pronounced
impact, marked by a peak daily rainfall of 31.25 mm on D7. The dynamic
trajectory of the cyclone resulted in varying levels of impact on different
states, showcasing the complexity of cyclonic behavior. The spatial
distribution analysis further underscores the extensive coverage area of
411.76 thousand square kilometers, emphasizing the far-reaching impact of
Cyclone Biparjoy.

The temporal progression of the cyclone’s influence, as revealed through this
study, sheds light on the resilience and adaptability of disaster management
strategies. The localized intensity and variability observed in states like
Gujarat, Rajasthan, Madhya Pradesh, and Uttar Pradesh highlight the need for
tailored response measures to address the diverse impacts of cyclones on
different regions. The short-lived reprieve from rainfall on D9 and the
subsequent resurgence underscore the challenges in predicting and managing
cyclonic behavior.

This study’s insights provide a foundation for future research and contribute
to the development of more targeted disaster preparedness and response
strategies. The nuanced understanding of cyclonic dynamics, as presented in
this research, is crucial for enhancing the overall resilience of communities
and infrastructure in the face of such meteorological events.

##  7 Conclusion

Our study highlighted the significant impact of tropical cyclone Biparjoy on
precipitation over India and the need for further research to better
understand the dynamics and implications for regional precipitation. This
study urges further investigations to improve predictive models, enhance
disaster readiness, and develop specific responses. Such efforts are vital for
building resilience against the diverse impacts of tropical cyclones in the
changing climate scenario.

## Acknowledgement

The author would like to express my heartfelt gratitude to the India
Meteorological Development (IMD) and the Indian Institute of Information
Technology, Allahabad (IIIT Allahabad) for their invaluable support and
contributions to this journal. Their guidance, resources, and assistance have
been instrumental in the successful completion of this research.

## Abbreviations

AMSU

    Advanced Microwave Sounding Unit
BMTPC

    Building Materials and Technology Promotion Council
BoB

    Bay of Bengal
CCA

    Connected Component Analysis
CS

    Cyclonic Storm
D

    Depression
DOA

    Date of Acquisition
ESCS

    Extremely Severe Cyclonic Storm
GEO

    Geosynchronous-Earth-Orbit
GPM

    Global Precipitation Measurement
IMD

    India Meteorological Department
IMERG

    Integrated Multi-satellitE Retrievals for GPM
IR

    Infrared
JAXA

    Japan Aerospace Exploration Agency
JTWC

    Joint Typhoon Warning Centre
LEO

    Low-Earth-Orbit
MSW

    Maximum Sustained Wind
NASA

    National Aeronautics and Space Administration
NDMA

    National Disaster Management Authority
NIO

    North Indian Ocean
NWP

    Numerical Weather Prediction
ROI

    Region of Interest
TC

    Tropical Cyclones
TCs

    Tropical Cyclones

## References

  * Ghosh and Mistri (2023) S. Ghosh, B. Mistri,  Cyclone-induced coastal vulnerability, livelihood challenges and mitigation measures of matla–bidya inter-estuarine area, indian sundarban,  Natural Hazards 116 (2023) 3857–3878.
  * Ali (1996) A. Ali, Vulnerability of Bangladesh to climate change and sea level rise through tropical cyclones and storm surges, Springer, 1996.
  * Michaels et al. (2006) P. J. Michaels, P. C. Knappenberger, R. E. Davis,  Sea-surface temperatures and tropical cyclones in the atlantic basin,  Geophysical Research Letters 33 (2006).
  * Jakobsen and Madsen (2004) F. Jakobsen, H. Madsen,  Comparison and further development of parametric tropical cyclone models for storm surge modelling,  Journal of Wind Engineering and Industrial Aerodynamics 92 (2004) 375–391.
  * Hubbert et al. (1991) G. D. Hubbert, G. J. Holland, L. M. Leslie, M. J. Manton,  A real-time system for forecasting tropical cyclone storm surges,  Weather and Forecasting 6 (1991) 86–97.
  * Needham et al. (2015) H. F. Needham, B. D. Keim, D. Sathiaraj,  A review of tropical cyclone-generated storm surges: Global data sources, observations, and impacts,  Reviews of Geophysics 53 (2015) 545–591.
  * Hoque et al. (2017) M. A.-A. Hoque, S. Phinn, C. Roelfsema, I. Childs,  Tropical cyclone disaster management using remote sensing and spatial analysis: A review,  International journal of disaster risk reduction 22 (2017) 345–354.
  * Ray et al. (2021) K. Ray, R. Giri, S. Ray, A. Dimri, M. Rajeevan,  An assessment of long-term changes in mortalities due to extreme weather events in india: A study of 50 years’ data, 1970–2019,  Weather and Climate Extremes 32 (2021) 100315.
  * Marchok et al. (2007) T. Marchok, R. Rogers, R. Tuleya,  Validation schemes for tropical cyclone quantitative precipitation forecasts: Evaluation of operational models for us landfalling cases,  Weather and forecasting 22 (2007) 726–746.
  * Rogers et al. (2009) R. Rogers, F. Marks, T. Marchok,  Tropical cyclone rainfall,  Encyclopedia of hydrological sciences 3 (2009).
  * BISWAS et al. (2013) H. BISWAS, P. Kundu, D. Pradhan,  A case study for cyclone ‘aila’for forecasting rainfall using satellite derived rain rate data,  Mausam 64 (2013) 77–82.
  * Kidder et al. (2005) S. Q. Kidder, J. A. Knaff, S. J. Kusselson, M. Turk, R. R. Ferraro, R. J. Kuligowski,  The tropical rainfall potential (trap) technique. part i: Description and examples,  Weather and Forecasting 20 (2005) 456–464.
  * Bhatla et al. (2020) R. Bhatla, R. Raj, R. Mall, Shivani,  Tropical cyclones over the north indian ocean in changing climate,  Techniques for disaster risk management and mitigation (2020) 63–76.
  * Jiang and Zipser (2010) H. Jiang, E. J. Zipser,  Contribution of tropical cyclones to the global precipitation from eight seasons of trmm data: Regional, seasonal, and interannual variations,  Journal of climate 23 (2010) 1526–1543.
  * Lin et al. (2015) Y. Lin, M. Zhao, M. Zhang, Tropical cyclone rainfall area controlled by relative sea surface temperature. nat. commun., 6, 6591, 2015.
  * Singh et al. (2019) K. Singh, M. Mandal, P. K. Bhaskaran,  Impact of radiance data assimilation on the prediction performance of cyclonic storm sidr using wrf-3dvar modelling system,  Meteorology and Atmospheric Physics 131 (2019) 11–28.
  * Osuri et al. (2012) K. K. Osuri, U. Mohanty, A. Routray, M. A. Kulkarni, M. Mohapatra,  Customization of wrf-arw model with physical parameterization schemes for the simulation of tropical cyclones over north indian ocean,  Natural Hazards 63 (2012) 1337–1359.
  * Mohapatra et al. (2013) M. Mohapatra, D. Nayak, R. Sharma, B. Bandyopadhyay,  Evaluation of official tropical cyclone track forecast over north indian ocean issued by india meteorological department,  Journal of Earth System Science 122 (2013) 589–601.
  * Mohapatra et al. (2012) M. Mohapatra, B. Bandyopadhyay, A. Tyagi,  Best track parameters of tropical cyclones over the north indian ocean: A review,  Natural Hazards 63 (2012) 1285–1317.
  * Tang and Matyas (2018) J. Tang, C. Matyas,  A nowcasting model for tropical cyclone precipitation regions based on the trec motion vector retrieval with a semi-lagrangian scheme for doppler weather radar,  Atmosphere 9 (2018) 200.
  * Kidder et al. (2000) S. Q. Kidder, M. D. Goldberg, R. M. Zehr, M. DeMaria, J. F. Purdom, C. S. Velden, N. C. Grody, S. J. Kusselson,  Satellite analysis of tropical cyclones using the advanced microwave sounding unit (amsu),  Bulletin of the American Meteorological Society 81 (2000) 1241–1260.
  * Lonfat et al. (2004) M. Lonfat, F. D. Marks Jr, S. S. Chen,  Precipitation distribution in tropical cyclones using the tropical rainfall measuring mission (trmm) microwave imager: A global perspective,  Monthly Weather Review 132 (2004) 1645–1660.
  * Luitel et al. (2018) B. Luitel, G. Villarini, G. A. Vecchi,  Verification of the skill of numerical weather prediction models in forecasting rainfall from us landfalling tropical cyclones,  Journal of hydrology 556 (2018) 1026–1037.
  * Mohapatra et al. (2012) M. Mohapatra, G. Mandal, B. Bandyopadhyay, A. Tyagi, U. Mohanty,  Classification of cyclone hazard prone districts of india,  Natural hazards 63 (2012) 1601–1620.
  * NASA GES DISC (2019) NASA GES DISC, GPM IMERG late precipitation L3 1 day 0.1 degree x 0.1 degree V06, 2019. URL: https://disc.gsfc.nasa.gov/datasets/GPM_3IMERGHHL_06/summary.
  * Steinwand et al. (1995) D. Steinwand, J. Hutchinson, J. Snyder,  Map projections for global and continental data sets and an analysis of pixel distortion caused by reprojection. photogramm. eng,  Rem. S 61 (1995) 1–487.
  * Ayoobi and Tangestani (2017) I. Ayoobi, M. H. Tangestani,  Evaluating the effect of spatial subsetting on subpixel unmixing methodology applied to aster over a hydrothermally altered terrain,  International journal of applied earth observation and geoinformation 62 (2017) 1–7.

  *[TC]: Tropical Cyclones
  *[IMERG]: Integrated Multi-satellitE Retrievals for GPM
  *[TCs]: Tropical Cyclones
  *[NIO]: North Indian Ocean
  *[BoB]: Bay of Bengal
  *[IMD]: India Meteorological Department
  *[JTWC]: Joint Typhoon Warning Centre
  *[NWP]: Numerical Weather Prediction
  *[AMSU]: Advanced Microwave Sounding Unit
  *[NASA]: National Aeronautics and Space Administration
  *[BMTPC]: Building Materials and Technology Promotion Council
  *[NDMA]: National Disaster Management Authority
  *[GPM]: Global Precipitation Measurement
  *[LEO]: Low-Earth-Orbit
  *[GEO]: Geosynchronous-Earth-Orbit
  *[IR]: Infrared
  *[ESCS]: Extremely Severe Cyclonic Storm
  *[CS]: Cyclonic Storm
  *[MSW]: Maximum Sustained Wind
  *[D]: Depression
  *[DOA]: Date of Acquisition
  *[JAXA]: Japan Aerospace Exploration Agency
  *[ROI]: Region of Interest
  *[CCA]: Connected Component Analysis


#  The energy, temperature, and heat capacity in discrete classical dynamics

Søren Toxvaerd  Department of Science and Environment, Roskilde University,
Post box 260, DK-4000 Roskilde, Denmark<EMAIL_ADDRESS>
###### Abstract

Simulations of objects with classical dynamics are in fact a particular
version of discrete dynamics since almost all the classical dynamics
simulations in natural science are performed with the use of the simple
”Leapfrog” or ”Verlet” algorithm. It was, however, Newton who in $Principia$,
$Proposition$ $I$ in 1687 first formulated the discrete algorithm, which much
later in 1967 was rederived by L. Verlet. Verlet also formulated a first-order
approximation for the velocity $\textbf{v}(t)$ at time $t$, which has been
used in simulations since then. The approximated expressions for
$\textbf{v}(t)$ and the kinetic energy lead to severe errors in the
thermodynamics at high densities, temperatures, strong repulsive forces, or
for large discrete time increments used in the discrete ”Molecular Dynamics”
(MD) simulations. Here we derive the exact expressions for the discrete
dynamics, and show by simulations of a Lennard-Jones system, that these
expressions now result in equality between temperatures determined from the
kinetic energies, and the corresponding configurational temperatures
determined from the expression of Landau and Lifshitz, derived from the
forces.

10.1103/PhysRevE.00.005300

##  I Introduction

Simulations of objects with classical dynamics are in fact a particular
version of discrete dynamics since almost all the classical dynamics
simulations in natural science are performed with the use of a simple
algorithm, first formulated by Newton Newton1687 ; Toxvaerd2020 . A new
position in Newton’s discrete dynamics, $\textbf{r}_{i}(t+\delta t)$, at time
$t+\delta t$ of an object $i$ with the mass $m_{i}$ is determined by the force
$\textbf{f}_{i}(t)$ acting on the object at the discrete position
$\textbf{r}_{i}(t)$ at time $t$ together with the position
$\textbf{r}_{i}(t-\delta t)$ at $t-\delta t$ as

$m_{i}\frac{\textbf{r}_{i}(t+\delta t)-\textbf{r}_{i}(t)}{\delta
t}=m_{i}\frac{\textbf{r}_{i}(t)-\textbf{r}_{i}(t-\delta t)}{\delta t}+\delta
t\textbf{f}_{i}(t),$ (1)

where the velocities $\textbf{v}_{i}(t+\delta t/2)=(\textbf{r}_{i}(t+\delta
t)-\textbf{r}_{i}(t))/\delta t$ and $\textbf{v}_{i}(t-\delta
t/2)=(\textbf{r}_{i}(t)-\textbf{r}_{i}(t-\delta t))/\delta t$ and
corresponding momenta are constant in the time intervals in between the
discrete positions. Newton begins Principia by postulating Eq. (1) in
Proposition I, and he obtained his second law as the limit $lim_{\delta
t\rightarrow 0}$ of the equation. For Newton’s derivation of his second law
and his discrete algorithm in Proposition I see Ref. 3, Chap. 2 Toxvaerd2023 .

The algorithm, Eq. (1), is usually presented as the Leapfrog algorithm

$\displaystyle\textbf{v}_{i}(t+\delta t/2)=\textbf{v}_{i}(t-\delta t/2)+\delta
t/m_{i}\textbf{f}_{i}(t)$ (2) $\displaystyle\textbf{r}_{i}(t+\delta
t)=\textbf{r}_{i}(t)+\delta t\textbf{v}_{i}(t+\delta t/2),$ (3)

where the positions are obtained from the discrete values of the velocities.
The rearrangement of Eq. (1) gives the Verlet algorithm Verlet1967

$\textbf{r}_{i}(t+\delta t)=2\textbf{r}_{i}(t)-\textbf{r}_{i}(t-\delta
t)+\delta t^{2}\textbf{f}_{i}(t)/m_{i}.$ (4)

The algorithm is used in almost all discrete ”Molecular Dynamics” (MD)
simulations and with Verlet’s expression for the velocity at the discrete time
$t$

$\textbf{v}_{0}(t)=\frac{\textbf{r}(t+\delta t)-\textbf{r}(t-\delta
t)}{2\delta t}.$ (5)

The approximation (5) for the velocity $\textbf{v}(t)$ is the first term in a
symmetrical Taylor expansion from the position r(t) on the classical analytic
trajectory, but the discrete dynamics trajectory is not analytic. Furthermore,
there are times for an analytic trajectory where a Taylor expansion is very
slowly converging, e.g. at particle collisions or for fast vibrations of atoms
in molecules with covalent bonds. The velocity is, however, as Newton noticed
in $Principia$ constant in times except at the discrete times where it is
exposed to a force impulsive (see next section), and it is given by Eq. (2).
Verlet’s inaccurate expression for the velocity results in a corresponding
inaccurate measure of the energy, temperature, and heat capacity.

The purpose of the article is to show that Newton’s discrete dynamics has the
same properties as his analytic Classical Mechanics, and to correct the
expressions for the energy, temperature, and heat capacity in MD simulations
Toxvaerd2023 .

In Sec. II we first present Newton and Verlet’s formulation of the discrete
classical dynamics and derive the energy invariance. The connections between
the ”Verlet” expressions for the energy, temperature, and heat capacity and
the corresponding exact Newtonian discrete dynamic expressions are derived in
Section III. In Section IV we determine the differences between the
traditional (Verlet) values of the temperature and heat capacities, and the
corresponding Newtonian values for a system of $N=2000$ particles with
Lennard-Jones (LJ) interactions. Section V summarizes the results.

##  II Discrete Newtonian dynamics

The classical discrete dynamics between $N$ spherically symmetrical objects
with masses $m^{N}=m_{1},m_{2},..m_{i},..,m_{N}$ and positions
r${}^{N}(t)=$r${}_{1}(t)$, r${}_{2}(t),..,$r${}_{i}(t),..$r${}_{N}(t)$ is
obtained by Eq. (1). Let the force $\textbf{f}_{i}(t)$ on object $i$ be a sum
of pairwise forces $\textbf{f}_{ij}(t)$ between pairs of objects $i$ and $j$

$\textbf{f}_{i}(t)=\sum_{j\neq i}^{N}\textbf{f}_{ij}(t).$ (6)

###  II.1 Newton’s formulation of discrete classical dynamics

Newton begins Principia by postulating Eq. (1) in Proposition I. The English
translation of Proposition I is

Of the Invention of Centripetal Forces.
PROPOSITION I. Theorem I.
The areas, which revolving bodies describe by radii drawn to an immovable
centre of force do lie in the same immovable planes, and are proportional to
the times in which they are described.
For suppose the time to be divided into equal parts, and in the first part of
time let the body by its innate force describe the right line AB. In the
second part of that time, the same would (by Law I.), if not hindered, proceed
directly to c, along the line Bc equal to AB; so that the radii AS, BS, cS,
drawn to the centre, the equal areas ASB, BSc, would be described. But when
the body is arrived at B, suppose that a centripetal force acts at once with a
great impulse, and, turning aside the body from the right line Bc, compels it
afterwards to continue its motion along the right line BC. Draw cC parallel to
BS meeting BC in C; and at the end of the second part of the time, the body
(by Cor. I of Laws) will be found in C, in the same plane with the triangle
ASB. Join SC, and, because SB and Cs are parallel, the triangle SBC will be
equal to the triangle SBc, and therefore also to the triangle SAB. By the like
argument, if the centripetal force acts successively in C, D, E, & c., and
makes the body, in each single particle of time, to describe the right lines
CD, DE, EF, & c., they will all lie in the same plane; and the triangle SCD
will be equal to the triangle SBC, and SDE to SCD, and SEF to SDE. And
therefore, in equal times, equal areas are described in on immovable plane:
and, by composition, any sums SADS, SAFS, of those areas, are one to the other
as the times in which they are described. Now let the number of those
triangles be augmented; and their breadth diminished in infinitum; and (by
Cor. 4, Lem III) their ultimate perimeter ADF will be a curve line: and
therefore the centripetal force, by which the body is perpetually drawn back
from the tangent of this curve, will act continually; and any described areas
SADS, SAFS, which are always proportional to the times of description, will,
in this case also, be proportional to those times. Q. E. D.

Proposition I in Principia is illustrated with a figure (Fig. 1). The
assumption ..suppose that a centripetal force acts at once with a great
impulse.. in Proposition I is highlighted here. It is the central part in
$Principia$. Newton’s second law is derived from his discrete dynamics, where
the time, forces, and the change in momenta are quantized and the velocities
are constant in between the discrete times. Newton’s discrete dynamics in
Proposition I is in fact Classical Quantum Dynamics Toxvaerd2023 .

Nobel Laureate T. D. Lee was perhaps the first to suggest that difference
equations are to be preferred in the foundation of dynamics, and with a
classical difference equation as the classical limit path of Feynman’s quantum
paths Lee1987 . This classical limit path is Newton’s discrete dynamics in
Proposition I. Our classical world is of course analytic. The difference
between analytic and discrete dynamics is proportional to the square of
$\delta t=$ Planck time ($\approx 10^{-43}$ sec), and negligible. But one can
say that Newtonian discrete dynamics simulations with time increments $\delta
t\approx 10^{-14}$ sec, used in many MD simulations is in the same
”universality class” as the analytic dynamics, and exact in the same manner as
classical analytic dynamics. If so it places Newtonian discrete dynamics as a
remarkably useful scientific tool.

Figure 1: Newton’s figure at $Proposition$ $I$ in Principia, with the
formulations of the discrete dynamics. The discrete positions are A:
$\bf{r}$$(t-\delta t)$; B: $\textbf{r}(t)$; C: $\textbf{r}(t+\delta t)$, etc..
The deviation cC from the straight line ABc (Newton’s first law) is caused by
a force from the position S at time $t$.

###  II.2 The Verlet algorithm

Loup Verlet (1931-2019) published in 1967 the article  Computer ”Experiments”
on Classical Fluids. I. Themodynamical Properties of Lennard-Jones Molecules
Verlet1967 , where Newton’s discrete algorithm with the formulation by Eq (4)
was introduced. The algorithm was derived by a forward and backward Taylor
expansion Toxvaerdappendix

$\displaystyle\textbf{r}(t+\delta t)=\textbf{r}(t)+\delta
t\frac{\partial\textbf{r}(t)}{\partial t}+\frac{1}{2}\delta
t^{2}\frac{\partial^{2}\textbf{r}(t)}{\partial t^{2}}+\frac{1}{6}\delta
t^{3}\frac{\partial^{3}\textbf{r}(t)}{\partial t^{3}}+..$
$\displaystyle\textbf{r}(t-\delta t)=\textbf{r}(t)-\delta
t\frac{\partial\textbf{r}(t)}{\partial t}+\frac{1}{2}\delta
t^{2}\frac{\partial^{2}\textbf{r}(t)}{\partial t^{2}}-\frac{1}{6}\delta
t^{3}\frac{\partial^{3}\textbf{r}(t)}{\partial t^{3}}+..,$ (7)

and the algorithm was obtained from the sum $\textbf{r}(t+\delta
t)+\textbf{r}(t-\delta t)$ and $\delta t^{2}\partial^{2}\textbf{r}(t)/\partial
t^{2}=\frac{\delta t^{2}}{m}\textbf{f}(t)$. All the odd terms in the sum
cancel, and the Verlet algorithm was formulated as a four-order time symmetric
predictor of the positions at the analytic trajectories. The velocity is
correspondingly derived from the difference $\textbf{r}(t+\delta
t)-\textbf{r}(t-\delta t)$, all even terms cancel, an it is

$\textbf{v}(t)=\frac{\textbf{r}(t+\delta t)-\textbf{r}(t-\delta t)}{2\delta
t}+\frac{1}{6}\frac{\delta t^{2}}{m}\textbf{f}^{\prime}(t)+..$ (8)

The scientific community and Verlet were first much later aware, that it
actually was Newton who first published the geometric formulation of the
algorithm in $Proposition$ $I$ Nauenberg2018a .

###  II.3 The energy invariance in discrete Newtonian dynamics

Newton’s algorithm is a symmetric time centered difference whereby the
dynamics is time reversible and symplectic. The conservation of momentum and
angular momentum for a conservative system follows directly from Newton’s
third law for the conservative system with the forces
$\textbf{f}_{ij}(t)=-\textbf{f}_{ji}(t)$ between objects $i$ and $j$, but the
energy invariance is not so obvious.

The energy in a system with analytic dynamics is the sum of potential energy
$U(\textbf{r}^{N}(t))$ and kinetic energy $K(t)$, and it is time invariant for
a conservative system. The kinetic energy in the discrete dynamics at time $t$
is, however, not well-defined since the velocities change at $t$.
Traditionally one uses Verlet’s first-order expression for the velocity at
time $t$, Eq. (5), obtained by the symmetric Taylor expansion, and

$\displaystyle K_{0}(t)=\sum_{i}^{N}\frac{1}{2}m_{i}\textbf{v}_{0,i}(t)^{2}$
(9) $\displaystyle E_{0}(t)=U(\textbf{r}^{N}(t))+K_{0}(t)$ (10)

The energy $E_{0}$ obtained by using the approximation (5) and with
$K(t)=K_{0}(t)$ for the kinetic energy fluctuates with time, although it is
constant when averaged over long time intervals.

The velocities are, however, constant in between the discrete times in
Newton’s discrete dynamics, and the energy invariance can be obtained by
considering the energy in the time interval $[t-\delta t/2,t+\delta t/2]$, and
dividing the interval into two sub-intervals $[t-\delta t/2,t]$ and
$[t,t+\delta t/2]$. The energy invariance, $E_{\textrm{D}}$ in Newton’s
discrete dynamics (D) can then be obtained by considering the change in
kinetic energy, $\delta K_{\textrm{D}},$ the change in potential energy,
$\delta U_{\textrm{D}},$ and the work, $W_{\textrm{D}}$ done by the forces in
the time interval $[t-\delta t/2,t+\delta t/2].$

The loss in potential energy, $-\delta U_{\textrm{D}}$ is defined as the work
done by the forces for a change of positions Goldstein . The work,
$W_{\textrm{D}}$ done in the time interval by the discrete dynamics from the
position $(\textbf{r}_{i}(t)+(\textbf{r}_{i}(t-\delta t))/2$ at $t-\delta t/2$
to the position $(\textbf{r}_{i}(t+\delta t)+\textbf{r}_{i}(t))/2$ at
$t+\delta t/2$, and with the change $\delta\textbf{r}_{i}$ of the position
$\delta\textbf{r}_{i}=\textbf{(}\textbf{r}_{i}(t+\delta
t)-\textbf{r}_{i}(t-\delta t))/2$ is Toxvaerd2023

$-\delta
U_{\textrm{D}}=W_{\textrm{D}}=\sum_{i}^{N}\textbf{f}_{i}(t)(\textbf{r}_{i}(t+\delta
t)-\textbf{r}_{i}(t-\delta t))/2.$ (11)

By rewriting Eq. (4) to

$\textbf{r}_{i}(t+\delta t)-\textbf{r}_{i}(t-\delta
t)=2(\textbf{r}_{i}(t)-\textbf{r}_{i}(t-\delta t))+\frac{\delta
t^{2}}{m_{i}}\textbf{f}_{i}(t),$ (12)

and inserting in Eq. (11) one obtains an expression for the total work in the
time interval

$-\delta
U_{\textrm{D}}=W_{\textrm{D}}=\sum_{i}^{N}[(\textbf{r}_{i}(t)-\textbf{r}_{i}(t-\delta
t))\textbf{f}_{i}(t)+\frac{\delta t^{2}}{2m_{i}}\textbf{f}_{i}(t)^{2}].$ (13)

The mean kinetic energy $K_{\textrm{D}}$ of the discrete dynamics in the time
interval $[t-\delta t/2,t+\delta t/2]$ is

$\displaystyle K_{\textrm{D}}=$
$\displaystyle\frac{1}{2}\sum_{i}^{N}\frac{1}{2}m_{i}[\frac{\textbf{(}\textbf{r}_{i}(t+\delta
t/2)-\textbf{r}_{i}(t))^{2}}{\delta(t/2)^{2}}+\frac{\textbf{(}\textbf{r}_{i}(t)-\textbf{r}_{i}(t-\delta
t/2))^{2}}{\delta(t/2)^{2}}]$
$\displaystyle=\frac{1}{2}\sum_{i}^{N}\frac{1}{2}m_{i}[\frac{\textbf{(}\textbf{r}_{i}(t+\delta
t)-\textbf{r}_{i}(t))^{2}}{\delta
t^{2}}+\frac{\textbf{(}\textbf{r}_{i}(t)-\textbf{r}_{i}(t-\delta
t))^{2}}{\delta t^{2}}],$ (14)

with the change

$\displaystyle\delta K_{\textrm{D}}=$
$\displaystyle\sum_{i}^{N}\frac{1}{2}m_{i}[\frac{\textbf{(}\textbf{r}_{i}(t+\delta
t)-\textbf{r}_{i}(t))^{2}}{\delta
t^{2}}-\frac{\textbf{(}\textbf{r}_{i}(t)-\textbf{r}_{i}(t-\delta
t))^{2}}{\delta t^{2}}].$ (15)

By rewriting Eq. (4) to

$\textbf{r}_{i}(t+\delta
t)-\textbf{r}_{i}(t)=\textbf{r}_{i}(t)-\textbf{r}_{i}(t-\delta t)+\frac{\delta
t^{2}}{m_{i}}\textbf{f}_{i}(t)$ (16)

and inserting the squared expression for $\textbf{r}_{i}(t+\delta
t)-\textbf{r}_{i}(t)$ in Eq. (15), the change in kinetic energy is

$\delta K_{\textrm{D}}=\sum_{i}^{N}[(\textbf{r}_{i}(t)-\textbf{r}_{i}(t-\delta
t))\textbf{f}_{i}(t)+\frac{\delta t^{2}}{2m_{i}}\textbf{f}_{i}(t)^{2}].$ (17)

Figure 2:  The relative mean difference $<(T_{D}-T_{0})/T_{D}>$ between
$T_{D}$ and $T_{0}$ as a function of the density for a LJ fluid system at
$T_{D}=1$.

The energy invariance from the a discrete change of time from $t-\delta t/2$
to $t+\delta t/2$ in Newton’s discrete dynamics is expressed by Eqn. (13), and
Eq. (17) as Toxvaerd2023

$\delta E_{\textrm{D}}=\delta U_{\textrm{D}}+\delta K_{\textrm{D}}=0.$ (18)

##  III The kinetic energy, temperature, and heat capacity in a Newtonian
discrete conservative system

The constant velocities in the time intervals in between a force impulse at
time $t$ are related. It is easy to derive the relation Toxvaerd2013

$\textbf{v}_{0,i}(t)^{2}=\frac{1}{2}\textbf{v}_{i}(t+\delta
t/2)^{2}+\frac{1}{2}\textbf{v}_{i}(t-\delta t/2)^{2}-\frac{1}{4}(\frac{\delta
t}{m_{i}}\textbf{f}_{i}(t))^{2}$ (19)

from Eq. (4) between the square of the first term in Verlet’s Taylor expansion
$\textbf{v}_{0,i}(t)^{2}$ , Eq. (5), for the velocity at time $t$ and Newton’s
exact expression for the square of the constant velocities in the time
intervals $[t-\delta t,t]$ and $[t,t+\delta t]$. The corresponding expression
for the kinetic energy, $K_{0}(t)$, and the traditional value for the
temperature used in MD simulations

$T_{0}=\frac{<2K_{0}(t)>}{N_{f}}$ (20)

used in MD simulations for a system with $N_{f}$ degrees of freedom is less
than the mean kinetic energy, $K_{D}$, and the temperature $T_{\textrm{D}}$.
In the discrete time interval $[t-\delta t/2,t+\delta t/2]$ the relation is

$K_{0}(t)=K_{D}(t)-\sum_{i}^{N}\frac{1}{8}\frac{\delta
t^{2}}{m_{i}}\textbf{f}_{i}(t)^{2},$ (21)

and the systematic difference $T_{D}-T_{0}$ in the time interval by using
Verlet’s first-order approximation is

$T_{D}(t)-T_{0}(t)=\sum_{i}^{N}\frac{1}{4}\frac{\delta
t^{2}}{m_{i}}\textbf{f}_{i}(t)^{2}/N_{f}.$ (22)

Newton’s discrete dynamics depends purely on the positions and the forces. The
momenta are not dynamic variables, but only ”book-keeping” expressions in the
dynamics Toxvaerd2023 . Therefore, the configurational temperature

$T_{\textrm{Conf}}(t)=<\nabla^{2}U(\textbf{r}^{N})>=\frac{\sum_{i}^{N}\textbf{f}_{i}(t)^{2}}{-\sum_{i}^{N}\nabla\cdot\textbf{f}_{i}(t)}.$
(23)

could be a more relevant expression for the temperature since it depends
purely on the forces. Eq. (23) is derived from the average of the Laplacian of
the potential energy $\nabla^{2}U(\textbf{r}^{N})$, and is obtained from
canonical averaging in the configurational phase space Landau ; Rugh1997 . The
value of $T_{\textrm{Conf}}(t)$ fluctuates with time, but its mean value
obtained from long simulations agrees with the corresponding temperatures,
$T_{D}$ obtained from the kinetic energy.

The heat capacity $C_{V}=3/2R+C_{V}^{i}$ consists of two terms, the term from
the kinetic energy and the term $C_{V}^{i}$ from the interactions. The heat
capacity can be determined either from numerical differentiation of the energy
$U(T,\rho)$ or from the mean square fluctuation $<K(t)^{2}>-\bar{K}^{2}$ of
the kinetic energy Lebowitz1967 . The term
$K_{D}(t)-K_{0}(t)=\sum_{i}^{N}\frac{\delta t^{2}}{8m}\textbf{f}_{i}^{2}$ in
the kinetic energy $K_{0}(t)$ is important at high densities with hard
particle collisions, and it affects $C_{V}^{i}$. The ratio

$\Delta
C_{V}^{i}=\frac{<(K_{D}(t)-\bar{K}_{D})^{2}>-<(K_{0}(t)-\bar{K}_{0})^{2}>}{<(K_{D}(t)-\bar{K}_{D})^{2}>}$
(24)

is a measure of the relative error in $C_{V}^{i}$ by using Eq. (9) for the
kinetic energy.

Figure 3:  Temperature fluctuations in a $NVE$ system of $N=2000$ LJ particles
at the density $\rho=1.40$, temperature $T_{D}=1$, and with $\delta t=$0.010.
With red is $T_{D}(t)$; with blue is the traditional expression $T_{0}(t)$ for
the temperature, and with thin green is the configurational temperature
$T_{Conf}(t)$.

##  IV Simulations.

Verlet used Lennard-Jones (LJ) forces to simulate argon and later his force
field (FF) has been a standard used in many simulations, e.g., simulations of
organic molecules to describe the atom-atom FF between the atoms in the
molecules together with stronger FF, e.g. for covalent bonds. The errors
$T_{D}-T_{0}$, $E_{D}-E_{0}$, and in the heat capacities by using Verlet’s
first order approximation for the velocity are biggest for strong force
fields, and especially at collisions where the FF is strongly repulsive. So
the errors for systems with strong intramolecular FF and by using Eq. (20) are
bigger than for an LJ system.

LJ systems with $N=2000$ LJ particles were simulated at $T_{D}=1.00$ for
different densities, temperatures, strength of the repulsive forces, and time
increments. The difference between $T_{D}$ and $T_{0}$ at $T_{D}=1.00$ as a
function of density for $10^{6}$ time steps are shown in Fig. 2 for different
time increments $\delta t$ and repulsive forces given by the exponent $n$ in
the $n-6$ LJ potential Heyes1998  . The differences increase with density,
temperature, the strength of the repulsive forces, or the discrete time step
$\delta t$, but they are relatively small at state points with low densities.
But the differences are significant for state points with high densities,
systems with stronger forces, and for large time increments $\delta t$.

Most MD simulations are with a thermostat ($NVT$ and $NPT$ ensemble MD) and
with the thermostat temperature $T(thermostat)=T_{0}$. We have performed $NVT$
simulations with a Nosé-Hover thermostat Hoover ; Toxvaerd2024  and with
$T(thermostat)=T_{0}$ as well as $T(thermostat)=T_{D}$. The results from the
$NVT$ simulations agree with the $NVE$ data. So the $NVE$ results in Fig. 2
are also valid for the other ensemble simulations.

The fluctuations of the temperatures $T_{D}(t),T_{0}(t)$, and $T_{Conf}(t)$ in
the time interval $t\in[4000\delta t,4300\delta t]$ for the MD are shown in
Fig. 3. The temperatures are obtained for $N=2000$ LJ particles at
$(T_{D},\rho)=(1.00,1.40)$, and with $\delta t=0.010$. The red curve is for
three hundred time steps for $T_{D}(t)$, and the green curve is $T_{Conf}(t)$.
The mean temperatures $T_{D}$ and $T_{Conf}$ from $10^{6}$ time steps are
approximately equal, $T_{D}=0.999\pm 0.011$ and $T_{Conf}=1.008\pm 0.020$, and
higher than the mean $T_{0}=0.954\pm 0.011$ for $T_{0}(t)$ shown in blue. The
temperatures in Fig. 3 show, that the temperature $T_{D}$ in Newtonian
discrete dynamics is consistent with the otherwise obtained configurational
temperature $T_{Conf}$, but is not consistent with $T_{0}$.

Figure 4:  The relative difference in the time interval $t\in[4000\delta
t,4300\delta t]$ between the fluctuations in the square of the excess kinetic
energies $\Delta C_{V}^{i}$ for the LJ system at $(T,\rho)=(1.00,1.40)$ and
with $\delta t=0.010$. The mean of the difference of the square excess kinetic
energies (Eq. 24) is -0.12. The inset shows the two excess square kinetic
energies per particle: with red is $K_{D}(t)^{2}-\bar{K}_{D}^{2}$ and with
blue is $K_{0}(t)^{2}-\bar{K}_{0}^{2}$. Figure 5:  The relative difference
between the fluctuations in the square of the excess kinetic energies $\Delta
C_{V}^{i}$ for the LJ system at $T=1.00$ for different densities and time
increments. The corresponding relative differences in temperatures are shown
in Fig. 2.

The relative difference between fluctuations in the square of the excess
kinetic energies $\Delta C_{V}^{i}$ for a few hundred time steps is shown in
Fig. 4 for the LJ system, and with the corresponding fluctuations
$K_{D}(t)^{2}-\bar{K}_{D}^{2}$ and $K_{0}(t)^{2}-\bar{K}_{0}^{2}$ in the inset
in the figure. The fluctuations are for the LJ system at
$(T_{D},\rho)=(1.00,1.40)$ with the corresponding temperature fluctuations
shown in Fig. 3. At first glance, the fluctuations are synchronous (Fig. 3 and
the inset in Fig. 4), but a closer examination reveals that there are
significant differences and of the order of twelve percent for the state point
with relatively high density $\rho=1.40$ and for the time increment $\delta
t=0.010$.

The relative difference in $C_{V}^{i}$ from $10^{6}$ time steps for the system
at $T_{D}=1.00$ and at different densities and time increments is shown in
Fig. 5. (The corresponding relative temperature differences are shown in Fig.
2.) The differences between $C_{V}^{i}$ obtained by fluctuations in $K_{D}(t)$
and $K_{0}(t)$, respectively, increases with increasing density, temperature,
strength of the repulsive forces, or time increment.

##  V  Conclusion.

Newton published his discrete algorithm in Principia Proposition I , and the
algorithm is now used everywhere in natural science disciplines ranging from
simulations of planetary systems to simulations of atoms, molecules, and
systems with complex proteins. L. Verlet rederived in 1967 the algorithm from
a forward and backward Taylor expansion Toxvaerdappendix  and performed his MD
simulations with Eq. (5) for the velocity and Eq. (9) for the kinetic energy.
Discrete dynamics with Newton’s algorithm has been further developed and is
now a standard tool used in natural science, and MD was honored with a Nobel
Prize in 2013.

Newton’s discrete dynamics has the same properties as his analytic Classical
Mechanics, it is time reversible, symplectic, and has the same invariances for
a conservative system Toxvaerd2023 . Furthermore, there exists a shadow
Hamiltonian where the positions of the discrete dynamics are located on the
analytical trajectories for the shadow Hamiltonian Toxvaerd2023 ; Toxvaerd1994
; Toxvaerd2012 . This means that there is no qualitative difference between
the analytical and the discrete dynamics, and MD with the Newton-Verlet
algorithm is the exact generation of positions for the discrete dynamics.

Verlet’s approximative expression for the velocity at time $t$, Eq. (5), which
is used in MD simulations does not affect the discrete dynamics and the
discrete positions, but it is an unnecessary approximation. It is the first
term in his Taylor expansion, and the corresponding expression for the kinetic
energy, Eq. (9) ought to be replaced with the exact expression, Eq. (14).
Doing so one achieves an agreement between the kinetic and configurational
temperature (Fig. 3). Potential energy fluctuations $\delta U$ should be
obtained by Eq. (11) or (13), and velocities and their time correlations by
the (constant) velocities in the respective time intervals. Most MD
simulations are canonical $NVT$ ensemble or $NPT$ ensemble simulations with
the thermostat temperature $T(thermostat)=T_{0}$, and not only are the
$NVT_{0}$ and $NPT_{0}$ temperatures wrong so are also the values of
$C_{V}^{i}$. The errors caused by using Eq.(5) for the velocity and Eq. (9)
for the kinetic energy might only be a few percent for systems with relatively
weak forces. But the errors increase according to Eq. (22) with the strength
of the forces and the time increment $\delta t$ used in the simulations, and
can lead to severe errors for systems with strong forces such as full-atomic
models of e.g. biomolecules with fast intramolecular vibrations.

The biggest challenge with an MD simulation of a real system, however, is to
formulate a correct force field, and also for this reason it is important to
use the correct expression for the kinetic energy, whereby one avoids
introducing a systematic error in the calculated energies. The public software
used for MD simulations should be corrected.

## References

  * (1) I. Newton, PHILOSOPHIÆ NATURALIS PRINCIPIA MATHEMATICA. LONDINI, Anno MDCLXXXVII.
  * (2) S. Toxvaerd, Eur. Phys. J. Plus 135, 267 (2020).
  * (3) S. Toxvaerd, Comprehensive Computational Chemistry 3, 329 (Elsevier, 2023).
  * (4) L. Verlet, Phys. Rev. 159, 98 (1967).
  * (5) T. D. Lee, J. Stat. Phys. 46, 843 (1987).
  * (6) See Appendix in Toxvaerd2023 .
  * (7) M. Nauenberg, Am. J. Phys. 86, 765 see Reference No. 7. (2018).
  * (8) H. Goldstein, Classical Mechanics,( Addison-Wesley Press Second Ed. 1980), Chap. 1.
  * (9) S. Toxvaerd. J. Chem. Phys. 139, 224106 (2013).
  * (10) L. D. Landau and E. M. Lifshitz, Statistical Physics (Pergamon, Oxford, 1958), Eq.(33.14).
  * (11) H. H. Rugh, Phys. Rev. Lett. 78, 772 (1997).
  * (12) J. L. Lebowitz, J. K. Percus and L. Verlet Phys. Rev., 153, 250 (1967).
  * (13) D. Heyes and J. Powels, Mol. Phys. 95, 259 (1998).
  * (14) W. G. Hoover, Phys. Rev. A, 31, 1695 (1985).
  * (15) S. Toxvaerd, submitted to Mol. Phys.
  * (16) S. Toxvaerd Phys. Rev. E, 50, 2271 (1994).
  * (17) S. Toxvaerd, O. J. Heilmann, J. C. Dyre, J. Chem. Phys. 136, 224106 (20012).