# Large-scale Kinetic Simulations of Colliding Plasmas within a Hohlraum of
Indirect Drive Inertial Confinement Fusions

Tianyi Liang  Institute for Fusion Theory and Simulation, School of Physics,
Zhejiang University, Hangzhou 310058, China  Dong Wu<EMAIL_ADDRESS>Key
Laboratory for Laser Plasmas and School of Physics and Astronomy, and
Collaborative Innovation Center of IFSA, Shanghai Jiao Tong University,
Shanghai 200240, China  Xiaochuan Ning  Institute for Fusion Theory and
Simulation, School of Physics, Zhejiang University, Hangzhou 310058, China
Lianqiang Shan  Science and Technology on Plasma Physics Laboratory, Research
Center of Laser Fusion, CAEP, Mianyang 621900, China  Zongqiang Yuan  Science
and Technology on Plasma Physics Laboratory, Research Center of Laser Fusion,
CAEP, Mianyang 621900, China  Hongbo Cai  Institute of Applied Physics and
Computational Mathematics, CAEP, Beijing 100094, China  Zhengmao Sheng
<EMAIL_ADDRESS>Institute for Fusion Theory and Simulation, School of
Physics, Zhejiang University, Hangzhou 310058, China  Xiantu He  Institute for
Fusion Theory and Simulation, School of Physics, Zhejiang University, Hangzhou
310058, China

###### Abstract

The National Ignition Facility has recently achieved successful burning plasma
and ignition using the inertial confinement fusion (ICF) approach. However,
there are still many fundamental physics phenomena that are not well
understood, including the kinetic processes in the hohlraum. Shan et al.
[Phys. Rev. Lett, 120, 195001, 2018] utilized the energy spectra of neutrons
to investigate the kinetic colliding plasma in a hohlraum of indirect drive
ICF. However, due to the typical large spatial-temporal scales, this
experiment could not be well simulated by using available codes at that time.
Utilizing our advanced high-order implicit PIC code, LAPINS, we were able to
successfully reproduce the experiment on a large scale of both spatial and
temporal dimensions, in which the original computational scale was increased
by approximately 7 to 8 orders of magnitude. When gold plasmas expand into
deuterium plasmas, a kinetic shock is generated and propagates within
deuterium plasmas. Simulations allow us to observe the entire progression of a
strong shock wave, including its initial formation and steady propagation.
Although both electrons and gold ions are collisional (on a small scale
compared to the shock wave), deuterium ions seem to be collisionless. This is
because a quasi-monoenergetic spectrum of deuterium ions can be generated by
reflecting ions from the shock front, which then leads to the production of
neutrons with unusual broadening due to beam-target nuclear reactions. This
work displays an unprecedented kinetic analysis of an existing experiment,
shedding light on the mechanisms behind shock wave formation. It also serves
as a reference for benchmark simulations of upcoming new simulation codes and
may be relevant for future research on mixtures and entropy increments at
plasma interfaces.

##  I INTRODUCTION

The recent experiments carried out at the National Ignition Facility (NIF) [1,
2, 3, 4, 5, 6] have successfully validated the feasibility of controllable
inertial confinement fusion (ICF). This achievement is considered a noteworthy
milestone in the endeavor to attain inexhaustible and environmentally friendly
sources of energy. This remarkable achievement can be attributed to the
collaborative efforts involving various technological advancements and
physical insights. Some notable contributions include the utilization of high-
density-carbon (HDC) capsules in low-gas-fill hohlraums [7] and the
implementation of the “BigFoot” (BF) [8] scheme. Moreover, a comprehensive set
of scaling laws [9] and evaluation metrics [10] based on the analysis of
experimental data have been consolidated, serving as crucial tools in
comprehending the underlying physical mechanisms of ICF and providing guidance
for the development of pertinent experimental designs. However, from a
theoretical standpoint, there remain numerous foundational physics concepts
that are not yet comprehensively understood, particularly those that are
closely tied to kinetic and non-equilibrium processes.

Figure 1: Subgraph (a) depicts the schematic diagram of the colliding plasmas
within a hohlraum. It illustrates various phenomena such as collisional
interaction, shock formation, collisionless ion reflection, and nuclear
reactions. Subgraph (b) represents the simulation setup.

In vacuum or near-vacuum hohlraums [11], as depicted in Fig. 1 (a), the
expansion of high-Z plasma from the inner wall and its collision with the
blow-off from the capsule or the low-density-fill gas [12, 13], can lead to
the dominance of kinetic effects. This is due to the ion-ion mean free paths
being larger than the size of the interaction regions. The existence of these
kinetic effects has been substantiated through experiments and simulations
[14, 15, 16], and their influence on the implosion process of ICF has been
partially elucidated [17]. To attain high-gain laser fusion, it is imperative
to thoroughly comprehend and regulate any notable physics phenomena that may
arise during ICF processes. Therefore, the study of kinetic effects and their
impact on energy depositions, implosion symmetries, and other inherent plasma
properties holds significant importance.

Recently, in an experiment conducted by Shan et al. [14], a kinetic colliding
plasma was observed within a hohlraum of indirect drive ICF. This observation
was made by measuring the energy spectra of neutrons. The width of the spectra
in this experiment is unusually large. To explain this, they conducted
collisionless particle-in-cell (PIC) simulations to investigate the
interactions between the gold plasma and the low-density-fill gas. It was
found that the kinetic shock wave reflects upstream ions and results in beam-
target fusion. This phenomenon could, in turn, explain the broadening of the
neutron energy spectrum. However, the PIC simulation is somewhat ideal due to
the absence of collisional effects and nuclear reactions. Moreover, the scale
of the PIC simulation was only several micrometers, much smaller than the
scale of real physics, which is hundreds of micrometers in the experiment.

As the conventional PIC method is constrained by computational resources,
making it unsuitable for large-scale simulations, scholars have put forward
hybrid approaches that combine the advantages of both fluid and kinetic PIC
methods [18, 19, 20]. While in the context of the fluid component, the authors
of this study still rely on empirical coefficients, such as the flux limiter
and electron-ion coupling coefficients, for their treatment. Additionally, the
authors make use of various approximations in different forms. The utilization
of the hybrid approach assumes that the electrons must adhere to the condition
of fluid approximation, wherein the mean free paths of the electrons should be
smaller than the spatial resolution of the simulation. But in hohlraums or
other laboratory astrophysics experiments, where the electron density ranges
$10^{18}\>\rm{cm^{-3}}$ to $10^{20}\>\rm{cm^{-3}}$, the behavior of electrons
cannot be accurately described using fluid dynamics. As a result, hybrid
approaches may not be applicable in these scenarios. To handle a wide range of
time scales, space scales, and densities, we developed a high-order implicit
multidimensional PIC method. With an appropriate arrangement of space and
time, the proposed method can significantly minimize numerical errors.
Additionally, the utilization of a higher-order interpolation method can
effectively reduce numerical noises. This reduction in numerical noise is
particularly useful for simulating large-scale kinetic processes using the PIC
method.

In this paper, the experiment conducted by Shan et al. [14] was simulated and
analyzed using the high-order implicit PIC code LAPINS. The simulation
performed was analyzed in a large-scale set-up that closely resembles the real
experiment. In our simulations, not only the density distributions and ion
reflections associated with the shock wave were figured out, but also the
neutron spectra that arise from the nuclear reactions between the reflected
ions and the capsule were self-consistently generated. Considering the
collision effects, the expansion of gold plasma manifests as a central
rarefaction wave, instead of an isothermal one simulated in the reference
[14]. When the gold plasma expands into the deuterium plasma, it initiates the
generation of a kinetic shock wave, which subsequently propagates within the
deuterium plasma. The shock wave under consideration is situated in the
intermediate region between collisionless and collisional regimes. The primary
mechanisms accountable for dissipation in this particular context involve both
the reflections of upstreaming deuterium ions and collisional processes.
Deuterium ions exhibiting quasi-monoenergetic spectra can be produced through
the collisionless electrostatic shock wave mechanism, resulting in the
generation of neutrons with significantly broadened energy distribution due to
beam-target nuclear reactions.

The paper is organized as follows. A brief introduction to the high-order
implicit PIC method and the pairwise nuclear reaction algorithm of the LAPINS
code is described in Section II. In Section III, the simulation results are
presented and analyzed in detail. Finally, the discussion and conclusion are
displayed in Section IV.

##  II  Simulation methods

The LAPINS code has undergone significant advancements over the years,
including the incorporation of collision effects [21], the consideration of
quantum degeneracy [22], and the adoption of a hybrid-kinetic approach [23].
In this section, a brief introduction to the module of nuclear reactions [24],
and the high-order implicit PIC method [25] utilized in our simulation.

###  II.1  IMPLICIT PIC MODEL

In the LAPINS code, a high-order implicit multidimensional PIC method has been
devised to effectively address the complexities of astrophysics and dense
plasmas. The spatial-temporal arrangement is established by employing Yee’s
algorithm in conjunction with a leapfrog algorithm to simulate the propagation
of electromagnetic fields and the advancement of particles. Specifically, the
charge density is positioned at the Yee cell centers, while the electric
fields and current density are staggered upwards to the cell faces.
Additionally, the magnetic fields are located at the cell edges. This
arrangement considers the discretization of the Faraday and Ampere equations.
It can be demonstrated that the equation $\nabla\times B=0$ is consistently
fulfilled when it is initially valid.

Our field solver algorithm efficiently tackles the numerical instabilities
commonly encountered in explicit PIC methods when using relaxed time steps and
grid resolution. The algorithm used in this study effectively addresses the
problem of numerical cooling, a common issue in standard implicit PIC methods,
by employing a pseudo-electric-field approach. The violation of Gauss’s law in
the cell at time $t^{n+1}$ is denoted as $F^{n+1}=\nabla\cdot
E^{n+1}-\rho^{n+1}$. It is used to calculate the additional pseudo electric
field $E^{n+1}_{\text{psd}}=d\Delta t\nabla F^{n+1}$, with the introduction of
a customized dimensionless number $d$. Consequently, in each time increment,
the overall electric field $E^{n+1}$ exerted on particles needs to be updated
by adding the pseudo electric field, which can be expressed as
$E^{n+1}=E^{n+1}+E^{n+1}_{\text{psd}}$.

The particle pusher algorithm is a combination of the standard Boris particle
pusher and the Newton-Krylov iteration method. This algorithm has the
potential to greatly enhance precision accuracy, surpassing the standard Boris
particle pusher by several orders of magnitude. Additionally, it offers a
substantial reduction in iteration time compared to the pure Newton-Krylov
method. For further information, readers are encouraged to consult our recent
publication [25].

###  II.2  MODEL OF NUCLEAR REACTION

We have successfully developed a model for pairwise nuclear reactions
involving weighted particles at relativistic energies [24]. It is worth
mentioning that the particle-pairing routine used for fusion reactions is
identical to the one employed for binary collisions, known as the Takizuka-Abe
algorithm. In every spatial cell, pairs of particles engaging in nuclear
reactions are randomly selected. Energy and momentum exchanges are calculated
for every pair of particles and performed in the center-of-momentum (CM)
frame. Our model is additionally applicable in the context of the relativistic
regime. In the CM frame, the probability of the reaction, denoted as
$P_{\text{ab}}$, is determined by the calculation of

$P_{\text{ab}}=n_{\text{min}}\sigma\rm{v}_{rel,CM}\gamma_{\text{CM}}\Delta t,$
(1)

where $n_{\text{min}}$ is the minimum density between particle species $a$ and
$b$, $\sigma_{\text{ab}}$ denotes the cross-section of nuclear reaction, and
$\Delta t$ is the simulation time step, which is multiplied by a factor of
$\gamma_{\text{CM}}$ when taken into account in the CM frame. The nuclear
reaction cross sections used in the program are obtained from the
International Atomic Energy Agency (IAEA) and are stored within the code as
Legendre polynomial coefficients. To calculate these values, the relative
velocity between the two particles in the CM frame is required and provided by

$\rm{v}_{\text{rel,CM}}=\left|\frac{\mathbf{v}_{a,CM}-\mathbf{v}_{b,CM}}{1-\mathbf{v}_{a,CM}\cdot\mathbf{v}_{b,CM}}\right|.$
(2)

Therefore, the yield of nuclear reaction for each pair of macro-particles is
determined as

$Y_{\text{ab}}=w_{\text{min}}P_{\text{ab}},$ (3)

where $w_{\text{min}}$ represents the minimum weight between macro-particles
$a$ and $b$. To enhance the calculation accuracy, we have introduced a
variable parameter $F_{\text{multi}}$ to increase the probability of the
reaction while decreasing the weight of the products, denoted as
$w_{\text{min}}\to w_{\text{min}}/F_{\text{multi}}$ and $P_{\text{ab}}\to
P_{\text{ab}}F_{\text{multi}}$, while keeping $Y_{\text{ab}}$ unaltered. For
details of the nuclear fusion scheme, one can refer to our recent paper and
other relevant publications [24].

##  III  SIMULATION RESULTS

Figure 2: The subgraphs (a)-(c) depicted in the left part illustrate the
temporal evolution of the density profiles of electrons, deuterium ions
($\rm{D^{+}}$), and gold ions ($\rm{Au^{50+}}$), respectively. The black
dashed line depicted the trajectory of the shock front, characterized by a
velocity of $V_{sh}=662\>\rm{km/s}$, while the dotted line represents the
trajectory with double the velocity, $2V_{sh}$. The subgraphs (d)-(f) on the
right side depict the phase space distributions of gold ions and deuterium
ions at different times. The white-black colormap and colorful colormap are
used to represent the phase density of gold ions and deuterium ions,
respectively.

To conduct a more in-depth analysis of the experimental data [14], the LAPINS
code is employed, utilizing parameters that closely resemble those of the
experimental setup. This simulation involves the collision interaction between
the ablated gold plasma and deuterium gas, as depicted in Fig. 1 (b). The
simulation box has a length of $L=3\>\rm{mm}$, which is resolved by $3000$
cells. Each cell contains $1000$ macro-particles for electrons and $400$
macro-particles for various ions. Simulation times in this study are on the
order of nanoseconds, suggesting that the scale of the simulation is
approximately 7 to 8 orders of magnitude larger than that of traditional PIC
codes, which typically operate on a temporal scale of picoseconds and a
spatial scale of micrometers.

The gold plasma is located on the left side of the simulation box. The
ionization degree of gold ions is fixed at $Z=50$, while the initial ion
temperature is $T_{\text{Au}}=100\>\rm{eV}$. The initial electron temperature
of the gold plasma is $T_{e1}=3000\>\rm{eV}$, while the electron density is
$n_{e1}=Zn_{\text{Au}}=1.0\times 10^{21}\>\rm{cm^{-3}}$. Full-ionized
deuterium plasma is situated on the right side of the simulation box, with a
density of $n_{e0}=n_{\text{D}}=2.0\times 10^{19}\>\rm{cm^{-3}}$ and
temperature of $T_{\text{D}}=T_{\text{e0}}=100\>\rm{eV}$. Electrostatic shock
structure is generated from the expansion of gold plasma into deuterium
plasma. The gold plasma is commonly recognized as the downstream region
located behind the shock, whereas the deuterium plasma is considered the
upstream region situated ahead of the shock. An absorbing layer is designated
to function as the capsule characterized by a high density of $2.0\times
10^{21}\>\rm{cm^{-3}}$ and a low temperature of $1\>\rm{eV}$. The speed of
sound in deuterium plasma is
$c_{s}=\sqrt{{kT_{e0}}/{m_{\text{D}}}}=126.4\>\rm{km/s}$. It is important to
note that there is no initial drift velocity for any species in the plasmas.

In the initial stages of the expansion, the hotter electrons move faster than
the gold ions. This differential motion generates an electrostatic field
characterized by charge separation, commonly referred to as the sheath
electrostatic field. This electric field plays an important role in both the
expansion of the gold plasma and the reflection of the deuterium plasma.
Plasma expansion is a phenomenon wherein the internal energy of plasma is
converted into kinetic energy. When considering the collisional effect, the
expansion behaves as a centered rarefaction wave [26] instead of an isothermal
one, as the deuterium plasma acts as a piston. During the process of
expansion, the velocity of the gold plasma exhibits a linear increase until it
reaches its maximum value. It is illustrated in Fig. 2 (d)-(f), which depict
the phase space distributions of gold ions and deuterium ions at various time
intervals. The collision frequency between deuterium and gold ions,
$\nu_{\text{D,Au}}=4/3\pi^{1/2}({Z_{\text{Au}}Z_{\text{D}}}/{4\pi\varepsilon_{0}})^{2}n_{\text{D}}{m_{\text{D}}^{-1/2}T_{\text{D}}^{-3/2}}\ln\Lambda$,
is estimated to be in the range of $10^{12}\sim 10^{13}/\rm{s}$, and the mean
free path is approximately $l_{\text{D,Au}}\approx 0.15\>\rm{\mu m}$. In the
context of the penetration region, the presence of intense collisions between
deuterium ions and gold ions hinders the entry of deuterium ions into the
downstream region.

Therefore, during the initial formation of the shock wave, the shock strength
denoted as $\delta=n_{\text{Dmax}}/n_{\text{D}}$ and referred to as the
density compression ratio, rapidly reaches a value of 6. Subsequently, after
approximately 1 ns of simulation time, it gradually decreases to 3.5 and
remains stable. The decrease can be attributed to the expansion of the dense
deuterium plasma in the shock region to the upstream deuterium plasma. In the
stable stages, the shock strength $\delta$ remains below the upper limit
determined by the Rankine-Hugoniot equations, indicating that $\delta$ is less
than 4.

Figure 3: The temporal evolution of the energy spectra of deuterium ions. The
total spatial energy spectra are presented in subgraph (a), while the time-
integrated energy spectra obtained at the location of $z=1200\>\rm{\mu m}$ are
displayed in subgraph (b).

One of the notable kinetic properties exhibited by electrostatic shock is the
reflection of upstream ions. The electric potential difference, $\Delta\phi$,
can reflect any ions with kinetic energies lower than it, as expressed by
$m_{i}v_{i}^{2}/2<e\Delta\phi$ [27, 28]. In our simulations, it is evident
that the reflection of deuterium ions can be categorized into two distinct
phases. Firstly, there is a pronounced acceleration of the sheath electric
field at the initial stage, which contributes to the high-energy portion of
the reflected ions. Secondly, there is a consistent reflection of the
electrostatic shock wave, with a speed that can reach up to twice the speed of
the shock wave ($2V_{\text{sh}}$) in the laboratory frame. The velocity of the
shock, denoted as $V_{\text{sh}}$ was determined based on the temporal
evolution in the density distribution of electrons and ions depicted in Fig. 2
(a)-(c). The measured velocity of the shock wave is
$V_{\text{sh}}=662\>\rm{km/s}$ and the Mach number is
$M=V_{\text{sh}}/c_{s}=5.25$.

In the phase space depicted in Fig. 2 (d)-(f), the reflected deuterium ions
are accelerated to velocities exceeding $2000\>\rm{km/s}$, with a steady
velocity of approximately $1300\>\rm{km/s}$, which is nearly twice
$V_{\text{sh}}$. In the theoretical framework of the KdV-B equation, the
phenomenon of steady reflection can be understood as a form of dissipation.
The collisional effect does not provide sufficient dissipation due to the low
density of the deuterium plasma. As a result, reflection compensates for
dissipation, which plays an important role in forming the shock as well as
dispersion and nonlinearity. When dissipation, dispersion, and nonlinearity
are balanced, the shock wave can propagate steadily, otherwise, it will
collapse or disappear [28].

Diagnosing the evolution of the total energy spectra of deuterium ions over
time, as depicted in Fig. 3 (a), it is evident that the energetic ions
generated through reflection exhibit substantial deviations from the initial
narrow Maxwell distribution. A time-integral diagnostic plane is established
at the location of $z=1200\>\rm{\mu m}$ to record the energy spectra of
reflected ions, which can arrive at the absorbing layer, as shown in Fig. 3
(b). Deuterium ions exhibiting quasi-monoenergetic spectra are observed,
reaching a maximum energy of approximately $50\>\rm{keV}$, with a full width
at half maximum (FWHM) $\Delta E_{\text{D}}=30\>\rm{keV}$. The spectra undergo
broadening when the shock wave passes through the diagnostic plane. When these
reflected ions are deposited inside the capsule, significant low-mode
asymmetry may occur in the implosion process [29], as well as the beam-target
nuclear reactions.

To simulate the deposition of reflected ions inside the capsule, an absorbing
layer located at $2000\>\rm{\mu m}$ is set to capture the energetic ions
accelerated by the shock wave. As the probability of the nuclear reaction
$\rm{D(D,n)^{3}He}$ is very small for reflected ions with a speed of several
hundred kilometers per second and a temperature of $1\>\rm{eV}$ in the plasma
layer, we have set the parameter $F_{\text{multi}}$ as $10000$ using more than
1 million deuterium macro-particles. The statistical and fitting data of the
neutron spectra are shown in Fig. 4 (b). The peak energy is $2.45\>\rm{MeV}$
with a FWHM of $\Delta E_{n}=0.3\>\rm{MeV}$. The relationship between FWHM and
$T_{\text{D}}$ in thermonuclear conditions can be described as $\Delta
E_{n}=2\sqrt{\ln{2}\langle E_{n}\rangle T_{\text{D}}}=82.5T_{\text{D}}^{1/2}$
[30], where $T_{\text{D}}$ is in keV and $\langle
E_{n}\rangle=2.45\>\rm{MeV}$. If we assume that this unusual broadening is due
to thermonuclear fusion, the estimated temperature ($T_{\text{D}}$) is
approximately $13\>\rm{keV}$, which is significantly higher than the
temperature achievable for the high-density absorbing layer ($1\>\rm{eV}$).

By utilizing double-differential cross-sections,
$\rm{d}^{2}\sigma/\rm{d}E\rm{d}\Omega$, the energies of emitted neutrons can
be observed as a function of scattering angle $\theta$ in the CM frame,
denoted as,

$\displaystyle
E_{n}=\frac{E_{\text{D}}m_{n}}{4m_{\text{D}}}\left[\frac{2m_{\text{D}}-m_{n}}{m_{n}}\left(2\frac{Q}{E_{\text{D}}}+1\right)\right.$
$\displaystyle\left.+2\sqrt{\frac{2m_{\text{D}}-m_{n}}{m_{n}}\left(2\frac{Q}{E_{\text{D}}}+1\right)}\cos{\theta}+1\right],$
(4)

where $E_{\text{D}}$ is the incident energy of deuterium ions, and
$Q=3.27\>\rm MeV$ is the released energy of $\rm{D(D,n)^{3}He}$ nuclear
reaction. The numerical results obtained from Eq. 4 for various incident
energies of deuterium ions are shown in Fig. 4 (c). When using the value of
our simulations, e.g. $E_{\text{D}}=50\>\rm{keV}$, we find that $\Delta
E_{n,\theta}=E_{n}(\theta=0)-E_{n}(\theta=\pi)=497\>\rm{keV}$, which is well
consistent with the spectra, as shown in Fig. 4 (b). The spectra exhibit
natural broadening due to beam-target reactions when the neutrons are counted
by integrating the full solid angle. This approach was used in our simulation.
As illustrated in Fig. 4 (a), during the experiment, the neutron spectrometer
was arranged in a narrow solid angle configuration to capture neutrons with
varying scattering angles $\theta$ emitted from the entire hohlraum. This
experimental setup is mathematically equivalent to our one-dimensional
simulations, where all neutrons are accounted for by integrating over the full
solid angle. When comparing the FWHM of the neutron spectra in our simulations
with that of the experiment data (represented by light blue circles in Fig. 4
(b)), a high level of consistency is observed.

Figure 4: Subgraph (a) is a top-view schematic diagram of the setup of a
neutron spectrometer in a hohlraum experiment. Subgraph (b) represents the
comparison of the neutron spectra obtained in the simulation with that of the
experimental data. The solid lines of different colors are the simulated data
over different electron temperatures $T_{e1}$, fitted with Gaussian profiles.
The experimental data presented by Shan et al. [14] is depicted as light blue
circles. Subgraph (c) represents the numerical results obtained from Eq. 4 for
various incident energies of deuterium ions, denoted as $E_{\text{D}}$.

##  IV DISCUSSION AND CONCLUSION

In a hohlraum experiment, the determination of the state of the gold plasma
poses a significant challenge, as it directly influences the properties of the
shock wave and the behavior of the reflected ions. To elucidate the impacts,
we manipulated the electron temperature $T_{e1}$ of the gold plasma. The phase
space distributions are shown in Fig. 5 (a)-(c). As $T_{e1}$ decreases, the
ions accelerated by the sheath electrostatic field have lower energy and
constitute a smaller portion of the reflected ions. The velocity of the shock
wave, $V_{\text{sh}}$ also decreases. According to the collisionless laminar
shock wave kinetic theory [31], a shock wave can arise from the collision
between two plasma slabs with different temperatures and densities. Ion
reflection occurs when the electrostatic potential surpasses the kinetic
energy of upstream ions. There exists a critical Mach numbers $M_{\text{cr}}$
for reflecting ions, that can be calculated numerically by [32],

$\displaystyle
M^{2}_{\text{cr}}=\frac{1}{1+\Gamma}\left[\frac{\sqrt{2}M_{\text{cr}}}{\sqrt{\pi}}+e^{\frac{M^{2}_{\text{cr}}}{2}}\text{Erfc}\frac{M_{\text{cr}}}{\sqrt{2}}-1\right.+\Gamma\Theta\times$
$\displaystyle\left.\left(\frac{\sqrt{2}M_{\text{cr}}}{\sqrt{\pi\Theta}}+e^{\frac{M^{2}_{\text{cr}}}{2\Theta}}\text{Erfc}\frac{M_{\text{cr}}}{\sqrt{2\Theta}}+\frac{4M^{3}_{\text{cr}}}{3\sqrt{2\pi\Theta^{3}}}-1\right)\right],$
(5)

where $\Gamma=n_{e1}/n_{e0}$, and $\Theta=T_{e1}/T_{e0}$, are the ratios of
upstream and downstream electron densities and temperatures. In the
simulation, the $\Gamma=51$ is fixed. When $M^{2}\gg 1$ and
$M^{2}\gg\sqrt{\Theta}$, there is an upper limit to the Mach number [31],

$M_{\text{max}}\simeq\frac{3(\Gamma+1)}{\Gamma}\sqrt{\frac{\pi\Theta}{8}}.$
(6)

The Mach numbers estimated in our simulations via different $\Theta$ are
presented in Fig. 5 (d), which is well located between two lines. The Mach
number can be mathematically represented as a function of $\Theta$ and
$\Gamma$. Consequently, the characteristics of the shock wave are primarily
determined by the distinct properties of electrons in the upstream and
downstream regions.

Figure 5: Subgraphs (a)-(c) depict the phase space distributions for gold ions
and deuterium ions concerning electron temperature $T_{e1}$ of the gold
plasma. The white-black colormap and colorful colormap represent the phase
density of gold ions and deuterium ions, respectively. Subgraph (d) depicts a
comparison of Mach numbers obtained from simulations (represented by blue
crosses) and theoretical calculation results from Eq. 5 and Eq. 6 (represented
by blue and red solid lines).

The corresponding neutron spectra are shown in Fig. 4 (b). As the energy of
the reflected ions decreases, the neutron energy spectra become narrower,
aligning with the theoretical prediction of the beam-target reaction. The FWHM
of the neutron spectra obtained from the experiments can serve as an indicator
for estimating the electron temperature of the gold plasmas.

To summarize, Shan et al. [14] conducted a study on a kinetic colliding plasma
within a hohlraum of indirect drive ICF by measuring the energy spectra of
neutrons. However, due to the typical large spatial-temporal scales involved,
this experiment could not be accurately simulated using available codes at
that time. The experiment was successfully replicated at large spatial and
temporal scales using our high-order implicit PIC code, LAPINS. The
computational scale of our simulations is approximately 7 to 8 orders of
magnitude greater than that of traditional PIC codes. When gold plasmas expand
into deuterium plasmas, a kinetic shock is generated and propagates within
deuterium plasmas. Through the utilization of simulations, we can observe the
complete progression of a strong shock wave, encompassing its initial
formation and subsequent steady propagation. Although both electrons and gold
ions are collisional, deuterium ions appear to be collisionless. The quasi-
monoenergetic spectra of deuterium ions can be generated by reflecting ions
from the shock front. This process leads to the production of neutrons with
unusual broadening due to beam-target nuclear reactions. This study provides
an unprecedented kinetic analysis of an existing experiment, which contributes
to our understanding of the mechanisms underlying the formation of shock
waves. It can be relevant for future research on mixtures and entropy
increments at plasma interfaces.

###### Acknowledgements.

This work is supported by National Natural Science Foundation of China (Grants
No. 12075204 and No. 11875235), Science and Technology on Plasma Physics
Laboratory Foundation of Chinese Academy of Engineering Physics, the Strategic
Priority Research Program of Chinese Academy of Sciences (XDA250050500), and
Shanghai Municipal Science and Technology Key Project (No. 22JC1401500). Dong
Wu thanks the sponsorship from Yangyang Development Fund.

## References

  * Abu-Shawareb  _et al._ [2022] H. Abu-Shawareb, R. Acree, P. Adams, _et al._ (Indirect Drive ICF Collaboration), Lawson criterion for ignition exceeded in an inertial fusion experiment, Phys. Rev. Lett. 129, 075001 (2022).
  * Zylstra  _et al._ [2022a] A. B. Zylstra, O. A. Hurricane, D. A. Callahan, _et al._ , Burning plasma achieved in inertial fusion, Nature 601, 542 (2022a).
  * Hartouni  _et al._ [2023] E. P. Hartouni, A. S. Moore, A. J. Crilly, _et al._ , Evidence for suprathermal ion distribution in burning plasmas, Nature Physics 19, 72 (2023).
  * Zylstra  _et al._ [2022b] A. B. Zylstra, A. L. Kritcher, O. A. Hurricane, _et al._ , Experimental achievement and signatures of ignition at the national ignition facility, Physical Review E 106, 025202 (2022b).
  * Atzeni  _et al._ [2022] S. Atzeni, D. Batani, C. N. Danson, _et al._ , Breakthrough at the nif paves the way to inertial fusion energy, Europhysics News 53, 18 (2022).
  * Zylstra  _et al._ [2021a] A. B. Zylstra, A. L. Kritcher, O. A. Hurricane, _et al._ , Record energetics for an inertial fusion implosion at nif, Physical Review Letters 126, 025001 (2021a).
  * Hopkins  _et al._ [2018] L. B. Hopkins, L. Divol, C. Weber, _et al._ , Increasing stagnation pressure and thermonuclear performance of inertial confinement fusion capsules by the introduction of a high-z dopant, Physics of Plasmas 25, 080706 (2018).
  * Baker  _et al._ [2020] K. L. Baker, C. A. Thomas, D. T. Casey, _et al._ , Hotspot parameter scaling with velocity and yield for high-adiabat layered implosions at the national ignition facility, Physical Review E 102, 023210 (2020).
  * Hurricane  _et al._ [2019] O. A. Hurricane, D. A. Callahan, P. T. Springer, _et al._ , Beyond alpha-heating: driving inertially confined fusion implosions toward a burning-plasma state on the national ignition facility, Plasma Physics and Controlled Fusion 61, 014033 (2019).
  * Zylstra  _et al._ [2021b] A. B. Zylstra, O. A. Hurricane, D. A. Callahan, _et al._ , Metrics for implosion performance with enhanced energy coupling on nif, Nuclear Fusion 61, 116066 (2021b).
  * Pape  _et al._ [2016] S. L. Pape, L. F. B. Hopkins, L. Divol, _et al._ , The near vacuum hohlraum campaign at the nif: A new approach, Physics of Plasmas 23, 056311 (2016).
  * Hopkins  _et al._ [2015a] L. B. Hopkins, N. Meezan, S. L. Pape, _et al._ , First high-convergence cryogenic implosion in a near-vacuum hohlraum, Physical Review Letters 114, 175001 (2015a).
  * Hopkins  _et al._ [2015b] L. F. B. Hopkins, S. L. Pape, L. Divol, _et al._ , Near-vacuum hohlraums for driving fusion implosions with high density carbon ablators, Physics of Plasmas 22, 056318 (2015b).
  * Shan  _et al._ [2018] L. Q. Shan, H. B. Cai, W. S. Zhang, _et al._ , Experimental evidence of kinetic effects in indirect-drive inertial confinement fusion hohlraums, Phys. Rev. Lett. 120, 195001 (2018).
  * Pape  _et al._ [2014] S. L. Pape, L. Divol, L. B. Hopkins, _et al._ , Observation of a reflected shock in an indirectly driven spherical implosion at the national ignition facility, Physical Review Letters 112, 225002 (2014).
  * Li  _et al._ [2010] C. K. Li, F. H. Séguin, J. A. Frenje, _et al._ , Charged-particle probing of x-ray-driven inertial-fusion implosions, Science 327, 1231 (2010).
  * Rinderknecht  _et al._ [2018] H. G. Rinderknecht, P. A. Amendt, S. C. Wilks, and G. Collins, Kinetic physics in icf: present understanding and future directions, Plasma Physics and Controlled Fusion 60, 064001 (2018).
  * Welch  _et al._ [2001] D. R. Welch, D. V. Rose, B. V. Oliver, and R. E. Clark, Simulation techniques for heavy ion fusion chamber transport, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 464, 134 (2001).
  * Cai  _et al._ [2021] H. B. Cai, X. X. Yan, P. L. Yao, and S. P. Zhu, Hybrid fluid-particle modeling of shock-driven hydrodynamic instabilities in a plasma, Matter and Radiation at Extremes 6, 035901 (2021).
  * Thoma  _et al._ [2017] C. Thoma, D. R. Welch, R. E. Clark, _et al._ , Hybrid-pic modeling of laser-plasma interactions and hot electron generation in gold hohlraum walls, Physics of Plasmas 24, 062707 (2017).
  * Wu  _et al._ [2017] D. Wu, X. T. He, W. Yu, and S. Fritzsche, Monte carlo approach to calculate proton stopping in warm dense matter within particle-in-cell simulations, Physical Review E 95, 023207 (2017).
  * Wu  _et al._ [2020] D. Wu, W. Yu, S. Fritzsche, and X. T. He, Particle-in-cell simulation method for macroscopic degenerate plasmas, Physical Review E 102, 033312 (2020).
  * Wu and Zhang [2023] D. Wu and J. Zhang, Head-on collision of large-scale high density plasmas jets: A first-principle kinetic simulation approach, Physics of Plasmas 30, 072711 (2023).
  * Wu  _et al._ [2021] D. Wu, Z. M. Sheng, W. Yu, S. Fritzsche, and X. T. He, A pairwise nuclear fusion algorithm for particle-in-cell simulations: Weighted particles at relativistic energies, AIP Advances 11, 075003 (2021).
  * Wu  _et al._ [2019] D. Wu, W. Yu, S. Fritzsche, and X. T. He, High-order implicit particle-in-cell method for plasma simulations at solid densities, Physical Review E 100, 013207 (2019).
  * Atzeni and Meyer-ter Vehn [2004] S. Atzeni and J. Meyer-ter Vehn,  _The physics of inertial fusion: beam plasma interaction, hydrodynamics, hot dense matter_ , Vol. 125 (OUP Oxford, 2004) pp. 134–194.
  * Sagdeev and Kennel [1991] R. Z. Sagdeev and C. F. Kennel, Collisionless shock waves, Scientific American 264, 106 (1991).
  * Balogh and Treumann [2013] A. Balogh and R. A. Treumann,  _Physics of Collisionless Shocks: Space Plasma Shock Waves_ (Springer New York, 2013) pp. 1–500.
  * Glenzer  _et al._ [2010] S. H. Glenzer, B. J. MacGowan, P. Michel, _et al._ , Symmetrie inertial confinement fusion implosions at ultra-high laser energies, Science 327, 1228 (2010).
  * Brysk [1973] H. Brysk, Fusion neutron energies and spectra, Plasma Physics 15, 611 (1973).
  * Sorasio  _et al._ [2006] G. Sorasio, M. Marti, R. Fonseca, and L. O. Silva, Very high mach-number electrostatic shocks in collisionless plasmas, Physical Review Letters 96, 10.1103/PhysRevLett.96.045005 (2006).
  * Fiuza  _et al._ [2013] F. Fiuza, A. Stockem, E. Boella, R. A. Fonseca, L. O. Silva, D. Haberberger, S. Tochitsky, W. B. Mori, and C. Joshi, Ion acceleration from laser-driven electrostatic shocks, Physics of Plasmas 20, 10.1063/1.4801526 (2013).


# Towards a relevant set of state variables to describe static granular
packings

Luis A. Pugnaloni<EMAIL_ADDRESS>Instituto de Física de Líquidos y
Sistemas Biológicos (CONICET La Plata, UNLP), cc. 565, 1900 La Plata,
Argentina.  Iván Sánchez  Departamento de Física y Matemática Aplicada,
Facultad de Ciencias, Universidad de Navarra, Irunlarrea S/N, 31080 Pamplona,
Spain.  Centro de Física, Instituto Venezolano de Investigaciones Científicas,
Apartado Postal 21827, 1020-A Caracas, Venezuela.  Paula A. Gago  Instituto de
Física de Líquidos y Sistemas Biológicos (CONICET La Plata, UNLP), cc. 565,
1900 La Plata, Argentina.  José Damas  Departamento de Física y Matemática
Aplicada, Facultad de Ciencias, Universidad de Navarra, Irunlarrea S/N, 31080
Pamplona, Spain.  Iker Zuriguel  Departamento de Física y Matemática Aplicada,
Facultad de Ciencias, Universidad de Navarra, Irunlarrea S/N, 31080 Pamplona,
Spain.  Diego Maza<EMAIL_ADDRESS>Departamento de Física y Matemática
Aplicada, Facultad de Ciencias, Universidad de Navarra, Irunlarrea S/N, 31080
Pamplona, Spain.

###### Abstract

We analyze, experimentally and numerically, the steady states, obtained by
tapping, of a 2D granular layer. Contrary to the usual assumption, we show
that the reversible (steady state branch) of the density–acceleration curve is
nonmonotonous. Accordingly, steady states with the same mean volume can be
reached by tapping the system with very different intensities. Simulations of
dissipative frictional disks show that equal volume steady states have
different values of the force moment tensor. Additionally, we find that steady
states of equal stress can be obtained by changing the duration of the taps;
however, these states present distinct mean volumes. These results confirm
previous speculations that the volume and the force moment tensor are both
needed to describe univocally equilibrium states in static granular
assemblies.

_Introduction:_ Finding out the appropriate set of macroscopic variables that
characterizes the equilibrium state of a given system is the first step in any
thermodynamic study callen . What these variables are for a sample of sand in
a box is still under discussion leiden . Twenty years ago, Edwards and
Oakshottedwards  put forward the idea that the number of grains $N$ and the
volume $V$ are the basic state variables that suffice to characterize a static
granular sample in equilibrium. The $NV$ granular ensemble was introduced as a
collection of microstates, where the sample is in mechanical equilibrium,
compatible with $N$ and $V$.

Experimentally, such equilibrium states are commonly obtained by tapping the
system at different dimensionless accelerations, $\Gamma$, which is the
standard parameter used to quantify the intensity of the external excitation
nowak ; nowak2 ; richard ; schroter ; ribiere . Hence, under the action of
repeated excitations at constant $\Gamma$, the system properties reach a
steady state where all observables fluctuate around well defined mean values.
Therefore, tapping provides the energy necessary to allow the system to
explore different configurations over which macroscopic observables can be
averaged. In the steady state, the mean volume $V$, or, more commonly, the
mean packing fraction, $\phi$ (defined as the percentage of the space occupied
by the particles), can be measured. It is now well established nowak ; nowak2
; richard ; schroter ; ribiere  that different steady states can be reached
reversibly by varying $\Gamma$. These experiments also suggest a monotonic
relation between $\Gamma$ and $\phi$. Nevertheless, $\phi$ does not only
depend on $\Gamma$. In harmonic pulses, if $\Gamma$ is kept constant while the
duration and the amplitude of the pulse are simultaneously modified, the
system evolves to different steady states. This fact was used to show that
equivalent equilibrium states pica ; schroter  can be generated by using
different combinations of pulse amplitude and duration. On the other hand, new
theoretical works snoeijer ; edwards2 ; blumenfeld ; henkes ; henkes2 ; tighe
suggest that the force moment tensor, $\Sigma$, ($\Sigma=V\sigma$, where
$\sigma$ is the stress tensor) must be added to the set of extensive
macroscopic variables (i.e., an $NV\Sigma$ ensemble) to describe adequately a
packing of grains.

Recent numerical simulations pugnaloni ; gago  have predicted also the
existence of a nonmonotonic relation between $\Gamma$ and $\phi$, if
sufficiently large values of $\Gamma$ are considered. A hint as to this
behavior can already be seen in some packings obtained in the laboratory
sibille ; although this has passed unnoticed by those studying equilibrium
properties. These results imply that equilibrium states of equal $V$ (or
$\phi$) could be obtained not only by changing the pulse duration, as in Ref.
schroter , but also by simply using different values of the excitation
parameter $\Gamma$. In what follows, we will show that these equal volume
equilibrium states are experimentally accessible. We will then use numerical
simulations to show that these states are, nevertheless, different. We have
found that two equilibrium states obtained with different excitations, even if
they have same $V$ and $N$, may present distinct values of $\Sigma$. This
result suggests the incomplete picture upon which research has been based for
over two decades when the properties of static granular matter have been
investigated; and takes us one step closer to completing the definition of the
state variables that are necessary to describe a static packing of grains.

Figure 1: (Color online)(a) Evolution of $\phi$ towards the steady state at
$\Gamma=7.0$ starting from the steady state corresponding to $\Gamma=1.5$. (b)
The evolution when switching back to $\Gamma=1.5$ movie1  (Arrows indicate the
tap where $\Gamma$ is switched). (c) Steady state packing fraction, $\phi$, as
a function of $\Gamma$ (green down triangles: steady state obtained in Exp A
through tapping at constant $\Gamma$ from an initial ordered state; orange up
triangles: steady state in Exp B obtained by an annealing protocol) nowak ;
movie1 . (d) Steady state relative packing fraction $\phi/\phi_{\rm{max}}$ as
a function of tapping intensity $\Gamma$ for the 3D cell ($\phi_{\rm{max}}$
corresponds to the value of $\phi$ obtained at the lowest $\Gamma$ studied).
Error bars represent the standard error of the mean estimated from the
standard deviation.

.

_Experimental results:_ A quasi 2D Plexiglass cell (width: 28 mm, height: 150
mm) was used to study the packing dynamics. The cell was filled with 1000
alumina oxide (Al2O3) spheres of diameter 1.000 $\pm$ 0.005 mm. The separation
between the Plexiglass sheets was 10% larger than the bead diameter in order
to minimize the particle-wall friction and prevent arching in the transversal
direction. The system was tapped with an electromagnetic shaker (Tiravib
52100) with a series of harmonic pulses of variable amplitude and constant
frequency $\nu=30$ Hz every three seconds. The tapping intensity was measured
with a piezoelectric accelerometer attached to the base of the cell. Although
in recent years alternative parameters have been proposed to quantify the
intensity of the taps dijksman ; ludewing , we employ the usual nondimensional
peak acceleration $\Gamma\equiv a_{peak}/g$ (where $g$ is the acceleration of
gravity). High resolution digital images of the packing were taken after each
tap. The center of each sphere was detected with an error of less than 2%. The
packing fraction of each image was obtained by considering each grain as a
disk of the corresponding effective diameter and then calculating the
percentage of the area covered by the disks in a rectangular area 10% smaller
than the size of the packing. Due to the space left between the Plexiglass
sheets, small overlaps can be observed in the 2D projection taken by the
photographs. This effect leads to estimations of $\phi$ above the 2D hexagonal
close packing value for the densest configurations studied. After a relative
short number of taps, the system reaches a steady state characterized by a
plateau in curves of $\phi$ vs. the number of taps. In contrast with the three
dimensional case, any stationary steady state is obtained after a few taps. We
checked that the mean value and standard deviation of the packing fraction
differ less than $0.1\,\%$ if $200$ or $2000$ taps of equilibration are
applied. As an example, in Fig. 1(a), we show that the steady state
corresponding to $\Gamma=7.0$ can be reached in a couple taps, even when the
initial packing is in an ordered configuration. Equilibrating the system back
at $\Gamma=1.5$ takes somewhat longer (about 100 taps) [see Fig. 1(b)] movie1
. Therefore, the sample initially prepared in a highly ordered configuration
was tapped 500 times for any given $\Gamma$ before taking averages over 100
taps, and $10$ independent runs were averaged. Alternatively, an annealing
protocol was used in which $\Gamma$ was decreased in discrete steps from the
highest values and tapped 500 times at each $\Gamma$ value without emptying
the cell.

In Fig. 1(c), we plot $\phi$ in the steady state as a function of $\Gamma$ for
different independent repetitions of the experiment. It can be seen that the
same results are obtained by equilibrating from an initial ordered structure
(Exp. A) or by following an annealing path (Exp. B). Although the fluctuations
of $\phi$ are large [see Fig, 1(a)], its mean value is well defined with a
small confidence interval [see error bars in Fig. 1(c)]. For low excitations,
$\phi$ decreases as $\Gamma$ is increased, in agreement with previous results
reported by several groups. However, beyond a certain value
$\Gamma_{\rm{min}}$, the packing fraction grows movie1 . The same trend is
observed if the tap frequency $\nu$ is changed (Exp. C) [see Fig 2(a)]. An
explanation for this behavior based on the formation of arches has been given
in pugnaloni .

Figure 2: (Color online). Results from the MD simulations of the soft disk
model for two different frequencies of the tapping pulse [blue squares:
$\nu_{1}=0.5(g/d)^{1/2}$, and red circles: $\nu_{2}=0.25(g/d)^{1/2}$]. (a)
Steady state packing fraction, $\phi$, as a function of the drive $\Gamma$. In
order to compare with the quasi 2D experiments (green down triangles, orange
up triangles and black stars) the vertical and horizontal axes have been
scaled with the characteristic values $\phi_{\rm{min}}$ and
$\Gamma_{\rm{min}}$, respectively. (b) Fluctuations, $\sigma_{\phi}$, of the
packing fraction in the steady state (measured by the standard deviation) as a
function of $\phi$. The arrows indicate the direction of increasing tapping
intensity. Error bars correspond to the standard error. (c) The trace,
Tr$(\Sigma)$, of the force moment tensor in the steady state as a function of
$\Gamma/\Gamma_{\rm{min}}$.

In order to assess whether the presence of a minimum in the $\Gamma$–$\phi$
curves is induced by the highly ordered crystal-like structures present in the
quasi 2D cell, we repeated the experiment with a 3D cell. This 3D setup
consists of a cell of the same height and width (as the quasi 2D cell) but 6
mm thick. In this case, the granular sample was made of polydisperse glass
beads 1.0 $\pm$ 0.2 mm in diameter. For the sake of simplicity, in the 3D
cell, the relative packing fraction was estimated from the height of the
granular layer. The stationary regime in the 3D setup was obtained after
$2\times 10^{4}$ taps. In Fig. 1(d) we show the steady state packing fraction
as a function of $\Gamma$ for this three-dimensional setup. Again, the same
trend as in the quasi 2D experiment is observed. This nonmonotonic dependence
challenges the idea that the volume defines the equilibrium state, unless
these equal volume states are proved to be equivalent mehta .

_Numerical evidences:_ In 2005, Edwards suggested that the stress tensor
should be included in the description of the equilibrium state of a static
granular sample along with $N$ and $V$ in the so-called full canonical
ensemble edwards2 . There is some consensus now snoeijer ; blumenfeld ; henkes
; henkes2 ; tighe  that the force moment tensor, $\Sigma\equiv V\sigma$ ,
allows for a microcanonical description henkes2 . However, no experiments or
simulations have shown that this variable is necessary to distinguish between
different states of thermodynamic equilibrium.

Since it is difficult to measure the force moment tensor in our experimental
setup, we use a realistic computer model of the quasi 2D setup pugnaloni ;
arevalo . We used a velocity Verlet algorithm to integrate the Newton
equations for 512 monosized disks in a rectangular box. The disk–disk and
disk–wall contact interaction comprises a linear spring–dashpot in the normal
direction and a tangential friction force that implements the Coulomb
criterion to switch between dynamic and static friction. Units are reduced
with the diameter of the disks, $d$, the disk mass, $m$, and the acceleration
of gravity, $g$. Details on the force equations and the interaction parameters
can be found elsewhere pugnaloni ; arevalo . Tapping is simulated by moving
the confining box in the vertical direction following a half sine wave
trajectory. The intensity of the excitation is controlled either through the
amplitude, $A$, or the frequency, $\nu$, of the sinusoidal trajectory; and it
is characterized by the parameter $\Gamma=A(2\pi\nu)^{2}/g$. We consider the
system has reached the steady state whenever a plateau in $\phi$ is observed
(with no visible trend in a plot of $\phi$ vs $\log(taps)$). Averages were
taken over 400 taps in the steady state and over 20 independent simulations
for each value of $A$ and $\nu$.

In Fig. 2(a), we show $\phi$ as a function of $\Gamma$ for this model.
Different values of $\Gamma$ were obtained by varying both $A$ and $\nu$. As
in the experiment, states of equal $\phi$ are generated at both sides of a
minimum packing fraction, $\phi_{\rm{min}}$. States of equal volume at each
side of $\phi_{\rm{min}}$ display slightly different volume fluctuations [Fig.
2(b)] suggesting that these states are not the same. In contrast to the
packing fraction, the trace of $\Sigma$ sigma  gently grows as $\Gamma$ is
increased [see Fig. 2(c)]. There is neither a minimum nor a maximum in the
$\Gamma$–$\Sigma$ curves. Therefore, states of equal $\phi$ at each side of
$\phi_{\rm{min}}$ present distinct $\Sigma$. Hence, states of equal volume at
each side of $\Gamma_{\rm{min}}$ are, in fact, different. This finding proves
that the volume alone cannot characterize the equilibrium states of a static
granular sample and suggests that the force moment tensor is a convenient
extra state variable.

Figure 3: (Color online). The loci of the generated equilibrium states in the
$\phi$–$\Sigma$ phase diagram. (a) The stationary packing fraction as a
function of the trace Tr$(\Sigma)$ for $\nu_{1}$ (blue squares) and $\nu_{2}$
(red circles). (b) The fluctuations, $\sigma_{\rm{Tr}(\Sigma)}$, versus
Tr$(\Sigma)$. Error bars as in Fig. 2. The shaded areas indicate the position
of $\phi_{\rm{min}}$.

In Fig. 3(a), the loci defined by the equilibrium states visited during the
numerical tapping experiments are plotted in a hypothetical $\phi$–$\Sigma$
thermodynamic phase space. It is clear that states of equal $\Sigma$ can
present different $\phi$ if prepared at different $\nu$. Thus, $\Sigma$ cannot
define the equilibrium state by itself. Both $V$ and $\Sigma$ have to be
specified to fully identify a given state. This result is confirmed by the
fluctuations of the force moment tensor, $\sigma_{\Sigma}$, [see Fig. 3(b)].
The fluctuations $\sigma_{\Sigma}$ for a given $\Sigma$ are not unique: they
depend on the frequency $\nu$ of the tapping. However, for the states obtained
with different excitations that display equal $\phi$ and $\Sigma$ (note that
such states would actually be a unique state), we find that $\sigma_{\phi}$
and $\sigma_{\Sigma}$ display also very similar values. Hence, the coincidence
of the mean and fluctuations of the state variables suggests that no other
extensive parameter would be necessary to describe such equilibrium states.

It is interesting to mention that the definition of an entropy $S(N,V,\Sigma)$
is particularly convenient since $\Sigma$ is not bounded as $V$ is. This
implies that $S$ can be a monotonic increasing function of $\Sigma$, which
permits a well behaved Legendre transformation to a canonical ensemble callen
. Hence, one can define a non-negative temperature-like quantity: the
angoricity edwards2 ; henkes2  —the proposed name for the inverse of the
conjugate variable to the force moment tensor. Let us finally recall that the
temperature-like quantity associated with the volume (so-called compactivity)
can present negative values due to population inversion brey .

_Conclusions:_ We have shown that steady states of static granular packings
obtained by tapping the system with different pulse strength and duration can
present different mean force moment tensor even if they correspond to the same
mean volume. Moreover, steady states that present the same mean force moment
tensor are distinguishable by their mean volume. To our knowledge, this is the
first experimental/numerical evidence that both extensive variables should be
included in an entropic formulation of the thermodynamics of such steady
states.

###### Acknowledgements.

This work was supported by project FIS2008-06034-C02-01 (Spain), PIUNA (Univ.
Navarra), CONICET (Argentina) and ANPCyT (Argentina).

## References

  * (1) H. B. Callen, Thermodynamics and an introduction to thermostatistics, 2nd Edition, Wiley-VCH, New York (1985).
  * (2) A special workshop on statistical mechanics of static granular media took place in Lorentz Center (Leiden, The Netherland) in July 2009. The main topic was the definition of state variables for equilibrium states in static samples. Presentations can be seen at http://www.lorentzcenter.nl/lc
/web/2009/340/info.php3?wsid=340

  * (3) S. F. Edwards and R. B. S. Oakeshott, Physica A 157, 1080 (1989).
  * (4) E. R. Nowak, J. B. Knight, M. L. Povinelli, H. M. Jeager and S. R. Nagel, Powder Tech. 94, 79 (1997).
  * (5) E. R. Nowak, J. B. Knight, E. Ben-Naim, H. M. Jaeger and S. R. Nagel, Phys. Rev. E 57, 1971 (1998).
  * (6) P. Richard, M. Nicodemi, R. Delannay, P. Ribière and D. Bideau, D, Nat. Mater. 4, 121 (2005).
  * (7) Ph. Ribière, P. Richard, P. Philippe, D. Bideau and R. Delannay, Eur. Phys. J. E 22, 249 (2007).
  * (8) M. Schröter, D. I. Goldman and H. L. Swinney, Phys. Rev. E 71, 030301(R) (2005).
  * (9) M. Pica Ciamarra, A. Coniglio and M. Nicodemi, Phys. Rev. Lett. 97, 158001 (2006).
  * (10) J. H. Snoeijer, T. J. H. Vlugt, W. G. Ellenbroek, M. van Hecke and J. M. J. van Leeuwen, Phys. Rev. E. 70, 061306 (2004).
  * (11) S. F. Edwards, Physica A 353, 114 (2005).
  * (12) R. Blumenfeld and S. F. Edwards, J. Phys. Chem. B 113, 3981 (2009).
  * (13) S. Henkes, C. S. O’Hern and B. Chakraborty, Phys. Rev. Lett. 99, 038002 (2007).
  * (14) S. Henkes and B. Chakraborty, Phys. Rev. E. 79, 061301 (2009).
  * (15) B. P. Tighe, A. R. T. van Eerd and T. J. H. Vlugt, Phys. Rev. Lett. 100, 238001 (2008).
  * (16) L. A. Pugnaloni, M. Mizrahi, C. M. Carlevaro and F. Vericat, Phys. Rev. E 78, 051305 (2008).
  * (17) P. A. Gago, N. E. Bueno and L. A. Pugnaloni, Gran. Matt. 11, 365 (2009).
  * (18) L. Sibille, T. Mullin and P. Poullain, Eur. Phys. Lett. 86, 44003 (2009).
  * (19) J. A. Dijksman and. M. van Hecke, Eur. Phys. Lett. 88, 44001 (2009).
  * (20) F. Ludewig, S. Dorbolo, T. Gilet, and N. Vandewalle, Eur. Phys. Lett. 84, 44001 (2008).
  * (21) EPAPS Document No. [number will be inserted by publisher]. Movie 1: A movie that shows the evolution of the packing during tapping at $\Gamma=7.0$ followed by tapping at $\Gamma=1.5$. Movie 2: A movie that shows the dynamics of a single tap at three very different values of $\Gamma$. For more information on EPAPS, see http://www.aip.org/pubservs/epaps.html.
  * (22) R. Arévalo, D. Maza and L. A. Pugnaloni, Phys. Rev. E 74, 021303 (2006).
  * (23) We have seen that the two invariants of $\Sigma$ (i.e. the trace and the determinant in 2D) are proportional to one another (i.e. $Det(\Sigma)=aTr(\Sigma)+b$). Moreover $a$ and $b$ are independent of the amplitude and duration of the taps used to reach a given $\Sigma$. Therefore, $Tr(\Sigma)$ suffices to characterize $\Sigma$ for our purposes.
  * (24) A. Mehta, Granular Physics, Cambridge University Press, Cambridge (2007), Chap. 13.
  * (25) J. J. Brey, A. Prados and B. Sanchez-Rey, Physica A 275, 310 (2000).