# Essentially non-hourglass and non-tensile-instability SPH elastic dynamics

Shuaihao Zhang<EMAIL_ADDRESS>Sérgio D.N. Lourenço<EMAIL_ADDRESS>Dong Wu<EMAIL_ADDRESS>Chi Zhang<EMAIL_ADDRESS>Xiangyu Hu
<EMAIL_ADDRESS>Department of Civil Engineering, The University of Hong
Kong, Pokfulam, Hong Kong SAR, China  School of Engineering and Design,
Technical University of Munich, 85748 Garching, Germany

###### Abstract

Since the tension instability was discovered in updated Lagrangian smoothed
particle hydrodynamics (ULSPH) at the end of the 20th century, researchers
have made considerable efforts to suppress its occurrence. However, up to the
present day, this problem has not been fundamentally resolved. In this paper,
the concept of hourglass modes is firstly introduced into ULSPH, and the
inherent causes of tension instability in elastic dynamics are clarified based
on this brand-new perspective. Specifically, we present an essentially non-
hourglass formulation by decomposing the shear acceleration with the Laplacian
operator, and a comprehensive set of challenging benchmark cases for elastic
dynamics is used to showcase that our method can completely eliminate tensile
instability by resolving hourglass modes. The present results reveal the true
origin of tension instability and challenge the traditional understanding of
its sources, i.e., hourglass modes are the real culprit behind inducing this
instability in tension zones rather that the tension itself. Furthermore, a
time integration scheme known as dual-criteria time stepping is adopted into
the simulation of solids for the first time, to significantly enhance
computational efficiency.

###### keywords:

Smoothed particle hydrodynamics; Hourglass modes; Tensile instability; updated
Lagrangian formulation; Elastic dynamics

††journal: Elsevier

##  1 Introduction

Smoothed particle hydrodynamics (SPH), original proposed by Lucy [1] and
Gingold and Monaghan [2] for simulating astrophysical problems, is a fully
Lagrangian particle-based method. In SPH, the physical quantities such
velocity, position, and stress are carried by each particle, and the motion of
particles is described in the Lagrangian framework, which is naturally well-
suited for simulating problems involving large deformations, especially
fracture and failure. Over the past 40 years, SPH has been extensively
developed and improved, and it has been successfully applied to simulate
various physical problems including fluid dynamics [3, 4], solid dynamics [5,
6], and fluid-solid interactions [7, 8].

Based on whether the particle configurations, which define the neighbors of
each particle, are updated during the simulation, the SPH methods for solid
dynamics can be classified into the total Lagrangian SPH (TLSPH) [9] and
updated Lagrangian SPH (ULSPH) [5, 10]. TLSPH is able to handle the elastic
and plastic dynamics efficiently as it saves the time required for updating
particle configurations. Compared with TLSPH, while ULSPH is able to cope with
material failure and fracture beyond elastic or plastic deformations more
naturally with updating particle configurations at each time step, it faces
two important drawbacks associated with elastic dynamics: one is the
persistent issue of tensile instability; the other is poorer efficiency due to
the computational effort and memory latency when the particle configurations
are updated frequently.

As shown in Fig. 1b, tensile instability in ULSPH is often associated with
particle clustering and numerical/artificial fractures, and was first studied
by Swegle et al. [11] in 1995. At that time, it was believed to be caused by
tension stress. Since then, different approaches have been proposed to address
this problem as reported in the literatures [12, 13, 14]. However, these
methods also have their own issues, such as failure to maintain conservation
properties, low computational efficiency, or limited applicability, being only
suitable for specific cases rather than universally applicable [15, 16, 17].
Later, in 2000, inspired by the repulsive interactions observed among closely
spaced atoms, Monaghan [10] introduced a small repulsive force (named the
artificial stress) between particles in SPH to prevent particle clustering and
then remove tensile instability. Based on the artificial stress and the signs
of principal stresses, Gray et al. [5] further developed this approach by
determining the parameters in artificial stress from the dispersion relation
for elastic waves.

Despite the broader recognition compared to many other approaches in
addressing tension instability, the artificial stress method still faces the
following limitations: (1) it has two parameters requiring case-dependent
tuning; (2) it may fail in scenarios where the deformation is significant or
when dealing with materials featuring high Poisson’s ratio [18, 19]; (3) its
extension for three-dimensional (3D) simulations is yet to be developed. The
third issue is due to the fact that the derivation of the artificial stress
term is based on a two-dimensional (2D) scenario [5]. To our best knowledge,
there are no documented instances utilizing the artificial stress in 3D
simulations.

Tensile instability is also found in SPH simulations of fluids, in which the
negative pressure, especially in vortical flows, leads to the generation of
artificial void regions. A commonly used strategy is applying a constant
background pressure [20] to keep positive pressure everywhere in the
simulations of incompressible flow with the weakly compressible SPH (WCSPH)
method. Although tensile instability is first found for SPH elastic dynamics,
its generally effective remedies are first found for SPH fluid dynamics [21,
22, 23]. Further analysis [20] suggested that tensile stability in SPH fluid
dynamics is highly relevant to zero-order consistency error, very often due to
the non-regular particle distributions which are typically generated by the
complex velocity gradient in flow field.

Such observations in flow simulations lead to a puzzle for ULSPH in solid
dynamics, where the velocity gradient is much more regular than that of a flow
field and theoretically should only generate very regular particle
distributions hence no tensile instability, except when very large deformation
or material failure happens. On the other hand, it is well known that TLSPH
elastic dynamics does not suffer from tensile instability as the ULSPH
counterpart. Such property, however, can be well explained by the non-updated
configuration obtained from the very regular initial particle distribution.

Figure 1: Illustration for (a) non-tensile-instability and non-hourglass
modes, (b) tensile instability and (c) hourglass modes in ULSPH simulations of
2D oscillating plates. The particles are colored with von Mises stress.

In this work, we propose to address these above-mentioned two drawbacks of
ULSPH in elastic dynamics. Firstly, we solve the puzzle why the supposed very
regular velocity field of elastic dynamics still suffers tensile instability.
For this, other than tension, we relate such issue with the hourglass modes
and shear stress, as shown in Fig. 1c, a numerical instability producing
”zigzag” particle distribution initially found in TLSPH elastic dynamics when
the deformation is very large [24]. Specifically, we argue that, at least for
the elastic dynamics without material failure, the previous consensus of
”tensile instability” in ULSPH is actually not caused by tension but hourglass
modes introduced by the nested formulation on the acceleration from shear
stress or the collocation of deformation and stress at the same particle
positions, just like it has been recently found in TLSPH [25].

Based on this argument, we develop a non-nested angular-momentum conservative
ULSPH formulation for computing the shear stress induced acceleration and
essentially eliminate the hourglass modes. We show that the new formulation is
tensile stable, i.e. without tensile instability, even when very large tension
and deformation are involved. Different from Ref. [5], the present solution
works for both 2D and 3D scenarios without resorting to case-dependent tuning.

Secondly, to improve computational efficiency, a dual-criteria time stepping
method [26] is incorporated into ULSPH simulations of elastic dynamics for the
first time. There are two time steps named the advection time step and the
much smaller acoustic time step. By updating of particle configurations only
in the advection time step, the frequency for updating particle configurations
can be minimized, and the calculation time can be reduced. It is worth
mentioning that the speed of sound in solid simulations is the true sound
speed, much higher than the artificial sound speed used in the WCSPH method
for fluids. This implies the dual-criteria time stepping scheme leads to a
significant enhancement in computational efficiency for solid simulations.

The remainder of this article is arranged as follows. The basic theory of
elastic dynamics is introduced in Section 2. The original formulation and the
present essentially non-hourglass and non-tensile-instability formulation for
elastic dynamics are described in Section 3 and Section 4 respectively. The
dual-criteria time stepping scheme for solid simulations is described in
Section 5, and a set of benchmark cases for elastic dynamics are then used to
validate the convergence, accuracy and stability of the proposed method in
Section 6. Section 7 draws the conclusion. For future in-depth research, all
the code used in this study has been open-sourced in the SPHinXsys repository
[27] at https://www.sphinxsys.org and https://github.com/Xiangyu-Hu/SPHinXsys.

##  2 Governing equations and constitutive relations

In a Lagrangian framework, the governing equations include mass and momentum
conservation for continuum mechanics are defined as

$\frac{\text{d}\rho}{\text{d}t}=-\mathbf{\rho}\nabla\cdot\mathbf{v}$ (1)
$\frac{\text{d}\mathbf{v}}{\text{d}t}=\frac{1}{\mathbf{\rho}}\nabla\cdot\bm{\sigma}+\mathbf{g}$
(2)

where ${\rho}$ is the density, $\mathbf{v}$ is velocity, ${t}$ is the time,
$\bm{\sigma}$ is the stress tensor, and $\mathbf{g}$ is the body force. The
total stress tensor $\bm{\sigma}$ can be divided into two terms, i.e., the
hydrostatic pressure and the shear stress, as shown below.

$\bm{\sigma}=-p\mathbf{I}+\bm{\sigma}^{s}$ (3)

where ${p}$ is the hydrostatic pressure, $\mathbf{I}$ is the identity matrix,
and $\bm{\sigma}^{s}$ is the shear stress. The pressure ${p}$ can be evaluated
from density based on an artificial equation of state [5].

$p=c_{0}^{2}(\rho-\rho_{0})$ (4)

where ${\rho_{0}}$ and ${\rho}$ are the initial and the current density
respectively. ${c_{0}}$ is the sound speed, which is expressed as [18]

$c_{0}=\sqrt{\frac{E}{3(1-2\nu)\rho_{0}}}$ (5)

where ${E}$ is the Young’s modulus, ${\nu}$ is the Poisson’s ratio of the
given material. The shear stress is the integral of the shear stress rate with
respect to time.

$\bm{\sigma}^{s}=\int_{0}^{t}\dot{\bm{\sigma}}^{s}\text{d}t$ (6)

For a linear elastic model, the shear stress rate is defined as

$\dot{\bm{\sigma}}^{s}=2G(\dot{\bm{\varepsilon}}-\frac{1}{d}tr(\dot{\bm{\varepsilon}})\mathbf{I})$
(7)

where G is the shear modulus. ${tr(\varphi)}$ indicates the trace of a
variable ${\varphi}$ and ${\dot{\varphi}}$ is the change rate with time for
the variable ${\varphi}$ (${\varphi}$ is an arbitrary variable). ${d}$
represents the space dimension, and ${d=2}$ and ${3}$ for 2D and 3D cases
respectively. ${\dot{\bm{\varepsilon}}}$ is strain rate, which is defined as

$\dot{\bm{\varepsilon}}=\frac{1}{2}(\nabla\mathbf{v}+(\nabla\mathbf{v})^{T})$
(8)

where ${\nabla\mathbf{v}}$ donates the velocity gradient, and superscript $T$
indicates the transpose of a tensor.

##  3 Original SPH formulation

According to Eq. (2) and Eq. (3), The acceleration related to volumetric
(hydrostatic pressure) and deviatoric part (shear stress) of the stress tensor
can be express as

$\dot{\mathbf{v}}^{p}=-\frac{1}{\mathbf{\rho}}\nabla p$ (9)
$\dot{\mathbf{v}}^{s}=\frac{1}{\mathbf{\rho}}\nabla\cdot\bm{\sigma}^{s}$ (10)

where $\dot{\mathbf{v}}^{p}$ and $\dot{\mathbf{v}}^{s}$ donate the velocity
change rate (acceleration) induced by hydrostatic pressure and shear stress
respectively. Then the total velocity change rate
$\dot{\mathbf{v}}$=$\dot{\mathbf{v}}^{p}$+$\dot{\mathbf{v}}^{s}$+$\mathbf{g}$.

A low-dissipation Riemann solver [28] is incorporated in the WCSPH to discrete
the continuity equation and the momentum equation for hydrostatic pressure.

$\frac{\text{d}\rho_{i}}{\text{d}t}=2\rho_{i}\sum_{j}\frac{m_{j}}{\rho_{j}}(U^{*}-\mathbf{v}_{ij}\mathbf{e}_{ij})\frac{\partial
W_{ij}}{\partial{r}_{ij}}$ (11)
$\frac{\text{d}\mathbf{v}_{i}^{p}}{\text{d}t}=-2\sum_{j}m_{j}\frac{P^{*}}{\rho_{i}\rho_{j}}{\nabla_{i}W_{ij}}$
(12)

Here, $W_{ij}$ represents $W({\mathbf{r}}_{i}-{\mathbf{r}}_{j},h)$, which is
the kernel function. ${\mathbf{r}}$ is particle position and ${h}$ is the
smoothing length. The subscripts ${i}$ and ${j}$ donate particle numbers, and
$m$ is the particle mass. $\mathbf{e}_{ij}$ is the unit vector pointing from
particle ${j}$ to particle ${i}$ and
$\mathbf{v}_{ij}=\mathbf{v}_{i}-\mathbf{v}_{j}$.
${\nabla_{i}W_{ij}}=\frac{\partial
W({r}_{ij},h)}{\partial{r}_{ij}}\mathbf{e}_{ij}$ is the derivative of kernel
function, and ${r}_{ij}=|\mathbf{r}_{i}-\mathbf{r}_{j}|$ is the distance
between two particles. $U^{*}$ and $P^{*}$, which are obtained from the low-
dissipation Riemann solver [27, 28], are the solutions of an inter-particle
Riemann problem along the unit vector pointing from particle $i$ to $j$.

Then the shear acceleration $\dot{\mathbf{v}}^{s}$ can be discretized by

$\frac{\text{d}\mathbf{v}_{i}^{s}}{\text{d}t}=\sum_{j}m_{j}\frac{\bm{\sigma}^{s}_{i}+\bm{\sigma}^{s}_{j}}{\rho_{i}\rho_{j}}\cdot{\nabla_{i}W_{ij}}$
(13)

Refer to Eq. (6)-Eq. (8), the velocity gradient needs to be discretized for
computing shear stress $\bm{\sigma}^{s}_{i}$ and $\bm{\sigma}^{s}_{j}$, which
has the following form [29]

${\nabla\mathbf{v}}=\sum_{j}\mathbf{v}_{ij}{\nabla_{i}W_{ij}}V_{j}$ (14)

where $V_{j}$ is the volume of particle $j$.

Combined Eq. (6)-Eq. (8) with Eq. (10), it is found that when calculating
shear acceleration $\dot{\mathbf{v}}^{s}$ within a time step, we first compute
velocity gradient ${\nabla\mathbf{v}}$ with Eq. (14), then update shear strain
rate ${\dot{\bm{\varepsilon}}}$, followed by computing shear stress rate
${\dot{\bm{\sigma}}^{s}}$, and subsequently calculate shear stress
$\bm{\sigma}^{s}$. Finally, the shear acceleration can be estimated by the
divergence of shear stress $\bm{\sigma}^{s}$ with Eq. (13). This nested
formulation of shear acceleration will trigger hourglass modes, which is the
origin of tensile instability in SPH simulation.

##  4 Essentially non-hourglass SPH formulation

The formulation in section 3 suffers from hourglass issues induced by zero
energy modes [30], which is characterized by a pattern of particle
displacement that does not correspond to rigid body motion but still results
in zero strain energy [30, 31, 32]. In a recent effective remedy for hourglass
model in TLSPH [25], the particle acceleration due to the divergence of shear
stress is directly obtained from a one-step Laplacian formulation of the
particle displacement other than the nested implementation of the 2nd-order
derivatives used in the original TLSPH. Actually, such non-nested SPH
formulation of Laplacian is widely used SPH fluid dynamics for computing the
viscous-force term in the Navier-Stokes equations [3, 4, 33], and is found
much stabler than the nested counterpart. Inspired by these previous
solutions, we apply a Laplacian operator to calculate the shear acceleration
directly in ULSPH simulations of elastic dynamics.

Firstly, we re-formulated the theoretical expression for shear acceleration by
combining Eq. (6)-Eq. (10)

$\dot{\mathbf{v}}^{s}=\frac{G}{\mathbf{\rho}}\nabla\cdot{\int_{0}^{t}{(\nabla\mathbf{v}+(\nabla\mathbf{v})^{T}-\frac{2}{d}\nabla\cdot\mathbf{v}\mathbf{I})}\text{d}t}$
(15)

By substituting the first Hamiltonian operator on the right-hand side into the
integral symbol, the following equation can be obtained

$\dot{\mathbf{v}}^{s}=\frac{G}{\mathbf{\rho}}{\int_{0}^{t}{(\nabla^{2}\mathbf{v}+\nabla\nabla\cdot\mathbf{v}-\frac{2}{d}\nabla\nabla\cdot\mathbf{v})}\text{d}t}$
(16)

When considering weakly compressible conditions, i.e.,
${\nabla\cdot\mathbf{v}\approx}$ 0, Eq. (16) can be simplified to

${\dot{\mathbf{v}}^{s}}\approx{\frac{G}{\mathbf{\rho}}{\int_{0}^{t}{(\nabla^{2}\mathbf{v})}\text{d}t}}$
(17)

The Laplacian operator needs to be discretized for calculating the shear
acceleration. Refer to the literature [29], $\nabla^{2}\mathbf{v}$ can be
discretized as

${\nabla^{2}\mathbf{v}}=2\sum_{j}\frac{\mathbf{v}_{ij}}{r_{ij}}\frac{\partial
W_{ij}}{\partial{r}_{ij}}V_{j}$ (18)

However, the aforementioned formula does not fulfill the requirement for
angular-momentum conservation. In other words, using this discretized form of
the second derivative of velocity to compute shear acceleration in Eq. (17)
cannot remove the effects of rigid rotation. Based on the research of Hu et
al. [34], we consider an angular-momentum conservative form for the second
derivative of velocity, to eliminate the contribution of rigid rotation to the
shear acceleration that should not have occurred in the first place, as shown
in Eq. (19).

${\nabla^{2}\mathbf{v}}=2\zeta{\sum_{j}\frac{\mathbf{e}_{ij}\cdot\mathbf{v}_{ij}}{r_{ij}}{\nabla_{i}W_{ij}}V_{j}}$
(19)

where $\zeta$ is a parameter related to the smoothing length $h$ and the type
of kernel function. $\zeta$ needs to be calibrated with numerical experiments.
Then the non-nested formulation of shear acceleration can be written as

${\dot{\mathbf{v}}^{s}}=2\zeta{\frac{G}{\mathbf{\rho}}{\int_{0}^{t}{\left(\sum_{j}\frac{\mathbf{e}_{ij}\cdot\mathbf{v}_{ij}}{r_{ij}}{\nabla_{i}W_{ij}}V_{j}\right)}\text{d}t}}$
(20)

This is the final form of shear acceleration without hourglass modes, and thus
can eliminate the tensile instability in essence. In Section 6.1, we will
provide a demonstration of the disparities between angular-momentum
conservative and non-conservative approaches in numerical computations. It
should be noted that the shear acceleration calculated at time step ${n}$ will
be used at the next time step ${n+1}$.

##  5 Dual-criteria time stepping

As we mentioned in Section 1, due to the necessity of updating particle
configurations at each computational time step, a persistently challenging
issue in ULSPH is its low computational efficiency. In this section, the dual-
criteria time stepping originally proposed for fluid simulations [26] is
introduced to solid simulations for the first time, to improve the calculation
efficiency by reducing the frequency for updating particle configurations,
while maintain high computation accuracy at the same time.

The dual-criteria time stepping strategy employs a larger advection time step
$\bigtriangleup t_{ad}$, and a smaller acoustic time step $\bigtriangleup
t_{ac}$. The particle configuration is updated in the advection time step
$\bigtriangleup t_{ad}$, which is defined as

$\bigtriangleup
t_{ad}=CFL_{ad}\frac{h}{\left\lvert\mathbf{v}\right\rvert_{max}}$ (21)

where $CFL_{ad}=0.2$, $\left\lvert\mathbf{v}\right\rvert_{max}$ is the maximum
particle advection speed and $h$ is the smoothing length. The acoustic time
step $\bigtriangleup t_{ac}$, involving the update of particle properties such
velocity and density, has the following form

$\bigtriangleup
t_{ac}=CFL_{ac}\frac{h}{c_{0}+\left\lvert\mathbf{v}\right\rvert_{max}}$ (22)

where $CFL_{ac}=0.4$ and $c_{0}$ is the sound speed.

Then the position-based Verlet scheme is applied for the acoustic time
integration [35]. The beginning of the acoustic time step is indicated by
superscript $n$, and the midpoint and new time step are donated by superscript
$n+\frac{1}{2}$ and $n+1$ respectively. In the Verlet scheme, the particle
position and density are firstly updated to the midpoint with

$\begin{cases}{\mathbf{r}}^{n+\frac{1}{2}}={\mathbf{r}}^{n}+\frac{1}{2}{\bigtriangleup
t_{ac}}{\mathbf{v}}^{n}\\\
{\rho}^{n+\frac{1}{2}}={\rho}^{n}+\frac{1}{2}{\bigtriangleup
t_{ac}}({\frac{\text{d}\rho}{\text{d}t}})^{n}\\\ \end{cases}$ (23)

Then the velocity is updated to the new time step after the particle
acceleration is determined.

${\mathbf{v}}_{n+1}={\mathbf{v}}_{n}+{\bigtriangleup
t_{ac}}({\frac{\text{d}\mathbf{v}}{\text{d}t}})^{n}$ (24)

Finally, the particle position and density are updated to the new time step by

$\begin{cases}{\mathbf{r}}^{n+1}={\mathbf{r}}^{n+\frac{1}{2}}+\frac{1}{2}{\bigtriangleup
t_{ac}}{\mathbf{v}}^{n+1}\\\
{\rho}^{n+1}={\rho}^{n+\frac{1}{2}}+\frac{1}{2}{\bigtriangleup
t_{ac}}({\frac{\text{d}\rho}{\text{d}t}})^{n+1}\\\ \end{cases}$ (25)

##  6 Numerical examples

In this section, several benchmark cases are tested, and the results are
compared with analytical solutions or the results from other numerical
studies, from both qualitative and quantitative perspectives. Specifically,
our results are compared with the results obtained using original nested
formulation and nested formulation with artificial stress [5], to illustrate
the calculation accuracy, stability, and robustness of the present method.
Several abbreviations for different SPH methods are defined to facilitate the
writing and reading of the article. The original SPH method is defined as
”SPH-OG”; the original SPH method with introducing artificial stress [5] is
represented by ”SPH-OAS”; the present essentially non-hourglass formulation is
expressed as ”SPH-ENOG” in the following text. There are two artificial
coefficients in the artificial stress term, and the selection of them refers
to the literature [5] in this section.

The 5th-order Wendland kernel [36] with a smoothing length of $h$ = 1.3$dp$
and the cut-off radius equals to 2.6$dp$, where $dp$ is the initial particle
spacing, is applied for all the cases in this paper. Based on the selected
kernel type and smoothing length, the coefficient $\zeta$ in Eq. (20) is set
to $0.7d$+2.1 (3.5 for 2D situations and 4.2 for 3D situations) according to
numerical experiments and is adopted throughout this study. All the physical
quantities in this article are presented in dimensionless form.

In this newly-developed non-hourglass formulation, we skip the calculation of
shear stress and the shear acceleration is computed directly. The shear stress
can be obtained separately by Eq. (6)-Eq. (8) if needed. The velocity gradient
in Eq. (8) is discretized by

${\nabla\mathbf{v}}=\sum_{j}\mathbf{v}_{ij}(\mathbf{B}_{i}{\nabla_{i}W_{ij}})V_{j}$
(26)

where $\mathbf{B}_{i}$ is the correction matrix for kernel gradient [12, 37,
38] and is defined as

$\mathbf{B}_{i}=\left({\sum_{j}\mathbf{r}_{ij}{\nabla_{i}W_{ij}}V_{j}}\right)^{-1}$
(27)

All the simulations in this section are run on a CentOS-8 system with 32
cores. The detail information of the CPU is ”64 Intel(R) Xeon(R) Gold 6226R
CPU @ 2.90GHz”.

###  6.1 2D oscillating plate

As shown in Fig. 2, a 2D plate with one edge fixed is firstly used to verify
the proposed method, and the results are compared with previous theoretical
[39] and numerical [5] solutions. The length and thickness of the plate are
$L$ and $H$ respectively, and the left part is fixed to produce a cantilever
plate. An observation point is set at the middle of the tail, and the vertical
displacement of the point is recorded as the amplitude. An initial velocity
$v_{y}$, which perpendicular to the plate strip, is applied with

$v_{y}(x)=v_{f}c_{0}\frac{f(x)}{f(L)}$ (28)

where the constant $v_{f}$ is an input parameter and

$\displaystyle f(x)=(\sin(kL)+\sinh(kL))(\cos(kx)-\cosh(kx))$ (29)
$\displaystyle-(\cos(kL)+\cosh(kL))(\sin(kx)-\sinh(kx))$

where $kL=1.875$ is determined by $\cos(kL)\cosh(kL)=-1$. The frequency
$\omega$ of the oscillating plate is theoretically given by

${\omega}^{2}=\frac{EH^{2}k^{4}}{12\rho_{0}(1-{\nu}^{4})}$ (30)

Figure 2: Model setup for 2D oscillating plate

The material and dimensional parameters in this case follow literatures [5,
18], i.e., density $\rho_{0}=1000$, Young’s modulus $E=2\times 10^{6}$,
Poisson’s ratio $\nu=0.3975$, $L=0.2$, and $H=0.02$.

As shown in Fig. 3(a), SPH-OG leads to hourglass and tensile instability
issues when simulating elastic deformation. Numerical fractures occur at the
beginning of the simulation ($t=0.05$), which is the tensile instability
addressed by previous researchers [5, 11]; the nonphysical zigzag particle
distribution and the non-uniform profile of von Mises stress indicate the
hourglass mode. Fig. 3(b) shows the results obtained by SPH-OAS, in which the
tensile instability can be suppressed. However, the hourglass still occurs and
became visually evident when $t=0.37$. This is because the error in the nested
formulation is in integral form, which gradually accumulates over time. The
results produced by the present SPH-ENOG are shown in Fig. 3(c). Clearly,
neither hourglass nor tensile instability appears even when the time $t=0.67$.
The particle distribution is still uniform, and the stress profile is smooth.

Figure 3: Evolution of particle configuration with time (t=0.05, 0.37 and
0.67) for (a) SPH-OG, (b) SPH-OAS, and (c) SPH-ENOG. The particles are colored
by von Mises stress.

The convergence of the present new formulation is validated, as shown in Fig.
4. Three cases with different resolutions ($H/dp=10$, $H/dp=20$ and $H/dp=30$)
are tested and the variation of amplitudes over time are illustrated. It can
be observed that, with the increase of resolution, the differences between
different solutions are deceasing, which is consistent with the results in
literatures [5, 18, 25] and indicates the convergence of this present
algorithm.

Figure 4: Convergence verification for the 2D oscillating plate with the
present SPH-ENOG. The curves show the changes of amplitude with time. Here,
$L$=0.2, $H$=0.02, and $\mathbf{v}_{f}$=0.05.

Then a stress testing is performed with a long-time simulation, to check the
stability of the current algorithm. As shown in Fig. 5, the simulation lasts
for over 30 oscillations, and the result from SPH-OAS is also illustrated for
comparison. Here, the two simulations are run with single time step [18], as
we want to minimize the accumulated integration error in long-time
simulations. The images above and below the curves in Fig. 5 respectively
represent the particle distribution obtained using SPH-ENOG and SPH-OAS at a
time around 10. As depicted in Fig. 5, with the proposed SPH-ENOG, the
particle and stress distribution still keep uniform until the end of the
simulation; while for SPH-OAS, the hourglass issue appears to be serious at
$t\approx 10$. Moreover, with the present SPH-ENOG, the amplitude only
decrease marginally at $t$ around 10 compared with the amplitude at $t=0$. The
slight decrease in amplitude over time is due to the numerical dissipation
introduced in the Riemann solver [28]. On the contrary, the SPH-OAS exhibits
rapid energy decay, thus it cannot be used for long-duration computations.

Figure 5: Test the long-term stability for the SPH-ENOG, the result is
compared with the SPH-OAS. Here, $L$=0.2, $H$=0.02, $H/dp=10$ and
$\mathbf{v}_{f}$=0.05. The particles are colored by von Mises stress.

Furthermore, the accuracy is checked and verified with theoretical solutions
and SPH-OAS from Gray et al.’s work [5]. As shown in Table 1, the first period
of oscillation is recorded for both SPH-ENOG and SPH-OAS with $L$=0.2,
$H$=0.02 and $H/dp=30$. Compared with the analytical solution, the errors of
the SPH-ENOG are at the same level as SPH-OAG, confirming the accuracy of the
current SPH-ENOG. It should be noted that the analytical solution are obtained
based on a thin plate model. If the thickness of the plate is reduced to
$H$=0.01, the error with $v_{f}$=0.001 for the SPH-ENOG also decrease to
around 0.6$\%$, which is in consistent with previous studies [5, 18].

Table 1: Comparison of oscillation periods $T$ obtained from the present SPH-ENOG, SPH-OAS and analytical solutions. Here, $L$=0.2, $H$=0.02 and $H/dp=30$. $v_{f}$ | 0.001 | 0.01 | 0.03 | 0.05
---|---|---|---|---
$T$ (Analytical) | 0.254
$T$ (SPH-ENOG) | 0.262 | 0.263 | 0.268 | 0.279
$T$ (SPH-OAS) | 0.273 | 0.273 | 0.275 | 0.278
Error (SPH-ENOG) |  3.1$\%$ |  3.5$\%$ |  5.5$\%$ |  9.8$\%$
Error (SPH-OAS) |  7.5$\%$ |  7.5$\%$ |  8.3$\%$ |  9.4$\%$

Simultaneously, we provided visual evidence to demonstrate that the tensile
instability occurs after hourglass modes. As shown in Fig. 6, when we simulate
the 2D oscillating plate using the SPH-OG, we first observe the appearance of
hourglass modes at the initial stage ($t=0.012$), and tension instability
gradually emerges afterwards at $t=0.021$.

Figure 6: Illustration for occurrences of hourglass modes and tension
instability in 2D oscillating plates with the SPH-OG. The hourglass issue
appears at the initial stages when (a) $t=0.012$, and the tensile instability
can be observed later at (b) $t=0.021$. The particles are colored by von Mises
stress.

Next, we present the effects of employing angular-momentum conservation (Eq.
(19)) and non-conservation (Eq. (18)) approaches in discretizing the second
derivative of velocity gradient and their implications on the results. As
shown in Fig. 7, it can be observed that when applying the non-conservative
form, the behavior of the oscillating plate does not oscillate as expected,
but exhibits a strong resistance to motion. This is because, the influence of
rigid rotation is not eliminated in the non-conservative form, which means
that rigid rotation can also produce significant shear forces, thereby
hindering the plate’s motion. Correspondingly, the behavior of the plate can
be correctly estimated with the angular-momentum conservative form,
demonstrating the necessity of adopting this conservative type.

Figure 7: Demonstration of the differences in the behavior of oscillating
plates when using angular-momentum conservative and non-conservative forms to
discretize the second derivative of velocity gradient: (a) Angular-momentum
conservative form; (b) Angular-momentum non-conservative form. Here, $L$=0.2,
$H$=0.02, $H/dp=10$ and $\mathbf{v}_{f}$=0.05. The particles are colored by
von Mises stress.

The performance of dual-criteria time stepping scheme is tested with the SPH-
ENOG, and the results from previous single time stepping method [18] is used
for comparison. As shown in Table 2, the simulation ends at physical time
$t=1$, and the wall-clock time spent for dual-criteria (donate as $T_{d}$) and
single-criteria (donate as $T_{s}$) time stepping is recorded. It can be seen
the computing efficiency of the dual-criteria time steeping is approximately
twice that of the previous method.

Table 2: Comparison of computational efficiency for dual-criteria ($T_{d}$) and single-criteria ($T_{s}$) time stepping scheme with 2D oscillating beams. Here, $L$=0.2, $H$=0.02 and $\mathbf{v}_{f}$=0.05. The simulation stops at physical time $t=1$ and the wall-clock time spent is recorded. $N_{p}$ represents the total particle number. $H/dp$ | 10 | 20 | 30
---|---|---|---
$N_{p}$ (k) | 1.6 | 5.8 | 12.6
$T_{d}$ | 25.8 | 66.3 | 130.1
$T_{s}$ | 43.0 | 139.5 | 333.0

###  6.2 3D oscillating plate

The oscillation of a 3D thin plate, with a simple support boundary condition
for all lateral edges, is considered in this section. As shown in Fig. 8, a
square plate with length $L=0.4$, width $=0.4$ and height $H=0.01$ is
constructed [8, 25, 40]. A simple boundary condition is applied to the
particles in the middle of the four lateral sides. Specifically, the
displacement of these particles along the $z$-direction is fixed. An
observation point is set at the center of the plate. The particles are
subjected to the initial velocity $v_{z}$

${v_{z}}(x,y)=\sin{\frac{m\pi x}{L}}\sin{\frac{n\pi y}{W}}$ (31)

where $m$ and $n$ donate integers controlling the vibration mode in $x$ and
$y$ directions respectively. The theoretical vibration period for the 3D thin
plate is given by

$T=\frac{2}{\pi}\left[\left({\frac{m}{L}}\right)^{2}+\left({\frac{n}{W}}\right)^{2}\right]^{-1}{\sqrt{\frac{\rho_{0}H}{D}}}$
(32)

where $D$ represents the flexural rigidity and is defined as

$D=\frac{EH^{3}}{12(1-\nu^{2})}$ (33)

The material parameters are set as follows: density $\rho_{0}=1000$, Young’s
modulus $E=1\times 10^{8}$, and Poisson’s ratio $\nu=0.3$.

Figure 8: Model setup for the 3D oscillating plate, and the observation point
locates at the center of the plate.

Firstly, the particle distribution and the profile of von Mises stress are
checked. Fig. 9 illustrates the deformed particle configuration with von Mises
stress contour obtained by the SPH-ENOG at time $t=0.01$ for vibration modes
$(m,n)=(1,1)$ and $(2,2)$. The SPH-ENOG can generate a smooth stress profile
without hourglass modes and tensile instability. Fig. 10 shows the evolution
of the amplitude (displacement along $z$-axis of the observation point) with
time for plates with different resolutions. As can be seen, with the increase
of resolution (decrease of initial particle spacing $dp$), the period and
amplitude of the curve gradually approach a certain value, indicating the
convergence of the SPH-ENOG. Moreover, the first periods of oscillations with
different $(m,n)$ values and resolutions are calculated and compared with
analytical solutions, to further validate the convergence and accuracy of the
present SPH-ENOG. As shown in Table 3, the periods converge rapidly with
increasing resolutions and agree well with the analytical solutions.

Figure 9: 3D oscillating plates obtained by the SPH-ENOG at time $t=0.01$ for vibration modes $(m,n)=(1,1)$ and $(2,2)$. The figures are colored by von Mises stress and the spatial particle discretization $H/dp=9$. Figure 10: Convergence verification for the 3D oscillating plate with the present SPH-ENOG. The curves show the changes of amplitude with time. Here, $(m,n)=(1,1)$. Table 3: Comparison of the first oscillation period $T$ obtained the present SPH-ENOG and analytical solutions. SPH-ENOG | $H/dp=3$ | $H/dp=5$ | $H/dp=9$ | Analytical
---|---|---|---|---
$(m,n)=(1,1)$ | 0.0572 | 0.0543 | 0.0529 | 0.0532
$(m,n)=(2,1)$ | 0.0231 | 0.0218 | 0.0212 | 0.0213
$(m,n)=(2,2)$ | 0.0148 | 0.0140 | 0.0136 | 0.0133

The performance of dual-criteria time stepping scheme is tested with the SPH-
ENOG for 3D oscillating plates, and the results from previous single time
stepping method is used for comparison. As shown in Table 4, the simulation
ends at physical time $t=0.1$, and the wall-clock time spent for dual-criteria
(donate as $T_{d}$) and single-criteria (donate as $T_{s}$) time stepping is
recorded. Clearly, when using the dual-criteria time stepping approach, the
computational time reduces to approximately one-third of the time required for
single-criteria time stepping.

Table 4: Comparison of computational efficiency for dual-criteria ($T_{d}$) and single-criteria ($T_{s}$) time stepping scheme with 3D oscillating beam. Here, $(m,n)=(2,2)$. The simulation stops at physical time $t=0.1$ and the wall-clock time spent is recorded. $N_{p}$ represents the total particle number. $H/dp$ | 3 | 5 | 9
---|---|---|---
$N_{p}$ (k) | 44.7 | 204.0 | 1179.4
$T_{d}$ | 122.6 | 981.3 | 10629.4
$T_{s}$ | 331.2 | 2683.1 | 30273.3

###  6.3 2D colliding rubber rings

The collision of two rubber rings is simulated in this section refer to
literatures [5, 10, 18]. As shown in Fig. 11, two rings with inner radius 0.03
and outer radius 0.04 are moving towards each other with the initial velocity
magnitude $v_{0}$ (the relative velocity of the two rings is 2$v_{0}$), and
the distance between the centers of the two rings is 0.09. The initial uniform
particle distribution is achieved by a level-set based pre-processing
technique [41]. When two rings collide with each other, a significant tensile
force will be generated. In this case, we will show that the numerical
fracture (i.e., tensile instability) produced by SPH-OG and the zigzag
particle/stress distribution (i.e., hourglass modes) produced by SPH-OG and
SPH-OAS do not occur in the present SPH-ENOG. The material parameters are set
as follows: density $\rho_{0}=1200$, Young’s modulus $E=1\times 10^{7}$, and
Poisson’s ratio $\nu=0.4$. The initial particle spacing is 0.001.

Figure 11: Model setup for 2D colliding rubber rings.

Fig. 12 shows the evolution of particle configuration for the SPH-OG, SPH-OAS
and the present SPH-ENOG when the initial velocity magnitude
$v_{0}=0.06c_{0}$. Clearly, the SPH-OG suffers from serious hourglass and
tensile instability at the beginning of the computation ($t=0.002$), and the
calculation process can barely continue. For the SPH-OAS, the tensile
instability can be suppressed, and the particle distribution is uniform at the
initial stage (t=0.002). However, with the passage of time, the zigzag
distribution of particle configuration and von Mises stress gradually becomes
apparent. While for the SPH-ENOG, the particle and stress distribution are
uniform during the whole calculation process, and the hourglass and tensile
instability issues can be completely removed.

Figure 12: Evolution of particle configuration with time ($t=0.002$, 0.005,
0.008 and 0.012) for 2D colliding rubber rings. The results are obtained by
different SPH methods, i.e., SPH-OG (left column), SPH-OAS (middle column),
and SPH-ENOG (right column). The initial velocity magnitude $v_{0}=0.06c_{0}$
and the figures are colored by von Mises stress.

Then we increase the initial velocity $v_{0}$ to $0.07c_{0}$, to test the
stability and robustness of the present SPH-ENOG. As shown in Fig. 13, the
hourglass and tensile instability become more significant for the SPH-OAS,
compared with Fig. 12 when $v_{0}=0.06c_{0}$. Particularly, after the
collision and rebound of the two rings ($t=0.012$), the zigzag distribution of
particles cannot be restored, indicating that this accumulated integration
error over time is significant. On the contrary, the present SPH-ENOG can
produce a uniform particle and stress distribution throughout the entire
process.

Figure 13: Evolution of particle configuration with time ($t=0.002$, 0.005,
0.008 and 0.012) for 2D colliding rubber rings. The results are obtained by
different SPH methods, i.e., SPH-OAS (left column) and SPH-ENOG (right
column). The initial velocity $v_{0}=0.07c_{0}$ and the figures are colored by
von Mises stress.

Furthermore, the initial velocity is increased to $v_{0}=0.08c_{0}$. It can be
seen from Fig. 14, not only the hourglass, but also the tensile instability
appears when $t\geqslant 0.005$ for the SPH-OAS. Fortunately, the present SPH-
ENOG performs well even at such large initial velocity, and the hourglass and
tensile instability can be perfectly eliminated, which suggests the stability
and robustness of the present SPH-ENOG.

Figure 14: Evolution of particle configuration with time ($t=0.002$, 0.005,
0.008 and 0.012) for 2D colliding rubber rings. The results are obtained by
different SPH methods, i.e., SPH-OAS (left column) and SPH-ENOG (right
column). The initial velocity $v_{0}=0.08c_{0}$ and the figures are colored by
von Mises stress.

###  6.4 3D colliding rubber balls

The 2D colliding rubber rings are extended to 3D to validate the proposed SPH-
ENOG for 3D scenarios. The initial setup follows Fig. 11, i.e., two hollow
rubber balls are moving towards each other with inner radius 0.03 and outer
radius 0.04. The initial distance between the centers of the two balls is 0.09
and the initial velocity magnitude for each ball is $v_{0}$. The initial
particle spacing $dp=0.001$ and the uniform particle distribution at the
beginning is realized by Yu et al.’s method [41]. The selection of material
parameters follows section 6.3.

Figure 15: Model setup for (a) 3D colliding hollow rubber balls and (b) the
half of each ball are showcased for proper visualization.

Fig. 16 shows the collision process of two balls at different times
($t=0.001$, 0.003, 0.005 and 0.007) with the present SPH-ENOG. The initial
velocity is set as $v_{0}=0.08c_{0}$. Half of each ball is displayed
separately here for proper visualization. Obviously, the distribution of von
Mises stress is smooth, and the particle configuration is uniform, which means
hourglass modes and tensile instability can be completely removed.

Figure 16: Evolution of particle configuration with time ($t=0.001$, 0.003,
0.005 and 0.007) for 3D colliding rubber balls. The results are obtained by
the present SPH-ENOG. The initial velocity magnitude $v_{0}=0.08c_{0}$ and the
figures are colored by von Mises stress.

###  6.5 2D rubber ball-plate interaction

Another interesting problem with large deformation and tension produced is
simulated in this section. As shown in Fig. 17a, a rubber ball with a radius
of 0.05 is used to impact a rubber plate [18]. The plate, fixed at both ends,
has dimensions of 0.5 in length and 0.02 in width. The center of the ball is
located 0.1 units away from the center of the plate, and the initial velocity
of the rubber ball is $v_{0}$. The material parameters of the rubber ball and
the target plate are the same [18], i.e., density $\rho_{0}=1200$, Young’s
modulus $E=1\times 10^{7}$, and Poisson’s ratio $\nu$ is set as 0.49 to
produce a large deformation, mimicking realistic rubber materials. The initial
particle spacing is 0.0025.

Figure 17: Model setup for (a) 2D rubber ball-plate interaction, and (b) 3D
rubber ball-plate interaction.

Fig. 18 illustrates the process of collision between the rubber ball and the
target plate with an initial velocity of 0.02$c_{0}$. Significant tensions
will be generated on the plate after the ball touches it. Serious hourglass
modes (zigzag particle distribution) and tensile instability (numerical
fractures) can be observed when using the SPH-OG (Fig. 18a). For the SPH-OAS
(Fig. 18b), the tensile instability can be suppressed but the hourglass issue
still exists. While with the present SPH-ENOG (Fig. 18c), the tensile
instability and hourglass modes can be removed simultaneously. The phenomena
are consistent with the results from the oscillating plate and colliding
rubber rings/balls described in Section 6.1 and Section 6.3.

Figure 18: 2D rubber ball-plate interaction with initial ball velocity
0.02$c_{0}$ at different times. The results are obtained by (a) SPH-OG, (b)
SPH-OAS and (c) SPH-ENOG respectively. The figures are colored by von Mises
stress.

Afterwards, we increase the initial velocity of the ball ($v_{0}=0.06c_{0}$)
to test the stability and applicability of the SPH-OAS and SPH-ENOG under a
more demanding condition. It can be seen that the hourglass issue still exists
for the SPH-OAS(Fig. 19a), and to make matters worse, the tensile instability
starts to appear at $t=0.012$ and become very severe later. However, our
method (Fig. 19b) still performs well under such condition, and both hourglass
modes and tension instability do not occur.

Figure 19: 2D rubber ball-plate interaction with initial ball velocity
0.06$c_{0}$ at different times. The results are obtained by (a) SPH-OAS and
(b) SPH-ENOG respectively. The figures are colored by von Mises stress.

To further challenge the proposed SPH-ENOG, the initial velocity is set to
$v_{0}=0.12c_{0}$ and extremely large tensions will be generated under such
condition. As shown in Fig. 20, surprisingly, the particle configuration still
keeps a uniform distribution and a smooth stress profile can be obtained,
implying the robustness and stability of the current SPH-ENOG.

Figure 20: 2D rubber ball-plate interaction with initial ball velocity
0.12$c_{0}$ at different times. The results are obtained by the present SPH-
ENOG and the figures are colored by von Mises stress.

###  6.6 3D rubber ball-plate interaction

We further consider the collision of the rubber ball and plate in 3D
situations, as shown in Fig. 17b. The plate with a size of $0.5\times
0.02\times 0.5$ is fixed on all four edges. The line connecting the center of
the ball and the center of the plate is perpendicular to the plane of the
plate, and the direction of initial velocity $v_{0}$ lies along this line. The
materials of the rubber ball and plate follows section 6.5, and the initial
particle spacing is 0.0025.

The results with the SPH-ENOG and $v_{0}=0.12c_{0}$ are shown in Fig. 21. The
profile of von Mises stress is smooth throughout the calculation process and
there are no occurrences of numerical fractures. This demonstrates the
capability of the present SPH-ENOG in eliminating hourglass modes and tensile
instability for 3D situations.

Figure 21: 3D rubber ball-plate interaction with initial ball velocity
0.12$c_{0}$ at different times. The results are obtained by the present SPH-
ENOG and the figures are colored by von Mises stress.

##  7 Conclusion remarks

This paper presents an essentially non-hourglass and non-tensile-instability
formulation for ULSPH by decomposing the shear acceleration into an angular-
momentum conservative form with the Laplacian operator. This newly-developed
method is applicable for both 2D and 3D scenarios without introducing case-
dependent tuning parameters. The convergence and accuracy of the present
method are verified through several fundamental test cases. Furthermore, the
stability and robustness of the non-hourglass and non-tensile-instability
formulation are confirmed through long-term simulations and simulations under
extreme conditions. Unlike previous methods that mitigate tension instability
through post-compensation and corrections, the current approach essentially
eliminates tension instability from a novel perspective, i.e., resolving the
hourglass mode in ULSPH.

This research also corrects a long-standing misconception regarding the causes
of tension instability, and clarifies its true origin. In the past, tension
instability was believed to be caused by tensions in elastic dynamics.
However, this research demonstrates that if hourglass modes are addressed at
its root, tension instability is naturally eliminated, even in scenarios
involving extremely large tensions (as shown in Fig. 20). This implies that
tension is not the true source of tension instability in elastic dynamics;
rather, it is the further development and exacerbation of hourglass modes in
the tensile region that causes it.

Last but not least, a dual-criteria time stepping scheme is employed to
increase the calculation efficiency. For 2D and 3D oscillating plates, the
computational time is reduced to approximately one-half and one-third,
respectively, compared to the original calculations using single-criteria time
stepping approach.

It should be noted that, although the present formulation is proposed and
validated for elastic dynamics, it is possible to extend the current method to
plastic dynamics for modelling material fracture and failure, as our framework
is developed in the ULSPH.

## CRediT authorship contribution statement

Shuaihao Zhang: Conceptualization, Methodology, Investigation, Visualization,
Validation, Formal analysis, Writing - original draft, Writing - review &
editing. Sérgio D.N. Lourenço: Supervision, Investigation, Writing - review &
editing. Dong Wu: Investigation, Methodology, Formal analysis, Writing -
review & editing. Chi Zhang: Methodology, Writing - review & editing. Xiangyu
Hu: Supervision, Investigation, Methodology, Writing - review & editing.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## Data availability

The code and data are available on GitHub.

## Acknowledgements

Sérgio D.N. Lourenço would like to express his gratitude to the Research
Grants Council Hong Kong for their sponsorship of this research under a
Collaborative Research Fund (C6006-20GF). Dong Wu, Chi Zhang and Xiangyu Hu
would like to express their gratitude to the German Research Foundation (DFG)
for their sponsorship of this research under grant number DFG HU1527/12-4. The
computations were performed using research computing facilities offered by
Information Technology Services, the University of Hong Kong.

## References

  * [1] L. B. Lucy, A numerical approach to the testing of the fission hypothesis, The astronomical journal 82 (1977) 1013–1024.
  * [2] R. A. Gingold, J. J. Monaghan, Smoothed particle hydrodynamics: Theory and application to non-spherical stars, Monthly notices of the royal astronomical society 181 (3) (1977) 375–389.
  * [3] J. P. Morris, P. J. Fox, Y. Zhu, Modeling low reynolds number incompressible flows using sph, Journal of computational physics 136 (1) (1997) 214–226.
  * [4] X. Y. Hu, N. A. Adams, A multi-phase sph method for macroscopic and mesoscopic flows, Journal of Computational Physics 213 (2) (2006) 844–861.
  * [5] J. P. Gray, J. J. Monaghan, R. Swift, Sph elastic dynamics, Computer methods in applied mechanics and engineering 190 (49-50) (2001) 6641–6662.
  * [6] G. R. Johnson, R. A. Stryk, S. R. Beissel, Sph for high velocity impact computations, Computer methods in applied mechanics and engineering 139 (1-4) (1996) 347–373.
  * [7] C. Antoci, M. Gallati, S. Sibilla, Numerical simulation of fluid–structure interaction by sph, Computers & structures 85 (11-14) (2007) 879–890.
  * [8] A. Khayyer, H. Gotoh, Y. Shimizu, Y. Nishijima, A 3d lagrangian meshfree projection-based solver for hydroelastic fluid-structure interactions, Journal of Fluids and Structures 105 (2021) 103342.
  * [9] R. Vignjevic, J. R. Reveles, J. Campbell, Sph in a total lagrangian formalism, CMC-Tech Science Press- 4 (3) (2006) 181.
  * [10] J. J. Monaghan, Sph without a tensile instability, Journal of computational physics 159 (2) (2000) 290–311.
  * [11] J. W. Swegle, D. L. Hicks, S. W. Attaway, Smoothed particle hydrodynamics stability analysis, Journal of computational physics 116 (1) (1995) 123–134.
  * [12] P. Randles, L. D. Libersky, Smoothed particle hydrodynamics: some recent improvements and applications, Computer methods in applied mechanics and engineering 139 (1-4) (1996) 375–408.
  * [13] G. R. Johnson, S. R. Beissel, Normalized smoothing functions for sph impact computations, International Journal for Numerical Methods in Engineering 39 (16) (1996) 2725–2741.
  * [14] C. Dyka, P. Randles, R. Ingel, Stress points for tension instability in sph, International Journal for Numerical Methods in Engineering 40 (13) (1997) 2325–2341.
  * [15] D. Mandell, C. Wingate, L. Schwalbe, Computational brittle fracture using smooth particle hydrodynamics, Tech. rep., Los Alamos National Lab.(LANL), Los Alamos, NM (United States) (1996).
  * [16] G. A. Dilts, Moving-least-squares-particle hydrodynamics—i. consistency and stability, International Journal for Numerical Methods in Engineering 44 (8) (1999) 1115–1155.
  * [17] P. Randles, L. Libersky, Normalized sph with stress points, International Journal for Numerical Methods in Engineering 48 (10) (2000) 1445–1462.
  * [18] C. Zhang, X. Y. Hu, N. A. Adams, A generalized transport-velocity formulation for smoothed particle hydrodynamics, Journal of Computational Physics 337 (2017) 216–232.
  * [19] L. Lobovskỳ, J. Křen, Smoothed particle hydrodynamics modelling of fluids and solids (2007).
  * [20] S. Litvinov, X. Hu, N. A. Adams, Towards consistence and convergence of conservative sph approximations, Journal of Computational Physics 301 (2015) 394–401.
  * [21] X. Hu, N. A. Adams, An incompressible multi-phase sph method, Journal of computational physics 227 (1) (2007) 264–278.
  * [22] R. Xu, P. Stansby, D. Laurence, Accuracy and stability in incompressible sph (isph) based on the projection method and a new approach, Journal of computational Physics 228 (18) (2009) 6703–6725.
  * [23] S. Adami, X. Hu, N. A. Adams, A transport-velocity formulation for smoothed particle hydrodynamics, Journal of Computational Physics 241 (2013) 292–307.
  * [24] G. C. Ganzenmüller, An hourglass control algorithm for lagrangian smooth particle hydrodynamics, Computer Methods in Applied Mechanics and Engineering 286 (2015) 87–106.
  * [25] D. Wu, C. Zhang, X. Tang, X. Hu, An essentially non-hourglass formulation for total lagrangian smoothed particle hydrodynamics, Computer Methods in Applied Mechanics and Engineering 407 (2023) 115915.
  * [26] C. Zhang, M. Rezavand, X. Hu, Dual-criteria time stepping for weakly compressible smoothed particle hydrodynamics, Journal of Computational Physics 404 (2020) 109135.
  * [27] C. Zhang, M. Rezavand, Y. Zhu, Y. Yu, D. Wu, W. Zhang, J. Wang, X. Hu, Sphinxsys: An open-source multi-physics and multi-resolution library based on smoothed particle hydrodynamics, Computer Physics Communications 267 (2021) 108066\.
  * [28] C. Zhang, X. Hu, N. A. Adams, A weakly compressible sph method based on a low-dissipation riemann solver, Journal of Computational Physics 335 (2017) 605–620.
  * [29] P. Espanol, M. Revenga, Smoothed dissipative particle dynamics, Physical Review E 67 (2) (2003) 026705.
  * [30] R. Vignjevic, J. Campbell, L. Libersky, A treatment of zero-energy modes in the smoothed particle hydrodynamics method, Computer methods in Applied mechanics and Engineering 184 (1) (2000) 67–85.
  * [31] R. Vignjevic, J. Campbell, Review of development of the smooth particle hydrodynamics (sph) method, in: Predictive Modeling of Dynamic Processes: A Tribute to Professor Klaus Thoma, Springer, 2009, pp. 367–396.
  * [32] J. Swegle, et al., An analysis of smoothed particle hydrodynamics, sandia report no, Tech. rep., SAND93-2513-UC-705 (1994).
  * [33] J. J. Monaghan, Smoothed particle hydrodynamics, Reports on Progress in Physics 68 (2005) 1703 – 1759.
  * [34] X. Hu, N. Adams, Angular-momentum conservative smoothed particle dynamics for incompressible viscous flows, Physics of Fluids 18 (10) (2006).
  * [35] C. Zhang, M. Rezavand, X. Hu, A multi-resolution sph method for fluid-structure interactions, Journal of Computational Physics 429 (2021) 110028.
  * [36] H. Wendland, Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree, Advances in computational Mathematics 4 (1995) 389–396.
  * [37] J. Bonet, S. Kulasegaram, A simplified approach to enhance the performance of smooth particle hydrodynamics methods, Applied Mathematics and Computation 126 (2-3) (2002) 133–155.
  * [38] C. Zhang, J. Wang, M. Rezavand, D. Wu, X. Hu, An integrative smoothed particle hydrodynamics method for modeling cardiac function, Computer Methods in Applied Mechanics and Engineering 381 (2021) 113847.
  * [39] L. D. Landau, E. M. Lifshitz, Course of theoretical physics, Elsevier, 2013.
  * [40] A. W. Leissa, Vibration of plates, Vol. 160, Scientific and Technical Information Division, National Aeronautics and Space Administration, 1969.
  * [41] Y. Yu, Y. Zhu, C. Zhang, O. J. Haidn, X. Hu, Level-set based pre-processing techniques for particle methods, Computer Physics Communications 289 (2023) 108744\.


# Correlation of Coronal Mass Ejection Shock Temperature with Solar Energetic
Particle Intensity

Manuel Enrique Cuesta Princeton University, Department of Astrophysical
Sciences, Princeton, NJ 08544, USA  D. J. McComas Princeton University,
Department of Astrophysical Sciences, Princeton, NJ 08544, USA  L. Y. Khoo
Princeton University, Department of Astrophysical Sciences, Princeton, NJ
08544, USA  R. Bandyopadhyay Princeton University, Department of Astrophysical
Sciences, Princeton, NJ 08544, USA  T. Sharma Princeton University, Department
of Astrophysical Sciences, Princeton, NJ 08544, USA  M. M. Shen Princeton
University, Department of Astrophysical Sciences, Princeton, NJ 08544, USA  J.
S. Rankin Princeton University, Department of Astrophysical Sciences,
Princeton, NJ 08544, USA  A. T. Cummings  Princeton University, Department of
Astrophysical Sciences, Princeton, NJ 08544, USA  J. R. Szalay Princeton
University, Department of Astrophysical Sciences, Princeton, NJ 08544, USA  C.
M. S. Cohen Space Research Lab, California Institute of Technology, Pasadena,
CA 91125, USA  N. A. Schwadron Princeton University, Department of
Astrophysical Sciences, Princeton, NJ 08544, USA  University of New Hampshire,
Durham, NH 03824, USA  R. Chhiber NASA Goddard Space Flight Center, Greenbelt,
MD 20771, USA  Department of Physics and Astronomy, Bartol Research Institute,
University of Delaware, Newark, DE 19716, USA  F. Pecora Department of Physics
and Astronomy, Bartol Research Institute, University of Delaware, Newark, DE
19716, USA  W. H. Matthaeus Department of Physics and Astronomy, Bartol
Research Institute, University of Delaware, Newark, DE 19716, USA  R. A. Leske
Space Research Lab, California Institute of Technology, Pasadena, CA 91125,
USA  M. L. Stevens Smithsonian Astrophysical Observatory, Cambridge, MA 02138,
USA

###### Abstract

Solar energetic particle (SEP) events have been observed by the Parker Solar
Probe (PSP) spacecraft since its launch in 2018. These events include sources
from solar flares and coronal mass ejections (CMEs). Onboard PSP is the
IS$\odot$IS instrument suite measuring ions over energies from $\sim$ 20
keV/nucleon to 200 MeV/nucleon and electrons from $\sim$ 20 keV to 6 MeV.
Previous studies sought to group CME characteristics based on their plasma
conditions and arrived at general descriptions with large statistical errors,
leaving open questions on how to properly group CMEs based solely on their
plasma conditions. To help resolve these open questions, plasma properties of
CMEs have been examined in relation to SEPs. Here we reexamine one plasma
property, the solar wind proton temperature, and compare it to the proton SEP
intensity in a region immediately downstream of a CME-driven shock for seven
CMEs observed at radial distances within 1 au. We find a statistically strong
correlation between proton SEP intensity and bulk proton temperature,
indicating a clear relationship between SEPs and the conditions in the solar
wind. Furthermore, we propose that an indirect coupling of SEP intensity to
the level of turbulence and the amount of energy dissipation that results is
mainly responsible for the observed correlation between SEP intensity and
proton temperature. These results are key to understanding the interaction of
SEPs with the bulk solar wind in CME-driven shocks and will improve our
ability to model the interplay of shock evolution and particle acceleration.

Solar energetic particles, Solar coronal mass ejection shocks, solar wind

††journal: ApJ

##  1 Introduction

Coronal mass ejections (CMEs) are solar events involving the expulsion of
large amounts of solar coronal material into interplanetary space. These
events are observed to have wide variability in their characteristics, such as
their composition, magnetic field strength and behavior, in-situ plasma
conditions, and shock strengths. Many of these characteristics have been
studied resulting in an array of expectations that vary from one CME to the
next (Cane et al., 1990; Gosling et al., 1991; Neugebauer et al., 1997;
Zurbuchen & Richardson, 2006; Richardson & Cane, 2010; Wu & Lepping, 2011; Chi
et al., 2016, and references therein). A recent investigation into the
turbulence properties of CME cores inside 1 au (Good et al., 2023) via Solar
Orbiter and Parker Solar Probe (PSP) observations revealed radial trends in
the magnetic turbulent energy and the inertial-to-dissipative range break-
point frequency, which is analogous to the dissipation scale where turbulent
energy cascaded from larger to smaller scales can dissipate as heat. These
authors further demonstrated the variability of other CME characteristics,
such as plasma beta (ratio of thermal to magnetic pressure) and correlation
length, leaving an open question on how to properly group CMEs based solely on
their plasma properties.

To better group CMEs, their plasma properties have been examined in relation
to solar energetic particles (SEPs). CMEs can drive an upstream forward
propagating shock if they are traveling faster than the upstream solar wind by
at least the local fast magnetosonic speed, where the upstream plasma is
unshocked and the downstream plasma is shocked. These CME-driven shocks become
a source for particle acceleration (Lee, 1997) with variable efficiency in
their acceleration processes (Kallenrode, 1996; Lario et al., 1998), where
turbulent diffusion is thought to produce the shock’s dominant particle
acceleration process (Fisk, 1971; Axford et al., 1977; Blandford & Ostriker,
1978; Bell, 1978; Gosling et al., 1979; Lee, 1983). Noticeable increases in
the suprathermal particle population (energies generally between $\sim$ 10 keV
to $\sim$ 1 MeV per nucleon (Mason & Gloeckler, 2012)) and of those with
higher energy at CME-driven shocks are the result of these local acceleration
processes (Reames, 1999). CME-driven shock speeds have been found to be
correlated with peak SEP flux at the shock (Cane et al., 1990; Reames, 2000;
Kahler, 2001) for higher energies ($>10$ MeV), with considerable scatter. This
correlation can be explained by the association of stronger shocks with more
intense CMEs resulting in a larger SEP flux of higher energy particles.
Nonetheless, the considerable scatter over several decades of SEP intensity
for a given range of CME speeds suggests the peak SEP flux is also dependent
on other factors besides the CME speed, such as the ambient SEP population
(Kahler, 2001). When considering the behavior of the peak energetic particle
(EP) intensity in the heliosheath, Zank et al. (2015) found that the EP
intensity peaks downstream of the heliospheric bow shock and the distance
between the shock and this peak increases with increasing particle energy.

One method used to characterize the relationship between SEPs and turbulence
is the partial variance of increments (PVI) (Greco et al., 2008), a measure of
turbulence intermittency or roughness of the magnetic field at a given
increment scale, which is found to be well-correlated with solar wind
temperature (Osman et al., 2011, 2012; Greco et al., 2012; Chasapis et al.,
2015; Yordanova et al., 2016; Qudsi et al., 2020). Additionally, Tessein et
al. (2013, 2015) and Bandyopadhyay et al. (2020) found a correlation between
PVI and SEP intensity, where others discovered an association between local
acceleration regions and magnetic-like islands formed due to turbulence
(Khabarova et al., 2015; Khabarova & Zank, 2017; Malandraki et al., 2019).
Although PVI and SEP intensity are correlated on average, there are instances
where SEP enhancements are observed at times with low PVI. These times with
larger SEP intensity, however, are bounded by high PVI, suggesting that
coherent structures play a role in the local modulation and trapping of
energetic particles (Tessein et al., 2015, 2016; Pecora et al., 2021).

Based on the outcomes of the observational studies above, one might expect SEP
intensity can be used to gauge the solar wind proton temperature, a connection
motivated by the relation between SEP intensity and shock strength. From the
perspective of theory and modeling, the influence of the upstream turbulence
level and presence of coherent structures are found to result in more
efficient diffusive-like spreading of particles in both position and velocity
space (Trotta et al., 2021). Furthermore, the level of upstream turbulence
induces shock-surface fluctuations and downstream turbulence structuring that
have been linked to temperature anisotropy downstream of the shock (Trotta et
al., 2023). However, there are no studies to date that have directly
investigated the relationship between the bulk proton temperature immediately
downstream of CME-driven shocks and the intensity of suprathermal particle
populations locally accelerated by turbulent diffusion near the shock
(although, see Dayeh et al. (2018) and the inferred results of Lario et al.
(2017)). The shock strength and turbulence strength near the shock are found
to be well-correlated, where turbulence is known to affect SEP transport
(Jokipii, 1982, 1987; Dröge, 1994; Guo et al., 2021). This natural interaction
between SEPs and turbulence near the shock may also be reflected in the proton
temperature, since intermittency results in increased heating of the bulk
material at the targeted scales.

In this paper, we examine the relationship of proton temperature and particle
intensity through the analysis of SEP and in-situ plasma measurements by PSP
(Fox et al., 2016) for seven CME-driven shocks. In Section 2, we explain the
data and methods used to reach the results in Section 3. We discuss the
results and their impact on the understanding of SEP interactions with the
bulk solar wind in shocks in Section 4. Finally, we summarize the results in
Section 5.

##  2 Data and Methods

We utilize publicly available measurements of solar energetic particles via
PSP/IS$\odot$IS EPI-Lo (McComas et al., 2016), magnetic field data via
PSP/FIELDS (Bale et al., 2016), and proton temperatures, bulk velocity, and
proton density via PSP/SWEAP/SPAN-I (Kasper et al., 2016) from NASA’s Space
Physics Data Facility111Space Physics Data Facility can be accessed at
https://cdaweb.gsfc.nasa.gov/.. We identify possible shock candidates from the
FIELDS and SWEAP data and identify the relevant one by comparing its time to
the estimated shock arrival time. Estimated shock arrival times are found by
first considering their associated type-II radio bursts (Bale et al., 1999;
Knock et al., 2001; Mann & Klassen, 2005), measured by PSP/FIELDS, adjusted to
solar origin, denoted by time $t_{0}$. The time of solar origin is determined
by subtracting the travel time of light from the Sun to the radial distance of
PSP from the time of the observed radio-II burst. Next, we find the time
associated with the propagation of the CME-driven shock (estimated by the WSA-
ENLIL model (Odstrčil et al., 1996; Odstrčil & Pizzo, 1999a, b; Odstrčil,
2003; Odstrčil et al., 2004) for each shock separately) starting from a solar
origin to the radial distance of PSP, denoted as time $t_{1}$. Then we select
the nearest shock candidate to the time $t_{0}+t_{1}$. This gives us a total
of seven CME-driven shocks, overviews of which can be found in the Appendix
(Figure A.4), with vertical dashed red lines used to mark the corresponding
shock arrivals.

To reinforce the connection between SEP intensity and proton temperature via
turbulent processes, we compute the PVI of the magnetic field downsampled to 1
second cadence from the native ‘4_Sa_per_Cyc’ data product. The PVI is defined
in Greco et al. (2008) as

${\rm
PVI}(t,\tau)=\frac{\left|\Delta\bm{B}(t,\tau)\right|}{\sqrt{\langle\left|\Delta\bm{B}(t,\tau)\right|^{2}\rangle}}$
(1)

where $\Delta\bm{B}(t,\tau)=\bm{B}(t+\tau)-\bm{B}(t)$ is the increment of
magnetic field vector $\bm{B}$ at time $t$ and temporal increment scale
$\tau$, and the brackets $\langle\cdot\rangle$ represent an arithmetic mean.
We choose an increment scale of 10 seconds, which is shorter than that used in
Bandyopadhyay et al. (2020) due to the availability of higher resolution
magnetic field data and is in the inertial range far from the correlation
scale but larger than the dissipative scale.

The duration of the arithmetic mean used in the normalization of PVI varies
from event to event. To determine this duration, we first compute the auto-
correlation function:

$R(\tau)=\frac{\langle\bm{b}(t)\cdot\bm{b}(t+\tau)\rangle_{t}}{\langle\bm{b}(t)\cdot\bm{b}(t)\rangle_{t}}$
(2)

where $\bm{b}(t)=\bm{B}(t)-\langle\bm{B}(t)\rangle_{t}$ is the fluctuating
component of the magnetic field and $\langle\cdot\rangle_{t}$ is computed via
the Germano approach (Germano, 1992). The interval duration used in Equation
(2) is set at 16 hours centered about the estimated shock arrival time, which
is an appropriate interval duration at 1 au (Isaacs et al., 2015). The value
of $\tau_{c}$ such that $R(\tau=\tau_{c})$ drops below $e^{-1}$ is the
“e-folding” time or correlation time (Matthaeus et al., 1999). We then set the
duration of the arithmetic mean used in the normalization of PVI to
$10\cdot\tau_{c}$ centered about the estimated shock arrival time.

However, in this study, PSP encounters each of the CME-driven shocks within 1
au where the correlation scale is expected to be shorter, on average, which
presumably calls for shorter averaging intervals when computing magnetic
correlations. Therefore, assuming a 16 hour duration will lead to an
overestimated correlation scale and a larger averaging duration used for
normalizing the PVI. Although this is important with investigations that
directly depend on the correlation scale, an averaging duration for the PVI
normalization beyond 10 times the actual correlation scale is found to
insignificantly affect PVI enhancements (Servidio et al., 2011). Once the
PVI($\tau=10~{}{\rm seconds}$) is obtained, we compute the fraction of PVI
values, denoted as PVIfrac, greater than 2.5 only within the duration used for
time-integration and averaging of the energy-integrated intensity and proton
temperature, respectively. This threshold is selected since coherent
structures are associated with PVI values greater than 2.5 (Greco et al.,
2008, 2012, 2018).

To compute the SEP intensity of protons shown below in Section 3, we examine
proton intensities measured by the EPI-Lo instrument. EPI-Lo is comprised of 8
wedges, each with 10 apertures having different look directions (Hill et al.,
2017). Here we utilize all apertures that are unaffected by dust puncture
events on EPI-Lo (Szalay et al. (2020); Shen et al. (2023, in prep.)) to
account for the changing magnetic field direction (or changing pitch angle,
the angle between the observed particle’s incident direction and the
instantaneous magnetic field direction) during CME shock passage. The data
product used in this study is Channel P, which corresponds to protons with
triple coincidence measurements (Hill et al., 2017).

We first integrate the measured flux ($F$) over the energy ($E$) domain (100
keV to 500 keV), normalized by the summation of the energy bin widths to
properly account for energy bin weights. This results in an energy-integrated
flux ($F_{E}(t)$) defined as

$F_{E}(t)=\frac{\sum^{500{\rm keV}}_{100{\rm keV}}\left[\Delta E_{i}\cdot
F(t,E_{i})\right]}{\sum^{500{\rm keV}}_{100{\rm keV}}\left[\Delta
E_{i}\right]}$ (3)

where $\Delta E_{i}$ is the energy bin width associated with the measured flux
at energy $E_{i}$. We then integrate $F_{E}(t)$ over the time domain for a
duration equivalent to $3\cdot\tau_{c}$, which is expected to vary from event
to event. We choose this duration to (1) constrain the results to describe the
downstream region nearest to the shock and (2) include magnetic structures
large enough to capture the acceleration of and/or systematic trapping of
energetic particles. For event-to-event comparisons, we normalize time
integration by the duration of the integration. We also confirmed that the
qualitative results and strength of correlation between any pair of parameters
(proton temperature, proton SEP intensity, and PVIfrac; see Appendix, Table 1)
are not sensitive to our exact choice of time window; changing the durations
of the window from 0.5 to $5\cdot\tau_{c}$ only insignificantly affects the
results.

This results in a time and integrated measured flux in units of
${\rm(cm^{2}\cdot sr\cdot s\cdot keV/nuc)^{-1}}$, defined as

$F_{E,t}={\rm Intensity}=\frac{\sum^{t_{s}+3\tau_{c}}_{t_{s}}\left[\Delta
t\cdot F_{E}(t)\right]}{\sum^{t_{s}+3\tau_{c}}_{t_{s}}\left[\Delta t\right]}$
(4)

where $t_{s}$ is the time of the CME-driven shock. For the rest of the paper,
we interchangeably use the terms ‘flux’ and ‘intensity’ to describe the value
$F_{E,t}$ where it does not cause confusion. We compare the resulting proton
SEP intensity to the proton temperature $T_{p}$ measured by SPAN-I from
PSP/SWEAP in its native time resolution, averaged over the same duration of
$3\cdot\tau_{c}$ downstream of the CME-driven shock.

One challenge is that instrument uncertainties in SPAN-I proton temperatures
are unavailable. Therefore, we assume an uncertainty of each measurement of
20% (Kasper et al., 2016), which is a reasonable estimate for the radial
distances at which PSP observed the CMEs ($\sim$ 0.35 au to .80 au) and the
SPAN-I field-of-view throughout the averaging duration. Quality flags in the
SPAN-I data product are checked to ensure reasonable confidence in the
measured value. Also, we impose one condition on the data for the SEP event to
be considered for analysis: data availability throughout the passing of a CME-
driven shock by all three instrument suites (IS$\odot$IS/FIELDS/SWEAP).
Several CME-driven shocks, beyond those included here, observed by
PSP/IS$\odot$IS were rendered unusable for the present analysis due to data
gaps from instrument turn-offs and spacecraft rolls for satellite
repositioning during the shock encounter.

##  3 Results

Here we compile the shock intensity and proton temperature for seven different
CMEs, with event overviews provided in the Appendix, Figure A.4. In Figure 1,
we provide an example of the analysis for the Feb. 15, 2022, CME event. The
top panel shows energy-integrated proton intensities ($F_{E}(t)$) over
energies 100 keV to 500 keV, the middle panel shows bulk proton temperatures,
and the bottom panel shows PVI at a 10 second increment. The vertical dashed
red line corresponds to the estimated shock arrival time, for which each curve
is adjusted. The duration between the estimated shock arrival time and the
vertical solid black line corresponds to the time window used to compute
$F_{E,T}$ and average proton temperature (in this case, $3\cdot\tau_{c}$).
Parameters of average proton temperature ($\langle T_{p}\rangle$), proton SEP
intensity, and PVIfrac over a $3\cdot\tau_{c}$ duration for the seven CME
events are given in the Appendix, Table 2.

Figure 1: Energy-integrated proton intensity over energies 100 keV to 500 keV
(top), solar wind proton temperature (middle), and PVI($\tau=10~{}{\rm
seconds}$) (bottom) centered about a $3\cdot\tau_{c}$ duration marked between
the shock arrival (vertical dashed red line) and the vertical dashed black
line associated with the Feb. 15, 2022, CME event.

The compiled values of the intensity and proton temperature are shown in the
left panel of Figure 2 for the selected duration of $3\times\tau_{c}$. The
horizontal error bars represent the propagated error throughout the time and
energy integrated intensity, and the vertical error bars represent the
propagated error from the proton temperatures. We utilize an orthogonal
distance regression (ODR; Boggs & Rogers, 1990) power-law fit between the
proton SEP intensity and $\langle T_{p}\rangle$, which considers uncertainty
in dependent and independent variables, only to show the proportionality
between the two parameters (left panel of Figure 2). The fitted slope and
standard error are given in the corresponding legend. In the middle and right
panels of Figure 2, we examine the behavior of proton SEP intensity and proton
temperature as a function of PVIfrac, respectively. An arrow is drawn in the
middle and right panels of Figure 2 only to further demonstrate the
relationship between PVI and either $\langle T_{p}\rangle$ or SEP intensity;
regions with proportionately enhanced PVI are associated with regions of
hotter bulk plasma and larger SEP intensity.

Figure 2: Compiled solar wind $\langle T_{p}\rangle$, proton SEP intensity,
and PVIfrac each as a function of the other (black circles) in the region
immediately downstream of the appropriate CME-driven shock. In the left panel,
vertical error bars are propagated standard errors from averaged temperatures
and horizontal error bars are propagated uncertainties from time and energy
integrated proton flux. The solid black line is an ODR power-law fit to show
proportionality. The correlation coefficients (r) of the parameters in the
left, middle, and right panel are $0.95$, $0.81$, and $0.81$, respectively.
Arrows in the middle and right panels are to express a general trend.

Figure 2 reveals correlations between SEP intensity, bulk proton temperature,
and PVIfrac downstream of the shock over a duration of $3\cdot\tau_{c}$. While
there is only a small statistical sample via PSP observations due to
constraints on data availability, the correlation coefficients ($r$) of
$\log\langle T_{p}\rangle$ to $\log{\rm Intensity}$, $\log{\rm Intensity}$ to
$\log{\rm PVI}_{frac}$, and $\log\langle T_{p}\rangle$ to $\log{\rm
PVI}_{frac}$ are 0.95, 0.81, and 0.81, respectively, which all correspond to
statistically strong correlations (Taylor, 1990). In the Appendix, Table 1, we
provide the correlation coefficients for various durations downstream of the
CME-driven shock, which shows insignificant changes in the correlation
coefficients between any two parameters. The results here demonstrate that the
proton SEP intensity can be used to gauge the bulk proton temperature in the
downstream region of CME-driven shocks.

##  4 Discussion

We provide an illustration (Figure 3) to help visualize our explanations of
the observed results that we discuss below. Figure 3 has real data of a CME-
driven shock crossing by PSP for the Feb. 15, 2022, CME event accompanied by
an illustration depicting the structure of a CME and a sample trajectory
across its associated shock.

An interpretation of the results presented in Figure 2 involves the level of
turbulence and the shock strength. Suprathermal particles can be trapped
and/or modulated by the emergence of coherent structures, increasing the
amount of time they spend in the diffusive shock region, which are visualized
in the illustration of Figure 3. If these particles remain sufficiently
trapped, then one might expect a larger intensity of particles observed at
that shock. This property reflects the acceleration efficiency of the shock
and the level of turbulence in these regions, resulting in elevated energy
dissipation (heating) due to the emergence of coherent structures, as can be
measured via PVI (Qudsi et al., 2020, and references therein). In fact, the
PVI shown in Figures 1 and A.4 reveal enhanced values of PVI in the selected
region downstream of the shock, corroborating the correlation of PVI and SEP
intensity (Tessein et al., 2013, 2015; Bandyopadhyay et al., 2020) and/or the
effectiveness of coherent structures in locally modulating/trapping energetic
particles (Mazur et al., 2000; Ruffolo et al., 2003; Kittinaradorn et al.,
2009; Tessein et al., 2016; Khabarova et al., 2021; Pezzi et al., 2021; Pecora
et al., 2021). Additionally, the correlation coefficients between all possible
pairs of PVI, proton temperature, and SEP intensity indicate statistically
strong correlations. As a result, the SEP intensity enables an inference of
the proton temperature from the increased heating due to turbulence
intermittency.

Figure 3: Illustration of CME structure with PSP crossing its driven shock.
The top panel shows the measured proton flux integrated over energy (100 keV
to 500 keV) and proton temperature; the bottom panel shows the magnetic field
magnitude. The CME-driven shock encounter is labeled on PSP’s sketched
trajectory as an open circle, which corresponds to the open circle on the top
axis of the top panel. The trajectory shows PSP traversing the upstream
region, encountering the shock, and then observing the downstream plasma.
Sketch of CME structure is adapted from Cane (2000).

Supplementing the above interpretation, an increased intensity at the shock
can also be explained by the bulk temperature and its impact on the efficiency
of particle injection at shocks. High bulk temperatures increase the
likelihood for particles to have sufficient speed to move upstream and thereby
become injected into diffusive shock acceleration (Gloeckler & Geiss, 1998;
Chotoo et al., 2000; Zurbuchen et al., 2000). This dynamic can be visualized
in Figure 3 where particles move towards the upstream plasma and get picked up
into diffusive shock acceleration. This is due to the increased bulk
temperature raising the thermal core closer to the suprathermal energies,
enabling an increased seed population for acceleration and/or trapping within
the diffusion region. Given the current results, this implies a strong
interplay between solar wind temperature, injection efficiency, and SEP
intensity.

Considerations of other well-correlated properties of CME-driven shocks with
SEP intensity, such as shock normal relative to the magnetic field (Lario et
al., 2016; Giacalone, 2017; Chen et al., 2022), radial distance from the Sun
due to the weakening of such shocks (Salman et al., 2020), and connection to
the shock (Burlaga et al., 1981; Kilpua et al., 2011; Möstl et al., 2009) may
influence the present results. However, many of these other properties are
associated with the shock strength, which is an indicator of the efficiency in
acceleration processes due to the turbulent nature downstream of the shock.
Therefore, we only consider the individual connections of the intensity and
temperature to the level of turbulence in the downstream region.

##  5 Conclusion

In this paper, we investigate the relationship between shock proton
temperature with shock associated particle intensity. The intensity to
temperature correlation had previously only been inferred by existing
correlations of PVI, a measure of turbulence intermittency, to both intensity
and temperature, separately. Our results (Figure 2) suggest the proton SEP
intensity can be used to gauge whether one would expect a larger or smaller
proton temperature immediately downstream of a CME-driven shock. We compute
the PVI to help strengthen the explanations relating the SEP intensity to the
proton temperature, which reveals enhanced PVI in the downstream region.

The results presented here will likely improve the understanding and modeling
of SEP interactions with shocks and the downstream bulk plasma. Additional
spacecraft observations of SEPs and proton temperatures can be used to
significantly increase the statistical significance of these results and
possibly examine their radial profiles. The correction for radial weakening of
shocks and temperature cooling due to expansion will need to be considered
when examining the radial dependence of CME properties in general. Finally, as
the Sun transitions into solar maximum, we anticipate an increase of observed
CME-driven shocks in future PSP orbits that can be used to supplement the
results presented here in addition supporting other systematic studies of CME-
related phenomena.

##  6 Acknowledgements

We thank the IS$\odot$IS/FIELDS/SWEAP teams and everyone that made the PSP
mission possible. The IS$\odot$IS data and visualization tools are available
to the community at https://spacephysics.princeton.edu/missions-
instruments/PSP. PSP was designed, built, and is operated by the Johns Hopkins
Applied Physics Laboratory as part of NASA’s Living with a Star (LWS) program
(contract NNN06AA01C). R. Chhiber is particularly supported by a NASA LWS
grant (80NSSC22K1020). W. H. Matthaeus and F. Pecora are particularly
supported by a PSP/IS$\odot$IS subcontract (SUB0000165) and a NASA HGI grant
(80NSSC21K1765) to the University of Delaware.

##  Appendix A Supplementary Tables and Figures

We consider cases of different durations used in the time-integration and
averaging duration of the intensity and proton temperature, respectively. As
shown in Table 1, changing this duration insignificantly alters the strength
of correlation between any two parameters. We also provide the the correlation
time and the values of proton temperature, SEP intensity, and PVIfrac from
Figure 2 for each event, as shown in Table 2. These values in Table 2
correspond to a time-integration and averaging duration of $3\cdot\tau_{c}$
immediately downstream of the CME-driven shock.

Correlation | $0.5\cdot\tau_{c}$ | $1\cdot\tau_{c}$ | $2\cdot\tau_{c}$ | $3\cdot\tau_{c}$ | $4\cdot\tau_{c}$ | $5\cdot\tau_{c}$
---|---|---|---|---|---|---
$\log\langle T_{p}\rangle$ to $\log{\rm Intensity}$ | 0.94 | 0.92 | 0.94 | 0.95 | 0.96 | 0.95
$\log{\rm Intensity}$ to $\log{\rm PVI}_{frac}$ | 0.81 | 0.81 | 0.81 | 0.81 | 0.82 | 0.83
$\log\langle T_{p}\rangle$ to $\log{\rm PVI}_{frac}$ | 0.87 | 0.84 | 0.82 | 0.81 | 0.79 | 0.81
Table 1: Correlation coefficients for the pair of parameters given in the first column for different time-integration/averaging duration. Overall, the correlations do not significantly alter across different duration based on factors of the correlation time $\tau_{c}$. Event |  $\langle T_{p}\rangle$ [eV] |  SEP Intensity $\left[({\rm cm^{2}\cdot sr\cdot s\cdot keV/nuc})^{-1}\right]$ | PVIfrac |  $\tau_{c}$ [s]
---|---|---|---|---
2020-11-29 | 933 $\pm$ 194 | 237.8 $\pm$ 8.9 | 0.09 | 630
2021-05-28 | 262 $\pm$ 54 | 26.7 $\pm$ 0.9 | 0.02 | 3130
2022-02-15 | 492 $\pm$ 99 | 193.4 $\pm$ 2.1 | 0.11 | 2680
2022-03-21 | 15 $\pm$ 3 | 0.7 $\pm$ 0.1 | 0.02 | 8730
2022-07-09 | 23 $\pm$ 5 | 1.2 $\pm$ 0.2 | 0.01 | 2880
2022-08-28 | 431 $\pm$ 89 | 266.4 $\pm$ 3.0 | 0.20 | 1730
2023-01-03 | 93 $\pm$ 19 | 49.8 $\pm$ 1.4 | 0.01 | 1530
Table 2: Values of $\langle T_{p}\rangle$, SEP intensity, PVIfrac, and
$\tau_{c}$ for each event using a duration of $3\cdot\tau_{c}$ for time-
integration of the energy-integrated particle flux and averaging of $T_{p}$.
These values are reflected in Figure 2.

Figure A.4 provides an overview of the seven events examined in this study.
For each event, there are five panels that represent the EPI-Lo Channel P
omni-directional proton flux spectrogram, magnetic field, ${\rm
PVI}(\tau=10~{}{\rm seconds})$, magnitude of the bulk velocity, proton
density, and proton temperature in descending order. The time range includes
the estimated shock arrival time, marked by the vertical dashed red line. The
time window used for integration of the energy-integrated intensity and
averaging of the proton temperature is between the shock arrival time and the
time marked by the vertical solid red line. These event overviews provide
confirmation of data availability throughout these time periods.

Figure A.4: Overview plots for all seven CMEs. Panels are the EPI-Lo Channel P
omni-directional proton flux spectrogram, magnetic field, PVI at increment
scale $\tau=10~{}{\rm seconds}$, bulk velocity, proton density, and proton
temperature in descending order for each event. The vertical red dashed line
in all panels represents the estimated CME shock arrival time for that
corresponding event.

## References

  * Axford et al. [1977] Axford, W. I., Leer, E., & Skadron, G. 1977, in International Cosmic Ray Conference, Vol. 11, International Cosmic Ray Conference, 132
  * Bale et al. [1999] Bale, S. D., Reiner, M. J., Bougeret, J. L., et al. 1999, Geophys. Res. Lett., 26, 1573
  * Bale et al. [2016] Bale, S. D., Goetz, K., Harvey, P. R., et al. 2016, Space Sci. Rev., 204, 49
  * Bandyopadhyay et al. [2020] Bandyopadhyay, R., Matthaeus, W. H., Parashar, T. N., et al. 2020, ApJS, 246, 61
  * Bell [1978] Bell, A. R. 1978, MNRAS, 182, 147
  * Blandford & Ostriker [1978] Blandford, R. D., & Ostriker, J. P. 1978, ApJ, 221, L29
  * Boggs & Rogers [1990] Boggs, P. T., & Rogers, J. E. 1990, Contemporary mathematics, 112, 183
  * Burlaga et al. [1981] Burlaga, L., Sittler, E., Mariani, F., & Schwenn, R. 1981, J. Geophys. Res., 86, 6673
  * Cane [2000] Cane, H. V. 2000, Space Sci. Rev., 93, 55
  * Cane et al. [1990] Cane, H. V., von Rosenvinge, T. T., & McGuire, R. E. 1990, J. Geophys. Res., 95, 6575
  * Chasapis et al. [2015] Chasapis, A., Retinò, A., Sahraoui, F., et al. 2015, ApJ, 804, L1
  * Chen et al. [2022] Chen, X., Giacalone, J., & Guo, F. 2022, ApJ, 941, 23
  * Chi et al. [2016] Chi, Y., Shen, C., Wang, Y., et al. 2016, Sol. Phys., 291, 2419
  * Chotoo et al. [2000] Chotoo, K., Schwadron, N. A., Mason, G. M., et al. 2000, J. Geophys. Res., 105, 23107
  * Dayeh et al. [2018] Dayeh, M. A., Desai, M. I., Ebert, R. W., et al. 2018, in Journal of Physics Conference Series, Vol. 1100, Journal of Physics Conference Series, 012008
  * Dröge [1994] Dröge, W. 1994, ApJS, 90, 567
  * Fisk [1971] Fisk, L. A. 1971, J. Geophys. Res., 76, 1662
  * Fox et al. [2016] Fox, N. J., Velli, M. C., Bale, S. D., et al. 2016, Space Sci. Rev., 204, 7
  * Germano [1992] Germano, M. 1992, Journal of Fluid Mechanics, 238, 325
  * Giacalone [2017] Giacalone, J. 2017, ApJ, 848, 123
  * Gloeckler & Geiss [1998] Gloeckler, G., & Geiss, J. 1998, Space Sci. Rev., 86, 127
  * Good et al. [2023] Good, S. W., Rantala, O. K., Jylhä, A. S. M., et al. 2023, ApJ, 956, L30
  * Gosling et al. [1979] Gosling, J. T., Asbridge, J. R., Bame, S. J., & Feldman, W. C. 1979, in American Institute of Physics Conference Series, Vol. 56, Particle Acceleration Mechanisms in Astrophysics, ed. J. Arons, C. McKee, & C. Max, 81–99
  * Gosling et al. [1991] Gosling, J. T., McComas, D. J., Phillips, J. L., & Bame, S. J. 1991, J. Geophys. Res., 96, 7831
  * Greco et al. [2008] Greco, A., Chuychai, P., Matthaeus, W. H., Servidio, S., & Dmitruk, P. 2008, Geophys. Res. Lett., 35, L19111
  * Greco et al. [2012] Greco, A., Matthaeus, W. H., D’Amicis, R., Servidio, S., & Dmitruk, P. 2012, ApJ, 749, 105
  * Greco et al. [2018] Greco, A., Matthaeus, W. H., Perri, S., et al. 2018, Space Sci. Rev., 214, 1
  * Guo et al. [2021] Guo, F., Giacalone, J., & Zhao, L. 2021, Frontiers in Astronomy and Space Sciences, 8, 27
  * Hill et al. [2017] Hill, M. E., Mitchell, D. G., Andrews, G. B., et al. 2017, Journal of Geophysical Research (Space Physics), 122, 1513
  * Isaacs et al. [2015] Isaacs, J. J., Tessein, J. A., & Matthaeus, W. H. 2015, Journal of Geophysical Research (Space Physics), 120, 868
  * Jokipii [1982] Jokipii, J. R. 1982, ApJ, 255, 716
  * Jokipii [1987] —. 1987, ApJ, 313, 842
  * Kahler [2001] Kahler, S. W. 2001, J. Geophys. Res., 106, 20947
  * Kallenrode [1996] Kallenrode, M.-B. 1996, J. Geophys. Res., 101, 24393
  * Kasper et al. [2016] Kasper, J. C., Abiad, R., Austin, G., et al. 2016, Space Sci. Rev., 204, 131
  * Khabarova et al. [2015] Khabarova, O., Zank, G. P., Li, G., et al. 2015, ApJ, 808, 181
  * Khabarova et al. [2021] Khabarova, O., Malandraki, O., Malova, H., et al. 2021, Space Sci. Rev., 217, 38
  * Khabarova & Zank [2017] Khabarova, O. V., & Zank, G. P. 2017, ApJ, 843, 4
  * Kilpua et al. [2011] Kilpua, E. K. J., Jian, L. K., Li, Y., Luhmann, J. G., & Russell, C. T. 2011, Journal of Atmospheric and Solar-Terrestrial Physics, 73, 1228
  * Kittinaradorn et al. [2009] Kittinaradorn, R., Ruffolo, D., & Matthaeus, W. H. 2009, ApJ, 702, L138
  * Knock et al. [2001] Knock, S. A., Cairns, I. H., Robinson, P. A., & Kuncic, Z. 2001, J. Geophys. Res., 106, 25041
  * Lario et al. [2017] Lario, D., Decker, R. B., Roelof, E. C., et al. 2017, in Journal of Physics Conference Series, Vol. 900, Journal of Physics Conference Series, 012012
  * Lario et al. [1998] Lario, D., Sanahuja, B., & Heras, A. M. 1998, ApJ, 509, 415
  * Lario et al. [2016] Lario, D., Kwon, R. Y., Vourlidas, A., et al. 2016, ApJ, 819, 72
  * Lee [1983] Lee, M. A. 1983, J. Geophys. Res., 88, 6109
  * Lee [1997] —. 1997, Geophysical Monograph Series, 99, 227
  * Malandraki et al. [2019] Malandraki, O., Khabarova, O., Bruno, R., et al. 2019, ApJ, 881, 116
  * Mann & Klassen [2005] Mann, G., & Klassen, A. 2005, A&A, 441, 319
  * Mason & Gloeckler [2012] Mason, G. M., & Gloeckler, G. 2012, Space Sci. Rev., 172, 241
  * Matthaeus et al. [1999] Matthaeus, W. H., Smith, C. W., & Bieber, J. W. 1999, in American Institute of Physics Conference Series, Vol. 471, Solar Wind Nine, ed. S. R. Habbal, R. Esser, J. V. Hollweg, & P. A. Isenberg, 511–514
  * Mazur et al. [2000] Mazur, J. E., Mason, G. M., Dwyer, J. R., et al. 2000, ApJ, 532, L79
  * McComas et al. [2016] McComas, D. J., Alexander, N., Angold, N., et al. 2016, Space Sci. Rev., 204, 187
  * Möstl et al. [2009] Möstl, C., Farrugia, C. J., Biernat, H. K., et al. 2009, Sol. Phys., 256, 427
  * Neugebauer et al. [1997] Neugebauer, M., Goldstein, R., & Goldstein, B. E. 1997, J. Geophys. Res., 102, 19743
  * Odstrčil [2003] Odstrčil, D. 2003, Advances in Space Research, 32, 497
  * Odstrčil & Pizzo [1999a] Odstrčil, D., & Pizzo, V. J. 1999a, J. Geophys. Res., 104, 483
  * Odstrčil & Pizzo [1999b] —. 1999b, J. Geophys. Res., 104, 493
  * Odstrčil et al. [2004] Odstrčil, D., Riley, P., & Zhao, X. P. 2004, Journal of Geophysical Research (Space Physics), 109, A02116
  * Odstrčil et al. [1996] Odstrčil, D., Smith, Z., & Dryer, M. 1996, Geophys. Res. Lett., 23, 2521
  * Osman et al. [2011] Osman, K. T., Matthaeus, W. H., Greco, A., & Servidio, S. 2011, ApJ, 727, L11
  * Osman et al. [2012] Osman, K. T., Matthaeus, W. H., Hnat, B., & Chapman, S. C. 2012, Phys. Rev. Lett., 108, 261103
  * Pecora et al. [2021] Pecora, F., Servidio, S., Greco, A., et al. 2021, MNRAS, 508, 2114
  * Pezzi et al. [2021] Pezzi, O., Pecora, F., Le Roux, J., et al. 2021, Space Sci. Rev., 217, 39
  * Qudsi et al. [2020] Qudsi, R. A., Maruca, B. A., Matthaeus, W. H., et al. 2020, ApJS, 246, 46
  * Reames [1999] Reames, D. V. 1999, Space Sci. Rev., 90, 413
  * Reames [2000] Reames, D. V. 2000, in American Institute of Physics Conference Series, Vol. 516, 26th International Cosmic Ray Conference, ICRC XXVI, ed. B. L. Dingus, D. B. Kieda, & M. H. Salamon, 289–300
  * Richardson & Cane [2010] Richardson, I. G., & Cane, H. V. 2010, Sol. Phys., 264, 189
  * Ruffolo et al. [2003] Ruffolo, D., Matthaeus, W. H., & Chuychai, P. 2003, ApJ, 597, L169
  * Salman et al. [2020] Salman, T. M., Winslow, R. M., & Lugaz, N. 2020, Journal of Geophysical Research (Space Physics), 125, e27084
  * Servidio et al. [2011] Servidio, S., Greco, A., Matthaeus, W. H., Osman, K. T., & Dmitruk, P. 2011, Journal of Geophysical Research (Space Physics), 116, A09102
  * Shen et al. [2023] Shen, M. M., Szalay, J. R., Pokorny, P., et al. 2023, in prep.
  * Szalay et al. [2020] Szalay, J. R., Pokorný, P., Bale, S. D., et al. 2020, ApJS, 246, 27
  * Taylor [1990] Taylor, R. 1990, Journal of diagnostic medical sonography, 6, 35
  * Tessein et al. [2013] Tessein, J. A., Matthaeus, W. H., Wan, M., et al. 2013, ApJ, 776, L8
  * Tessein et al. [2016] Tessein, J. A., Ruffolo, D., Matthaeus, W. H., & Wan, M. 2016, Geophys. Res. Lett., 43, 3620
  * Tessein et al. [2015] Tessein, J. A., Ruffolo, D., Matthaeus, W. H., et al. 2015, ApJ, 812, 68
  * Trotta et al. [2021] Trotta, D., Valentini, F., Burgess, D., & Servidio, S. 2021, Proceedings of the National Academy of Science, 118, e2026764118
  * Trotta et al. [2023] Trotta, D., Pezzi, O., Burgess, D., et al. 2023, MNRAS, 525, 1856
  * Wu & Lepping [2011] Wu, C.-C., & Lepping, R. P. 2011, Sol. Phys., 269, 141
  * Yordanova et al. [2016] Yordanova, E., Vörös, Z., Varsani, A., et al. 2016, Geophys. Res. Lett., 43, 5969
  * Zank et al. [2015] Zank, G. P., Hunana, P., Mostafavi, P., et al. 2015, ApJ, 814, 137
  * Zurbuchen et al. [2000] Zurbuchen, T. H., Fisk, L. A., Schwadron, N. A., & Gloeckler, G. 2000, in American Institute of Physics Conference Series, Vol. 528, Acceleration and Transport of Energetic Particles Observed in the Heliosphere, ed. R. A. Mewaldt, J. R. Jokipii, M. A. Lee, E. Möbius, & T. H. Zurbuchen, 215–220
  * Zurbuchen & Richardson [2006] Zurbuchen, T. H., & Richardson, I. G. 2006, Space Sci. Rev., 123, 31