Vol.0 (20xx) No.0, 000–000

11institutetext: Center for Astrophscis, School of Physics and Materials
Science, Guangzhou University, Guangzhou 51006, China<EMAIL_ADDRESS>
22institutetext: Great Bay Center, National Astronomical Data Center,
Guangzhou, Guangdong, China, 510006
33institutetext: Department of Physics and Electronics, Rhodes University, PO
Box 94, Makhanda, 6140, South Africa
44institutetext: Shanghai Astronomical Observatory, Chinese Academy of
Sciences, Shanghai, China, 200030 Received 2021 November 22; accepted 2022
February 16

# Performance evaluation of baseline-dependent averaging based on full-scale
SKA1-LOW simulation

Qing-Wen Deng 1122 Feng Wang 1122 Hui Deng 1122 Ying Mei 1122 Jing Li 1122
Oleg Smirnov 33 Shao-Guang Guo 44

###### Abstract

The Square Kilometre Array (SKA) is the largest radio interferometer under
construction in the world. Its immense amount of visibility data poses a
considerable challenge to the subsequent processing by the science data
processor (SDP). Baseline dependent averaging (BDA), which reduces the amount
of visibility data based on the baseline distribution of the radio
interferometer, has become a focus of SKA SDP development. This paper
developed and implemented a full-featured BDA module based on Radio Astronomy
Simulation, Calibration and Imaging Library (RASCIL). Simulated observations
were then performed with RASCIL based on a full-scale SKA1-LOW configuration.
The performance of the BDA was systematically investigated and evaluated based
on the simulated data. The experimental results presented that the amount of
visibility data is reduced by about 50% to 85% for different time intervals
($\Delta t_{max}$). In addition, different $\Delta t_{max}$ have a significant
effect on the imaging quality. The smaller the $\Delta t_{max}$, the smaller
the degradation of the imaging quality.

###### keywords:

instrumentation: interferometers — methods: analytical — methods: numerical —
techniques: interferometric

##  1 Introduction

The Square Kilometre Array (SKA) (Braun 1996) is an ongoing international
project to build the world’s largest radio interferometric telescope with more
than one square kilometer potential collection area. With detailed design and
preparation well underway, the SKA represents a giant leap forward in
engineering to deliver a unique instrument.

The $uv$ distribution of a radio interferometer generally has a dense center
and a sparse edge. With the rotation of the Earth, each sampling point of the
radio interferometer draws an arc-shaped trajectory around the phase center on
the $uv$ plane. Short baselines have denser data than long baselines for the
same track length, and the difference can be even greater for a larger array.
Decorrelation can be avoided on the longer baselines when more samples are
averaged at the center than at the outer edges. At the same time, data
compression can be carried out on shorter baselines.

The baseline dependent averaging (BDA) was presented by Cotton (1986, 1999) to
reduce the visibility data volumes, which has been used by the MWA (Mitchell
et al. 2008) and is also used to shape the field of interest (Atemkeng et al.
2018). However, averaging visibilities over time and frequency will cause
image distortion, also called the smearing effect. The bandwidth smearing is
manifested as a position-dependent and radial convolution effect in the image
field. Time smearing is similar to bandwidth smearing but more complicated,
described as a loss in amplitude (Cotton 1986, 1999). Bridle & Schwab (1999)
conducted a mathematical analysis of these two smearing effects and found that
they cannot be effectively corrected by calibration or self-calibration
methods. It is recommended to design a comprehensive observation strategy to
reduce the impact to an acceptable level. Therefore, the short integration
time and small channel width are necessary for the long baselines for a radio
interferometer to suppress the smearing effect. In contrast, the resolution
requirement for time and frequency is relatively low on the short baselines.

The BDA does not change the channel bandwidth and integration time for long
baseline visibility data, but averaging of short baseline visibility data
corresponds to an increase in integration time and channel bandwidth.
Wijnholds et al. (2018) obtained the Cramer-Rao Bound of averaged visibilities
by estimating the number of the raw visibilities and compared it with the
covariance obtained by the error transfer formula. It is proven that BDA will
not cause other effects except for the approximately obtainable decorrelation
loss. Salvini & Wijnholds (2017) proposed the Compress-Expand-Compress method
to expand the visibilities to the required time resolution for calibration
after the first compression, and then perform the second compression, and
finally achieve a high compression ratio of the visibilities in time. However,
few previous literature presented the quantitative analysis of BDA on the
final image quality and storage costs.

As the SKA enters its construction phase and the SKA Science Regional Centers
(SRCs) construction white paper has been released, it becomes imperative to
research the BDA technique further and analyze its usability for the SKA1
scale. We wish to analyze and discuss this study systematically: 1. How much
space would be saved by using BDA technology for SKA1-LOW observations? 2. Is
there a significant degradation of the final dirty image with the BDA?

In the rest of this study, we first introduce the BDA algorithm and its
implementation. We then simulate full-scale SKA1-LOW observations and
investigate the BDA performance in Section 3. The discussions are described in
Section 4. The conclusions and future work are presented in the last section.

##  2 The algorithm and implementation of BDA

###  2.1 The BDA Algorithm

For a radio interferometer, a visibility function is obtained by correlating
the signal collected by two antennas of each baseline with the same time
interval $\delta t$ and frequency sampling interval $\delta f$. According to
the mathematical definition of BDA (Wijnholds et al. 2018), we can average the
raw visibilities and thus obtain the averaged visibilities. In data
processing, from $P$ receiving antennas $P^{2}$ correlations are assumed to be
collected in $K$ short-term integrations, either over time, frequency, or
both. The raw visibility data vector can be defined as

$\mathbf{r}=\left[r_{1},\cdots,r_{K}\right]^{T}\in\mathbb{C}^{KP^{2}\times 1}$
(1)

where $\mathbb{C}$ denotes the complex matrix. And the averaging process can
be described as

$\mathbf{r}_{\mathrm{ave}}=\mathbf{W}\mathbf{I}_{\mathrm{s}}\mathbf{r}\in\mathbb{C}^{M\times
1}$ (2)

where
$\mathbf{W}=\left(\mathbf{I}_{\mathrm{s}}\mathbf{I}_{\mathrm{s}}^{T}\right)^{-1}$
is a diagonal weighting matrix. $\mathbf{I}_{\mathrm{s}}\in\mathbb{R}^{M\times
KP^{2}}$ is a sparse matrix that determines which samples will be averaged
together, where $\mathbb{R}$ denotes the real matrix and $M$ is the length of
the average visibility data vector.

Suppose we ignore the correlation effects and assume that the values of the
visibility data averaged together are the same. In that case, the raw
visibilities can be obtained approximately from the averaged visibilities by

$\mathbf{r}=\mathbf{I}_{\mathrm{s}}^{T}\mathbf{r}_{\mathrm{ave}}$ (3)

The selection matrix $\mathbf{I}_{\mathrm{s}}$ is related to the averaging
intervals of time and frequency in the BDA. For a selected baseline $D$, the
averaging intervals can be calculated by the rounding ratio of that baseline
to the longest baseline $B_{max}$ as

$t=\delta t\frac{B_{max}}{D},~{}~{}~{}~{}f=\delta f\frac{B_{max}}{D}$ (4)

We used the rounding ratio of $B_{max}$ to $D$ to determine whether a BDA
processing is needed. On partial long baselines, when $\frac{B_{max}}{D}=1$,
the data will not be averaged. While on the short baselines,
$\frac{B_{max}}{D}>1$, indicating that more sampling data can be averaged. In
addition, $t$ is limited by the calibration time scale determined by the
environment and instrument for the interferometer. A larger averaging scale
will also make the smearing effects more serious. Therefore, it is necessary
to set reasonable upper limits for $t$ and $f$ in the implementation,
respectively.

###  2.2 Implementation

We implemented a full-featured BDA module based on the Radio Astronomy
Simulation, Calibration and Imaging Library (RASCIL)
111https://gitlab.com/ska-telescope/external/rascil. The RASCIL is a pure
Python software package suite for radio interferometer calibration and imaging
algorithms, especially for SKA data processing. Since the public release of
RASCIL, it has been widely used in data processing for some radio interference
telescopes(e.g., Wei et al. (2021) and so on). We developed a BDA module based
on RASCIL, which has been released online 222https://github.com/astronomical-
data-processing/ska-bda.

The flowchart of the implementation is shown in Figure 1. Figure 2 shows an
example of averaging process on a baseline, where the data will be reduced
from the original $21\times 14$ to $5\times 3$, by assuming that both $t$ and
$f$ are 5 according to the Equation 4. In the case of the shortest baselines,
$t$ and $f$ are also equal to the upper limits of the time-frequency interval
used in the BDA, defined as $\Delta t_{max}$ and $\Delta f_{max}$. In the
averaging process, we first calculate the position of the raw visibilities
corresponding to the flattened averaged visibilities. We then average the
visibilities based on the positions and the number of visibilities.

Figure 1: The flow chart of the BDA implementation. Figure 2: The averaging of
visibility data at a baseline when applying BDA.

We used the BlockVisibility class defined in the RASCIL. The visibility data
was stored using a multi-dimensional array, with dimensions including
baseline, polarization, time and frequency. To meet the requirements of BDA
performance profiling, we implemented BDA by using three underlying packages,
i.e., pure Python, Pandas, and Numba, respectively.

In the pure Python implementation, we grouped the data for averaging based on
Numpy. To optimize the performance, we tried to use Pandas, put all the
parameter data into a table when preprocessing, and then performed group-by
operations to complete all the calculations.

We also used Numba (Lam et al. 2015) to speed up the function and further
improved the processing performance. Numba is an on-the-fly compiler that
translates a subset of the Python and Numpy code used in the function into
efficient machine code, which can effectively improve the speed of the
program.

##  3 Performance assessments For BDA

###  3.1 Observational configuration

To more accurately evaluate the performance of the BDA, we used RASCIL to
simulate single-channel and one polarization visibility data observed by the
full scale SKA1-LOW telescope. During simulation, we used all SKA1-LOW 512
stations and set up 12 minutes of observations, of which 6 minutes on each
side of the zenith. The integration time is set to 0.9 seconds as required by
the array structure. The observing frequency is 100 MHz with the channel
bandwidth 1 MHz, and the phase center in the observation points to the right
ascension $15^{\circ}$ and the declination $-45^{\circ}$. We finally obtained
800 temporal sampling points on each baseline. The $uv$ distribution of
simulated observation is shown in Figure 3.

Figure 3: The $uv$ distribution of the SKA1-LOW observed in a single channel
at 100 MHz for 12 minutes.

###  3.2 Observation Simulation

With the observational configuration described above, we simulated
observations of point and extend sources separately. We selected the
corresponding sources from the GaLactic and Extragalactic All-sky Murchison
Widefield Array (GLEAM) survey catalog (Hurley-Walker et al. 2017) with an
imaging size of 32,768 multiplied by 32,768 pixel. The other is the M31 image
that observed by Very Large Array. The image has a pixel size of $512\times
512$ and a resolution of 1 arc second. We used the transform.resize() function
in skimage (van der Walt et al. 2014) to scale the M31 image to the same pixel
scale as the image generated by GLEAM model for this study. It should be clear
here that such a magnification of the original image is only necessary for the
simulation of the observation. Two dirty images for the cases are shown in
Figure 4.

Figure 4: Dirty images of the raw visibilities observed from the GLEAM model
and the M31 model images.

###  3.3 Evaluation Results

We invoke the BDA module to perform visibilities compression, decompression,
and imaging processing. Since BDA processing at frequency series has similar
results at time series, we only simulate and analyze temporal BDA in this
study. Also, the effects on visibility data and dirty images are investigated
by setting different upper limits on time integration in the BDA.

####  3.3.1 The compression ratio

The compression ratio (CR) is calculated by comparing the data volumes between
the averaged and raw visibilities, as
$\mathrm{CR}=\frac{\mathbf{r}_{\mathrm{ave}}}{\mathbf{r}}\times 100\%$.

We set different upper limits ($\Delta t_{max}$), i.e., 1, 2, 4, 8, 12, 16,
32, 48, 64, 128 and 256, to evaluate the compression ratio of the BDA. $\Delta
t_{max}$ also means the maximum number of samples being averaged together on
the shortest baselines. These different upper limits lead to different
compression ratios over the baseline length range, and some of these variation
curves are shown in Figure 6. Since the shorter baselines have larger data
volumes, we want to average the amount of data over the shorter baselines as
much as possible. We can also find that the increase of $\Delta t_{max}$ only
further compresses the data on the shorter baselines, while the volume
proportion of these data is decreasing in the total.

The final result of the compression ratio with different $\Delta t_{max}$ is
shown in Figure 6. With the increasing value of $\Delta t_{max}$, the
compression ratio reduces quickly and then becomes slow. Finally, a larger
$\Delta t_{max}$ does not significantly improve the final compression ratio.
It changes very little after the compression ratio reaches 15%, where $\Delta
t_{max}=48$. To avoid the more severe errors that may arise from a bigger
$\Delta t_{max}$, it is worth considering using a $\Delta t_{max}$ less than
48 in subsequent studies.

Figure 5: The distribution of visibilities with baseline length, and the
effect curves of the compression ratio.

Figure 6: The trend chart of the compression ratio with the maximum number of
samples in averaging $\Delta t_{max}$.

####  3.3.2 Imaging Quality Evaluation

In addition to the compression ratio, the impact on the image quality after
applying BDA is a significant issue. We defined the visibility data of the
simulated observation as raw visibilities, the visibility data processed by
the BDA as the averaged visibilities, and the final decompressed visibility
data as recovered visibilities. We first decompressed the averaged
visibilities by using a method similar to linear interpolation and then used
the recovered visibilities for the subsequent imaging processing.

To exclude the possible effects of different deconvolution methods on the
imaging results, we used dirty images to analyze the imaging quality. Due to
the difference between the recovered visibilities and the raw visibilities,
the brightnesses in the dirty images are not the same. The deviations may be
positive or negative relative to the dirty image of the raw visibilities. For
convenience, the absolute value of the deviation is used here, and its
distribution with the brightness is shown in Figure 7, where $\Delta t_{max}$
of the recovered visibilities is equal to 48.

Figure 7:  The distribution of the absolute value of the residual with
brightness, corresponding to the dirty image results of the GLEAM (left) and
the M31 (right) as observed model images, and $\Delta t_{max}$ = 48.

The brightness distribution of the dirty images is mainly concentrated around
0, and the lower limit of deviation increases with brightness, but the upper
limit does not change excessively. At the same time, the maximum value of the
residuals is small.

We used a pure noise image as a sky model for simulated observations and
performed the same BDA processing. In generating this noisy model image, the
same image size as the previous model was used, filled only with Gaussian
noise with a mean of 0 and a standard deviation of 0.1. Figure 9 shows the
result of the dirty image when $\Delta t_{max}$ is equal to 48, where the
maximum error is 0.0698. Moreover, when $\Delta t_{max}$ is equal to 2, the
maximum error is 0.0162. During this experiment, we tried to reduce the
overall amplitude of the noisy model image by a certain ratio. The
corresponding change in the dirty image is that both deviation and brightness
values are reduced by the same ratio, while the contours of Figure 9 do not
change much.

The statistical distribution of the residuals shows a Gaussian-like
distribution in the dirty image results when displayed in logarithmic form.
The mean and standard deviation of the residuals in this form are given in
Table 1, where Case 1 refers to the results of dirty images for the GLEAM, and
Case 2 refers to the one for the M31.

The standard deviations of these two cases are approximately the same and do
not change significantly with $\Delta t_{max}$. This indicates that the change
in residual is more like an overall shift, while the mean is the distance of
the shift.

Table 1: Mean and standard deviation of the logarithmic residual in two cases. |  | Case 1 | Case 2
---|---|---|---
$\Delta t_{max}$ | CR | Mean | Std. | Max | Mean | Std. | Max
2 | 50.16% | -6.5449 | 0.5857 | -3.1033 | -7.6400 | 0.5309 | -5.4398
4 | 28.56% | -6.0538 | 0.5570 | -2.6526 | -6.6431 | 0.5292 | -4.4348
8 | 20.06% | -5.8964 | 0.5514 | -2.5298 | -5.7456 | 0.5250 | -3.6300
12 | 17.72% | -5.8079 | 0.5485 | -2.4826 | -5.1231 | 0.5366 | -3.0849
16 | 16.66% | -5.7466 | 0.5460 | -2.4582 | -4.8532 | 0.5250 | -2.7622
32 | 15.35% | -5.6147 | 0.5401 | -2.4191 | -3.9763 | 0.5264 | -1.9117
48 | 15.02% | -5.5350 | 0.5378 | -2.3999 | -3.3577 | 0.5440 | -1.3557

Figure 9 shows the relationship between the means and $\Delta t_{max}$ for
these two cases, fitted with a logarithmic function for each. The same is that
a small $\Delta t_{max}$ corresponds to a small imaging error, while case 2
has a larger variation range of the means than case 1. This difference is
probably due to the different characteristics of the amplitude intensity
distribution of visibilities on the $uv$ plane in these two cases. Case 1 is
relatively uniform, while case 2 is more concentrated in the low-frequency
part.

Considering both the compression ratio and dirty image quality, a small
$\Delta t_{max}$ (e.g. $\Delta t_{max}$ = 12.) could meet the requirements of
common use. We also found that further compression over short baselines is the
cause of the relative error in the dirty images. A large $\Delta t_{max}$
implies large deviations in the recovered visibilities on short baselines. In
practice, it is difficult to invert a suitable $\Delta t_{max}$ from the
imaging results, while choosing a small one is feasible and safe.

Figure 8: The distribution of the absolute value of the residual with
brightness for the pure noise model results, and $\Delta t_{max}$ = 48.

Figure 9: The trend of the different mean values with $\Delta t_{max}$ in two
groups.

###  3.4 The Processing Performance

The processing performance of BDA is a fundamental metric. A series of tests
were performed on a Centos 7 server equipped with 32 processors (Intel Xeon
Gold 6226R), 2.9-GHz core frequency, and 1024 GB of RAM. The version of RASCIL
used to obtain the simulation data in the tests was v.0.1.11. Using a
12-minute single-channel simulation data of SKA1-LOW with a data volume of
12.5 GB, we profile the BDA module optimized by the Numba, pure Python code,
and Pandas, respectively. The performance results are presented in Table 2.
The BDA module implemented using the Numba has the best performance.

Table 2: The performance of the BDA implementation at different data volumes. Data volume | Numba | Pandas | Pure Python
---|---|---|---
(GB) | (seconds) | (seconds) | (seconds)
1.56 | 10.75 | 13.69 | 400.01
3.13 | 17.48 | 23.78 | 697.62
6.25 | 29.02 | 44.21 | 1413.21
12.50 | 58.79 | 86.18 | 2941.89
25.00 | 105.59 | 171.38 | 6523.35
50.00 | 243.12 | 367.48 | 16178.32

Therefore, we further tested the processing performance of the Numba optimized
code with a series of simulated data of different observation times and four
channels. The maximum observation time was 48 minutes, and the data volume
obtained was 143.75 GB. During the process of BDA, $\Delta t_{max}$ was set to
6, 12, and 24 respectively, and $\Delta f_{max}$ was 1. The performance result
is shown in Figure 10, fitted with the exponential function respectively.

Figure 10: The performance of the BDA implementation optimized with Numba at
different data volumes.

As the results are shown in Figure 10, the time consumption of BDA processing
is essentially linear with the amount of data to be processed. The processing
speed of BDA in the case of a single process is about 13 GB per minute.
Further improvements are expected under parallel computing conditions. This
speed is acceptable for SRCs data pre-processing. In addition, the processing
speed has little to do with the amount of $\Delta t_{max}$.

##  4 Discussions

Experimental results showed a significant decrease in the capacity of the
visibility data with different $\Delta t_{max}$, which indicates that BDA is
very valuable for SKA data processing, especially for the construction of
subsequent SRCs. The annual storage capacity of SRCs will be at least 5
petabytes (PB) per year in the beginning, which will increase to at least 1.7
Exabytes by 2028, with at least 700 PB online (Bolton 2019). BDA can compress
at least 50% of the short baseline visibility data, which is valuable for
reducing the cost of SRCs construction.

This study also gave the variation of the imaging quality at different $\Delta
t_{max}$. The experimental results provided an essential reference for
SKA1-LOW to carry out EoR and cosmic dawn. The MWA data were averaged
(Mitchell et al. 2008), but no specific details were given.

##  5 Conclusions

We implemented a new BDA module based on RASCIL, and this implementation is
created by designing functions using Numba. It has not introduced excessive
memory usage in the tests and completes the computation tasks faster than the
other modules. The speed of a single process is around 13 GB per minute. It
also performed well during the BDA processing of the simulation data.

Through the simulation of observing the GLEAM and M31 model images with
SKA1-LOW, we evaluate the performance of the BDA. According to the subsequent
analysis, the error due to BDA increases with the maximum upper limit of the
averaging interval on the short baselines. In contrast, the compression ratio
does not improve all the time, and the reduction in data volume remains at a
maximum ratio of approximately 85%. The smaller upper limit is sufficient for
the compression ratio, and the imaging error is reasonable. Overall, the BDA
technology will have applications in the face of massive SKA observation data
processing. The BDA can effectively reduce the storage space of visibility
data, as it is also valuable for the future construction and application of
SRCs.

###### Acknowledgements.

This work is supported by the National SKA Program of China (2020SKA0110300),
the Joint Research Fund in Astronomy (U1831204, U1931141) under cooperative
agreement between the National Natural Science Foundation of China (NSFC) and
the Chinese Academy of Sciences (CAS), the Funds for International Cooperation
and Exchange of the National Natural Science Foundation of China
(11961141001), the National Natural Science Foundation of China (No.11903009).
This work is also supported by Astronomical Big Data Joint Research Center,
co-founded by National Astronomical Observatories, Chinese Academy of Sciences
and Alibaba Cloud. We do appreciate the anonymous referee for valuable and
helpful comments and suggestions.

## References

  * Atemkeng et al. (2018) Atemkeng, M., Smirnov, O., Tasse, C., et al. 2018, MNRAS, 477, 4511
  * Bolton (2019) Bolton, R. C. 2019, SKA Regional Centre Requirements, Tech. rep., SKA-TEL-SKO-0000735, internal SKA document
  * Braun (1996) Braun, R. 1996, in Astrophysics and Space Science Library, Vol. 208, The Westerbork Observatory, Continuing Adventure in, ed. E. Raimond & R. Genee, 167
  * Bridle & Schwab (1999) Bridle, A. H., & Schwab, F. R. 1999, in Astronomical Society of the Pacific Conference Series, Vol. 180, Synthesis Imaging in Radio Astronomy II, ed. G. B. Taylor, C. L. Carilli, & R. A. Perley, 371
  * Cotton (1986) Cotton, W. D. 1986, in Synthesis Imaging, ed. R. A. Perley, F. R. Schwab, & A. H. Bridle, 123
  * Cotton (1999) Cotton, W. D. 1999, in Astronomical Society of the Pacific Conference Series, Vol. 180, Synthesis Imaging in Radio Astronomy II, ed. G. B. Taylor, C. L. Carilli, & R. A. Perley, 357
  * Hurley-Walker et al. (2017) Hurley-Walker, N., Callingham, J. R., Hancock, P. J., et al. 2017, MNRAS, 464, 1146
  * Lam et al. (2015) Lam, S. K., Pitrou, A., & Seibert, S. 2015, in Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC, LLVM ’15 (New York, NY, USA: Association for Computing Machinery)
  * Mitchell et al. (2008) Mitchell, D. A., Greenhill, L. J., Wayth, R. B., et al. 2008, IEEE Journal of Selected Topics in Signal Processing, 2, 707
  * Salvini & Wijnholds (2017) Salvini, S., & Wijnholds, S. J. 2017, in 2017 XXXIInd General Assembly and Scientific Symposium of the International Union of Radio Science (URSI GASS), 1
  * van der Walt et al. (2014) van der Walt, S., Schönberger, J. L., Nunez-Iglesias, J., et al. 2014, arXiv e-prints, arXiv:1407.6245
  * Wei et al. (2021) Wei, S.-L., Luo, K.-D., Wang, F., Deng, H., & Mei, Y. 2021, Research in Astronomy and Astrophysics, 21, 080
  * Wijnholds et al. (2018) Wijnholds, S. J., Willis, A. G., & Salvini, S. 2018, MNRAS, 476, 2029


#  Time-Reversible Thermodynamic Irreversibility :
One-Dimensional Heat-Conducting Oscillators
and Two-Dimensional Newtonian Shockwaves

William Graham Hoover and Carol Griswold Hoover
Ruby Valley Research Institute
Highway Contract 60, Box 601
Ruby Valley, Nevada 89833 ;

###### Abstract

We analyze the time-reversible mechanics of two irreversible simulation types.
The first is a dissipative one-dimensional heat-conducting oscillator exposed
to a temperature gradient in a three-dimensional phase space with coordinate
$q$, momentum $p$, and thermostat control variable $\zeta$. The second type
simulates a conservative two-dimensional $N$-body fluid with $4N$ phase
variables $\\{q,p\\}$ undergoing shock compression. Despite the time-
reversibility of each of the three oscillator equations and all of the $4N$
manybody motion equations both types of simulation are irreversible, obeying
the Second Law of Thermodynamics. But for different reasons. The irreversible
oscillator seeks out an attractive dissipative limit cycle. The likewise
irreversible, but thoroughly conservative, Newtonian shockwave eventually
generates a reversible near-equilibrium pair of rarefaction fans. Both problem
types illustrate interesting features of Lyapunov instability. This
instability results in the exponential growth of small perturbations, $\propto
e^{\lambda t}$ where $\lambda$ is a “Lyapunov exponent”.

Time Reversibility, Irreversibility, Heat Conduction, Shockwaves,
Thermodynamics’ Second Law

##  I Reversibility of Dissipative and Conservative Mechanics

Classical mechanics is time-reversible in the sense that a movie of the
motion, run backwards, obeys exactly the same motion equations as does the
original forward version. Classical mechanics is an excellent model for
conservative systems free of the real-life dissipative effects of friction,
viscosity, and heat conduction. In order to model dissipative phenomena on an
atomistic scale nonequilibrium molecular dynamics includes control variables
in the equations of motion. These variables use feedback to impose local
values of the temperature and pressure which drive nonequilibrium flows.
Thousands of implementations of this approach have been stimulated by Shuichi
Nosé’s pioneering 1984 workb1 ; b2 . We will explore the time reversibility of
an application of his work here.

Nosé’s 1984 papers generalize Hamiltonian mechanics with a frictional variable
$\zeta$ controlling the kinetic temperature $T$ of one or more particular
degrees of freedom :

$\zeta\longrightarrow kT=\langle\ p^{2}/m=mv^{2}\ \rangle\ .$

Here $k$ is Boltzmann’s constant, $m$ the mass of a particle, and $mv=p$ the
momentum of a controlled degree of freedom. For simplicity in what follows we
set both $k$ and $m$ equal to unity. Hoover applied Nosé’s idea to the
simplest special case, a one-dimensional harmonic oscillator, in 1985b3 ,
later extending that work in 1986 with Posch and Veselyb4 . Numerical
solutions of the thermostatted oscillator’s equations of motion (with $k$ and
$m$ unity now) ,

$\\{\ \dot{q}=p\ ;\ \dot{p}=-q-\zeta p\ ;\ \dot{\zeta}=p^{2}-T\ \\},$

can be reversed in either of two ways, by [1] changing the signs of the time
and timestep, $t$ and $dt$, or [2] by changing the signs of $p$ and the
control variable $\zeta$. Both ways simply reverse the time dependence of the
coordinate: $q(+t)\rightarrow q(-t)$.

##  II One-Dimensional Heat-Conducting Oscillator

Figure 1:  Two conservative tori are interlinked with a dissipative limit
cycleb6 . All three stationary solutions of the oscillator equations for
$T(q)=1+0.42\tanh(q)$ can be generated with initial values $(q,p,\zeta)$ =
(-2,-2,0), (-2.3,0,0) and $(+3.5,0,0)$. If instead $T(q)=1+\tanh(q)$ all
initial conditions lead to the limit cycle illustrated in Figures 2-4.

In 1997 Posch and Hoover generalized the oscillator problem, specifying a
coordinate-dependent temperature $T(q)=1+\epsilon\tanh(q)$. This temperature
profile has a maximum temperature gradient, $(dT/dq)=\epsilon$ at $q=0$.
Particular choices of $\epsilon$ generated a variety of $(q,p,\zeta)$ “strange
attractors” [fractal distributions in $(q,p,\zeta)$ space]. Figure 5, a cross-
section through a fractal attractor, gives an impression of the complicated
structures resulting from relatively simple ordinary differential equations.
The fractional dimensionalities of these attractor distributions were all
between 2 and 3. Other initial conditions or choices of $T(q)$ resulted in
one-dimensional limit cycles rather than fractalsb5 . More recently, in 2014,
Sprott and the Hoovers found initial conditions, $(q,0,0)$, with
$T(q)=1+0.42\tanh(q)$, which generate two distinct families of conservative
tori (with the initial values $q=-2.3$ and 3.5)b6 , The tori are interlinked
and coexisting stably with a one-dimensional dissipative limit cycle. The
cycle can be generated easily with initial values $(q,p,\zeta)=(-2,-2,0)$ and
a fourth-order Runge-Kutta timestep $dt=0.01$. Figure 2 in Reference 6 shows
the three interlinked phase-space structures. We reproduce it here as Figure
1.

Figure 2:  Stable and unstable limit cycles for a heat-conducting one-
dimensional harmonic oscillator. At the left the time dependence of coordinate
$q$, momentum $p$, and friction coefficient $\zeta$, purple, green, and blue
respectively, for the conducting oscillator with period $\tau=13.7494$. At the
right are three-dimensional plots of the attractive orbit (green), with
initial condition $\\{q,p,\zeta\\}=\\{0,0.1050726,0.1455481\\}$ and the
corresponding unstable repulsive orbit (red), with $p$ and $\zeta$ changed in
sign. The mean value of the (hot to cold) heat current is $\langle\ (p^{3}/2)\
\rangle=-0.74383$. The red and green circles indicate the four crossings of
the orbits with the $\zeta=0$ grey plane.

In this rich collection of one-dimensional limit cycles, two-dimensional tori,
and fractional-dimensional strange attractors the simplest special case is
arguably $T(q)=1+\tanh(q)$. The coordinate-dependent temperature varies from 0
to 2 as $q$ varies from $-\infty$ to $+\infty$. We believe that the basin of
attraction for this case is the entire three-dimensional phase space. To
support this idea we chose a square $200\times 200$ grid of $(q,p,\zeta)$
points in the $\zeta=0$ plane with $q$ and $p$ ranging from $-4.975$ to
$+4.975$ in steps of 0.05. For each of these 40,000 initial conditions we
generated an orbit of length 100$\tau$, and plotted the $(q,p)$ values
whenever $\zeta$ changed sign. Every one of these long orbits ended up
crossing at the two penetration points plotted in green in Figure 2. For each
orbit we used a timestep of $\tau/10000$ where the cycle period is
$\tau=13.7494$. Let us explore that solution in more detail, based on fourth-
order Runge-Kutta numerical simulations.

Figure 3:  Time dependence of the local Lyapunov exponent for the oscillator
during its 99th and 100th period of oscillation. The local value is purple.
The green time average of $\lambda_{1}(t)$ vanishes, corresponding to the lack
of averaged strain along the trajectory.

From the phase-space analog of Liouville’s continuity equation the mean value
of the friction coefficient $\zeta$ necessarily corresponds to the time-
averaged loss rate of phase volume, $\otimes=dqdpd\zeta$:

$\langle\ (\dot{\otimes}/\otimes)\ \rangle=\langle\ (\partial\dot{q}/\partial
q)+(\partial\dot{p}/\partial p)+(\partial\dot{\zeta}/\partial\zeta)\
\rangle=\langle\ 0-\zeta+0\ \rangle=-1.325\ .$

In a single period the comoving phase volume decreases by a factor of
$e^{\langle\ \zeta\ \rangle\tau}=e^{1.325\times 13.7494}\simeq 10^{8}$. The
maximum temperature gradient, $(dT/dq)$, is unity, at $q=0$. The mean heat
current, averaged over time, $\langle(p^{3}/2)\rangle$, is $-0.74383$, and the
net transport of kinetic energy $(p^{2}/2)$ is from right to left, consistent
with thermodynamics’ Second Law.

The Lyapunov exponents, three of them in a three-dimensional phase space,
measure the comoving expansion rates of the phase volume $\otimes$ :

$(\dot{\otimes}/\otimes)=\lambda_{1}(t)+\lambda_{2}(t)+\lambda_{3}(t)\ .$

The one-dimensional limit cycle’s largest Lyapunov exponent has an average
value of zero, as shown in Figure 3. The vanishing mean value of $\lambda_{1}$
corresponds to the averaged lack of relative motion of two adjacent trajectory
points along the attractive one-dimensional trajectory. $\lambda_{2}$ and
$\lambda_{3}$ have negative averages, $-0.110$ and $-1.215$, describing the
net rates of convergence of nearby trajectories in the two directions
perpendicular to the limit cycle.

For the special case $T(q)=1+\tanh(q)$ the longtime solution of the motion
equations forward in time is the unique attractive periodic orbit shown
(green) in the right panel of Figure 2. The period $\tau$ is 13.7494 with
$-2.616<q<+0.198$. The reversed orbit, with the same range of $q$ visited in
the opposite time direction, is half the attractor/repellor pair. The repellor
is only observable briefly due to its inherent Lyapunov instability,
proportional to $e^{\lambda_{1}(t)\times t}\simeq e^{+1.215t}$. We expect to
see the exponential growth of an original one-step roundoff error grow to
observability in just a few oscillator periods. We examine that next.

Figure 4:  Evolution of ten $(q,\zeta)$ orbits beginning at the initial upper
black circle with the conditions of Figure 2. At time $t=137.494$ the sign of
the timestep is reversed. The reversed trajectory (black) traces back
accurately for nearly three reversed orbits before making a rapid switch to
the stabilized repellor, following it to the lower black circle at $t=0$.

##  III Stability and Instability of Periodic Oscillator Orbits

Figure 4 illustrates the stability of the attractive limit-cycle orbit,
$\\{\lambda_{i}\\}=\\{0,-0.110,-1.215\\}$ and the instability of its time
reversal, $\\{\lambda_{i}\\}=\\{+1.215,+0.110,0\\}$ with a two-stage
simulation. First we follow ten circuits of the attractor in green, using a
million quadruple-precision fourth-order Runge-Kutta timesteps of
$dt=10^{-5}\tau=0.000137494$, implying a local integration error of order
$dt^{5}/120\simeq e^{-49}$ at each step. We then reverse time,
$+dt\rightarrow-dt$, converting the stable attractor to the unstable reversed
repellor. Unlike the stable attractor the repellor is unstable, with a
positive Lyapunov exponent $\lambda_{1}=1.215$. This suggests that the
reversal should become visibly inaccurate at a time of order
$49.3/1.215\simeq$ three or four periods. This agrees well with the reversed
black trajectory of Figure 4, which follows the reversed $(q,p,\zeta)$ orbit
(stabilized by the negative $dt$) for between six and seven periods, ending up
at the lower black circle :

$\\{\ q=0,\ p=-0.1050726,\ \zeta=-0.1455481\ \\}\ .$

The largest Lyapunov exponent is relatively easy to measure. Follow two
neighboring trajectories, the “reference” and the “satellite”, adjusting the
satellite location by rescaling its separation from the reference,
$\delta_{t+dt}\rightarrow\delta_{\rm o}$, at the conclusion of each time step.
This rescaling precisely counters the exponential growth which would occur in
the absence of rescaling. The local Lyapunov exponent follows from the
rescaling operation :

$\lambda_{1}dt=-\ln(\delta_{\rm after}/\delta_{\rm
o})\rightarrow\lambda_{1}(t)=(-1/dt)\ln(\delta_{t+dt}/\delta_{\rm o}).$

A convenient choice for $\delta_{\rm o}$ is $0.00001$. Figure 3 shows the time
variation of the largest Lyapunov exponent $\lambda_{1}$, which lies in the
range $-3.00<\lambda_{1}(t)<+3.48$ with a mean value of zero, corresponding to
the (lack of) growth rate of perturbations parallel to the trajectory.

##  IV Time Reversibility and Loschmidt’s Paradox

A classic physics puzzle addresses the surprising coexistence of macroscopic
irreversibility with microscopic time reversibility. In 1876 Loschmidt pointed
out that any solution of the equations of motion which is time reversible and
demonstrates the production of entropy can be made to violate the Second Law
of Thermodynamics by analyzing the reversed motion. Simply stated, time-
reversible mechanics necessarily violates the Second Law in one of the two
time directions. We have already seen that time-reversible Nosé-Hoover
mechanics, with control of the kinetic temperature, obeys the Second Law. This
is not only possible, but inevitable, for computational models of the heat-
conducting oscillator. With $dt$ positive, attractive distributions such as
the oscillator limit cycle are inevitably observed. Repellors are not, due to
their vanishing probabilities.

The exploration of a simple one-body time-reversible modelb7 , and its
relation to entropy production in a many-body system with heat flowb8 ,
clarified the reversibility paradox in 1987. The “Galton Staircase” pictures a
reversibly thermostatted particle in a downhill steady state, driven by a
periodic sinusoidal potential superimposed on a constant downhill field. As
the particle falls more than it climbs, the model generates a fractal
(fractional dimensional) phase-space distribution in its three-dimensional
phase space $(q,p,\zeta)$. The resulting zero-volume attractive fractal, when
reversed, corresponds to the extreme rarity of states violating the Second Law
of Thermodynamics by converting heat to work. The mirror-image repulsive
fractal, corresponding to an upward moving particle violating the Law has,
like the attractor, zero volume, but is repulsive and of zero probability. The
“attractor” is attractive, with probability one. This difference in behavior
occurs because Lyapunov instability is not time-reversible.

In the Galton staircase a particle travelling uphill, as described by the
repellor states, violates the Law by converting kinetic energy to potential.
Repulsion, coupled with zero volume, makes these fractal repellor states
unobservable. The conducting oscillator of Figures 2-4 offers a simple analog
for a particle transporting energy from hot to cold rather than transporting
mass through motion driven by a gravitational field. Both systems resolve
Loschmidt’s Paradox by introducing time-reversible variables controlling
temperature. It is the extreme unobservable rarity of repellor states, the
fractal set for uphill motion in the Galton Staircase, and here the repulsive
one-dimensional repellor limit cycle in three-dimensional space, that forces
motions to obey the Second Law. Staircase simulations reveal the exponential
growth of the separation from the fractal repellor and an irresistible
attraction to the repellor’s mirror-image attractor. Likewise the conducting
oscillator with $T=1+\tanh(q)$ follows the attractive limit cycle of the
Figures, rather than the cycle’s mirror image repulsive twin, which repels
rather than attracts.

Holian, Hoover, and Poschb8  stressed that similar irreversible behavior
occurs in the reversible simulations of thermostatted nonequilibrium manybody
systems. They described a heat-conducting system in contact with two
reservoirs, one hot and the other cold. Such a system loses phase volume when
it satisfies Fourier’s Law, transmitting heat from the hot reservoir (with an
entropy production $Q/T_{H})$ to the cold (with an entropy loss $Q/T_{C})$
which necessarily exceeds the gain. The result is a phase-volume loss
exponential in the time. Just as in the Galton Staircase mass-flow problem
heat flow from hot to cold results in fractal phase-space structures. Both
fractal types have zero volume, with zero probability of observing the
repellor and with inevitable longtime probability one for the attractor.

##  V An Illustrative Fractal for the Conducting Oscillator

To help visualize the attractors and repellors that characterize
nonequilibrium systems we consider here a fractal resulting from dissipation
controlled by a pair of control variables. $\xi$ controls the fourth moment of
the velocity distribution while $\zeta$ controls the second. Because fractal
distributions are difficult to visualize in their entirety they are typically
described by projections or cross sections. The four-dimensional conducting
oscillator introduced by Posch and Hooverb5  provides a variety of thought-
provoking fractal structures. As an example, for the same temperature profile
considered here, $T=1+\tanh(q)$, see Figure 5. The additional phase-space
dimension results from using two control variables rather than one. The
doubly-thermostatted oscillator requires the solution of four ordinary
differential equations:

$\\{\ \dot{q}=p\ ;\ \dot{p}=-q-\xi p^{3}-\zeta p\ ;\ \dot{\xi}=p^{4}-3p^{2}T\
;\ \dot{\zeta}=p^{2}-T\ \\}\ .$

These motion equations are fully ergodic for the special case in which $T=1$.
That is, all possible values of the four variables occur and with the known
distribution:

$4\pi^{2}{\rm prob}(q,p,\xi,\zeta)=e^{-(q^{2}+p^{2}+\xi^{2}+\zeta^{2})/2}\ .$

Figure 5:  This double cross section is generated by plotting $(q,p)$
trajectory points whenever the two thermostat values are near zero. The Figure
includes ten million points satisfying the condition
$\xi^{2}+\zeta^{2}<0.0001$. Fourth-order Runge-Kutta integration with
$dt=0.001$ was used.  Figure 6:  At the left we see five particle plots, at
times 0, 10, 20, 30, and 40 from bottom to top, during the twofold shock
compression $(0<t<20)$ and subsequent generation of rarefaction waves $(t>20)$
simulated by the headon collision of $2\times 840=1680$ soft disks. The
initial square-lattice coordinates (at bottom left) and horizontal velocities
$u_{p}=u_{s}/2=\pm 0.875$ (at bottom right) are selected so as to generate
twofold shock compression, doubling the density to generate a hot fluid,
reaching it with the coordinates and velocities in the middle view of Figure 6
and the bottom view of Figure 7. Initially the lattice symmetry is broken with
random displacements of $x\ $ and $y$ in the range from $-0.05$ to $+0.05$.
The top two plots show the growing twin rarefaction fans launched from the
edges of the hot dense fluid. At the right the evolving horizontal velocity
components are shown at the same five equally-spaced times, from 0 to 40. The
rarefaction fans at times 30 and 40 are expanding at about the speed of sound.
The smooth short-ranged repulsive pair potential governing the Newtonian
dynamics is $\phi(r<1)=(10/\pi)(1-r)^{3}$. The fourth-order Runge-Kutta
integrator with $dt=0.01$ was used for Figures 6 and 7.

##  VI Two-Dimensional Strong Shockwaves Do Not Reverse

We have seen that Nosé-Hoover temperature control provides a probabilistic
mechanism for irreversibility, the formation of attractor-repellor pairs in
phase space which stabilize the attractor and destabilize the repellor, both
through Lyapunov instability. We have recently detected a related mechanism,
but free of control variables, demonstrating the irreversibility of purely
Newtonian mechanical systems, illustrated here with an atomistic model
generating strong shockwaves. Shockwaves are localized regions, usually just a
few atomic diameters in width, within which density, pressure, energy, and
temperature all undergo substantial increases. We consider a model here where
the temperature changes by a factor of 100 and the density doubles. Shockwaves
are relatively easily treated computationally because they are bounded by
equilibrium thermodynamic states. To ease the computational burden we consider
the shockwave compression of two-dimensional soft-disk particles in two space
dimensions. The sudden compression occurs in the $x$ direction. The purely-
repulsive pair potential is chosen for simplicity,
$\phi(r<1)=(10/\pi)(1-r)^{3}$. The initial near-zero-pressure state is a
nearly perfect square lattice, with lattice spacing of unity. The small
initial displacements in the range $\pm 0.05$ correspond to a temperature of
order 0.001.

Consider the head-on collision of two $N$-body mirror-image zero-pressure
zero-energy blocks of material with opposite velocitiesb9 ; b10 ; b11 , shown
in Figure 6. Here $N=35\times 24=840$. With periodic boundaries in $y$ the two
colliding $N$-body blocks steadily convert their kinetic energy to heat. At
any stage in the simulation a reversed solution will show, briefly,
antithermodynamic behavior, converting some of the internal energy of the
stagnating blocks back into the original directed kinetic energy,
$(p_{x}^{2}/2)=(0.875^{2}/2)$ per particle.

Berni Alder and Marvin Ross emphasized the highly irreversible nature of
shockwavesb12  as follows : “the most irreversible way to go from one
thermodynamic state to another”. A simple example of this transformation is
illustrated in Figure 6, where $70\times 24=1680$ particles undergo twofold
compression and then expand to form a symmetric pair of rarefaction fans. The
initial condition for this example is a neighboring mirror-image pair of
colliding square-lattice blocks, both at the stress-free density of unity. The
left half travels rightward and the right half leftward. Periodic boundaries
are imposed in the $y$ direction, at the top and bottom of the two colliding
840-body blocks. The dynamics is purely Newtonian. The difference between the
steady shockwaves forward in time and the unsteady rarefaction (rather than
shocks) waves in the reversed time direction of Figure 7 shows that the
shockwaves are irreversible. In fact the irreversible Navier-Stokes equations
of motion predict that a reversed shockwave will immediately widen and slow,
transforming into an unsteady rarefaction fanb9 ; b10 .

Figure 7 was constructed by reversing the velocities of all particles in
Figure 6 at the time 20, the time of maximum twofold compression. Notice that
the snapshot second from the bottom of Figure 7, where the flow has been
reversed so that the configuration is only halfway to maximum compression,
resembles closely that second from the bottom in Figure 6, where the flow is
forward, and halfway to the time of maximum compression. This apparent
reversibility suggests that the initial single-step Runge-Kutta integration
error,

$dt=0.01\rightarrow 0.01^{5}/5!=e^{-27.8}\ {\rm or}\ dt=0.005\rightarrow
200^{-5}/5!=e^{-31.3}\ ,$

expands exponentially to become of order unity at $t=27.8/\lambda$ or
$31.3/\lambda$ in these two typical cases.

Attempting to confirm and elaborate this estimate we constructed reference and
satellite trajectories, for a range of timesteps from $n=1$ through $n=7$ :

$dt=0.1/2^{n}\ ;\ n=1\rightarrow dt=0.05\ ...\ n=7\rightarrow dt=0.00078125\
,$

rescaling their reference-to-satellite separation to 0.00001 at the end of
each timestep. Typical resulting values of the local-in-time Lyapunov exponent
$\lambda_{1}(t)\simeq 1.8$ are shown in Figures 8 and 9. The complete set of
6720 two-block motion equations is included in the Lyapunov calculation. The
signs of all the velocity components are changed at time 20 in Figure 8 and at
time 10 in Figure 9, with the particle coordinates pictured as the latter
figure’s inset.

Fluctuations of the local exponent $\lambda_{1}(t)$ can be substantially
reduced by smoothing, averaging the nearest 100 local values. In Figure 8 this
average is plotted in red. With $dt$ positive the exponent $\lambda_{1}(t)$ is
uniformly positive, changing sign with the velocity reversal. The reversed
reference and satellite trajectories attract for a while, as expected, but
only for a while, for a time of order 3 with $dt=0.01$ and 5 with
$dt=0.00078125$. Because integration errors are magnified exponentially, as
described by the Lyapunov exponent, $\lambda_{1}\simeq 1.8$, there is a
systematic timestep dependence of the time “$t_{rev}(dt)$”. The exponent is
reversed for a time proportional to the logarithm of the number of timesteps,
$\propto\ln(1/dt)$. Roughly speaking, $t_{rev}(dt)$ increases by about 0.4 for
each halving of the timestep. See Figure 9. A set of seven computations with
$0.05\leq dt\leq 0.0078125$ suggests the phenomenological relation

$e^{1.8t_{rev}(dt)}\simeq dt^{-?}\rightarrow 1.8t_{rev}\simeq-?\ln(dt)\ .$

In Figure 9 $t_{\rm rev}$ varies from about 3 to about 5 as the timestep
decreases by a factor of 64.

The logarithmic relationship between timestep and reversal time is consistent
with the exponential amplification of few-step integration errors. A similar
rough relationship holds for the limit-cycle problem of Figures 2-4. Capturing
a reversed trajectory with visual accuracy up to a time $t_{rev}(dt)$ requires
a simulation effort varying as $(t_{\rm rev}/dt)$.

We are pleased to offer a reward for further investigation of these problem
types: a special $1000 “Snook Prize”, in honor of our late colleague Ian
Snook, who died in 2013. Application for this Prize requires the submission of
an appropriate acceptable electronic manuscript addressing shockwave
reversibility to CMST.eu prior to year’s end of 2023.

Figure 7:  Velocity reversal following the twofold shock compression of 1680
soft disks. The initial velocities $\pm 0.875$ at the left and right, resulted
in a headon collision. The initial particle coordinates at the base of this
Figure correspond precisely to the maximum-compression middle view of Figure
6. The initial velocities were changed in sign so that the next-to-bottom
snapshot here corresponds closely to a reversal of the next-to-bottom view of
Figure 6. Soon after, Lyapunov instability, with $\lambda_{1}\simeq 2$,
prevents additional reversed configurations at time 20, 30, and 40. Just as in
Figure 6 particle positions appear at the left with corresponding horizontal
velocity components at the right.  Figure 8:  The local Lyapunov exponent for
the time-reversal of Figure 6. The reversal, at time 20, is indicated by a
filled red circle. The red curve corresponds to a smoothed exponent averaged
over unit time, 100 time steps with the step $dt=0.01$. Notice that the
Lyapunov exponent is of order $\pm 2$ in the shocked dense fluid or in the
cold initial solid, suggesting an observable lack of reversibility when the
amplified roundoff error reaches the amplitude of Lyapunov instability, a time
of order 14 for $dt=0.01$ : $0.01^{5}e^{\lambda_{1}(t)\times t}/120\simeq
1\rightarrow t\simeq 14\ .$

Figure 9:  The largest Lyapunov exponent, $\lambda_{1}\simeq 1.8$ for seven
simulations with $dt=0.1/2^{n}$ reflects the velocity sign change at time 10.
The accuracy of the reversal persists for a time $t_{\rm rev}$ varying as
$\ln(1/dt)$, consistent with the Lyapunov amplification of integration errors
to visible size. The inset snapshot corresponds to the $t=10$ configuration of
1680 soft-disk particles undergoing twofold shock compression up until the
reversal time indicated by the red dot.  Figure 10:  Smoothed-particle density
profiles, Newtonian at the left and reversed at the right. These profiles
correspond to the topmost configurations of Figures 6 and 7. Similar profiles
of $P_{xx}$ and $P_{yy}$ show that the pressure is nearly isotropic. The thick
black and red lines correspond to smoothed-particle ranges of 2 and 3
respectively.

##  VII Summary and Prognosis

Our first problem, the thermostatted generation of a heat current, stabilized
by Nosé-Hoover mechanics, demonstrates that such dissipative examples, despite
their time reversibility, stabilize mass, momentum, and energy flows
consistent with the Second Law of Thermodynamics. The mechanism for this
dissipative irreversibility is the formation of fractal phase-space
distributions with mirror-image fractal pairs of distributions. With these
pairs the attractor has probability one while the repellor has probability
zero, even for a small system, the heat-conducting oscillator. This same
explanation of thermostatted irreversibility holds for manybody systems, as
was suggested in References 7 and 8.

The purely-Newtonian shockwave problem illustrates a different, but related,
mechanism for irreversibility. Though “mechanism” is conceptual, its
realization is necessarily computational. Lyapunov instability destroys a
time-reversed system’s memory by magnifying computational noise. Soon the
amplified roundoff error becomes of order unity. Our shockwaves simulations
show that this time is of order $t\simeq$ 3 to 5, just a few collision times.
At that time the reversed flows are destabilized by the preponderance of
entropy-producing flows over entropy-reducing flows, giving a purely Newtonian
resolution of Loschmidt’s paradox. No thermostatting control variable is
required. But still it is likely that the probability of going forward with an
entropy-producing shock exceeds that of its reversal in a singular way. The
preponderance of states generating entropy, in the shock, over those which
would reduce it in the reversed shock, requires no modification of the
Newtonian motion equations.

In making the connection between computational simulations and continuum
mechanics smooth-particle averaging is an indispensable tool. To illustrate
this idea let us calculate density profiles $\rho(x)$ for the topmost
snapshots in Figures 6 and 7. Imagine the density of every particle to be
spread out in $x$ over a range $h$ according to a normalized weight function
$w(r)$ with finite range and two continuous derivatives:

$r<h\longrightarrow w(r,h)=[5/4hL_{y}][1-(r/h)]^{3}[1+3(r/h)]\ .$ $[\ {\rm
Lucy^{\prime}s\ weight\ function}\ ]$

Summing the weights contributed by each of the particles to every grid point
in the set $\\{x_{g}\\}$ produces the smoothed density profiles shown in
Figure 10 with $h=2$ and $h=3$. The grid spacing is unity so that every
particle contributes to four or 6 nearby grid points for these two values of
$h$. We see that neither profile matches the uniform density at the base of
Figure 6. The regular square-lattice structure near both ends of the hot
shocked fluid (near $x=\pm 20$, at the base of Figure 7 and the middle of
Figure 6) is less susceptible to the smoothing loss of memory due to Lyapunov
instability.

Smooth-particle averaging can be applied to any of the atomistic functions of
coordinates and momentum. Plots of the pressure tensor indicate isotropy with
$P_{xx}\simeq P_{yy}$. Research into the details of the atomistic distribution
functions could elucidate further the mechanism responsible for the
exponentially greater density of phase flows obeying the Second Law to those
flouting it. There remains much to do in understanding the failure of
Loschmidt’s cogent idea—questioning the ability of mechanics to provide an
understanding of thermodynamics.

##  VIII Acknowledgments

We thank Professor Krzysztof Witold Wojciechowski (Poznañ) for his interest
and support of this work. He spotted typos in the first version of this
manuscript, for which we are grateful. We also thank Professor Julien Clinton
Sprott (Wisconsin) for stimulation during our preparation of Reference 11.

## References

  * (1) S. Nosé, “A Unified Formulation of the Constant Temperature Molecular Dynamics Method”, The Journal of Chemical Physics 81, 511-519 (1984).
  * (2) S. Nosé, “A Molecular Dynamics Method for Simulations in the Canonical Ensemble”, Molecular Physics 52, 255-268 (1984).
  * (3) Wm. G. Hoover, “Canonical Dynamics. Equilibrium Phase-Space Distributions”, Physical Review A 31, 1695-1697 (1985).
  * (4) H. A. Posch, W. G. Hoover, and F. J. Vesely, “Canonical Dynamics of the Nosé Oscillator: Stability, Order, and Chaos”, Physical Review A 33, 4253-4265 (1986).
  * (5) H. A. Posch and W. G. Hoover, “Time-Reversible Dissipative Attractors in Three and Four Phase-Space Dimensions”, Physical Review E 55, 6803-6810 (1997).
  * (6) J. C. Sprott, W. G. Hoover, and C. G. Hoover, “Heat Conduction, and the Lack Thereof, in Time-Reversible Dynamical Systems: Generalized Nosé-Hoover Oscillators with a Temperature Gradient”, Physical Review E 89, 042914 (2014).
  * (7) W. G. Hoover, H. A. Posch, B. L. Holian, M. J. Gillan, M. Mareschal, and C. Massobrio, “Dissipative Irreversibility from Nosé’s Reversible Mechanics”, Molecular Simulation 1, 79-86 (1987).
  * (8) B. L. Holian, W. G. Hoover, and H. A. Posch, “Resolution of Loschmidt’s Paradox: The Origin of Irreversible Behavior in Reversible Atomistic Dynamics”, Physical Review Letters 59, 10-13 (1987).
  * (9) W. G. Hoover and C. G. Hoover, “Time-Symmetry Breaking in Hamiltonian Mechanics. Part II. A Memoir for Berni Julian Alder [1925-2020]”, Computational Methods in Science and Technology 26, 101-110 (2020).
  * (10) W. G. Hoover and C. G. Hoover, “Time-Symmetry Breaking in Hamiltonian Mechanics. Part III. A Memoir for Douglas James Henderson [1934-2020]”, Computational Methods in Science and Technology 26, 111-120 (2020).
  * (11) J. C. Sprott, W. G. Hoover, and C. G. Hoover, Elegant Simulations, (World Scientific, Singapore, 2023).
  * (12) M. Ross and B. J. Alder, “Shock Compression of Argon II. Nonadditive Repulsive Potential”, Journal of Chemical Physics 46, 4203-4210 (1967).