1 Introduction
The original SPH scheme, now called weakly compressible SPH (WCSPH), was developed in 1960s by R.Gingold and J. Monaghan [1]. The scheme used the leapfrog time integration and relied on an artificial term stabilizing the fluid in the presence of shocks. Many other stabilization strategies were conceived, including SPH [2] and methods based on Riemann solvers [3]. These approaches often provide results in reasonable agreement with experimental data or other numerical methods. Unfortunately, they often come at the cost of increased complexity and necessity to finetune additional numerical parameters. Failure in such finetuning can result in unrealistic dissipation, or large energy oscillations unless a convenient limiter is used [4].
There are essentially two different ways, how pressure can be calculated in WCSPH [5]
. A commonly used approach is to update the density iteratively, using the velocities obtained from the discretized balance of momentum. Alternatively, one can compute density directly as an SPH interpolation of particle masses. The latter approach is less common, but it leads to a symplectic structure of the system of SPH differential equations, which makes the WCSPH discretization completely stable even in the absence of viscosity or stabilization. Moreover, the WCSPH numerical scheme becomes globally reversible in time when also the fixedpoint arithmetic is used. For instance, the breaking dam simulation can be reversed to recover the initial conditions. The purpose of this paper is to show key advantages of the symplectic WCSPH compared to the standard approach.
However, as pointed out by D. Violeau, symplectic WCSPH shows unphysical behavior at free surfaces, where constant functions are not correctly interpolated. We suggest two ways how to correct this issue. The first method modifies the density by an appropriate integration constant, which is easy to implement. The second and more complex methods renormalizes the initial state by solving a certain nonlinear problem using Newton’s method. We also demonstrate that a careful implementation of the symplectic WCSPH scheme, using the corrected treatment of free surfaces and the fixedpoint (as opposed to the more usual floatingpoint) arithmetic leads to bitperfect global time reversibility. This can be considered an additional symmetrypreserving feature alongside the conservation of energy, mass, momentum and angular momentum.
Although the resulting code is energypreserving and reversible, it can still be considered dissipative in the sense that the Boltzmann entropy of the SPH particles grows. In other words, we observe the emergence of the second law of thermodynamics from purely reversible dynamics.
We demonstrate the various versions of WCSPH on dambreak and Gresho vortex benchmarks.
2 WCSPH for inviscid fluids
The standard weakly compressible spatial semidiscretization of Euler fluid equations in a uniform gravitational field is [8]:
(1a)  
(1b)  
(1c) 
where the pressure is a function of density via the baratropic formula,
(2) 
Constant stands for the referential density, is the gravitational acceleration, and is the numerical speed of sound (typically chosen as ten times the characteristic flow speed ). It is well known that the system of ordinary differential equations (1) conserves the total energy in the form
(3) 
where
Moreover, for , the momentum
and angular momentum
are conserved as well. Naturally, gravitational force can deliver some momentum to the fluid, while the equal opposite reaction of the fluid on Earth is considered grounded.
The WCSPH equations (1) represent a Hamiltonian system, which can be checked for instance using program [9]. But it is not a symplectic system because the mass density is treated as a state variable. If, on the other hand, the density would always be calculated from the current positions of the particles, then only the positions and momenta would be state variables and the system of equations would be symplectic. Symplecticity would make it easier to choose a geometric structurepreserving numerical integrator, for instance using the Verlet scheme.
A symplectic form of WCSPH is obtained by solving the ordinary differential equation for density (1a),
(4) 
instead of updating iteratively by 1a. However, this closed expression for density (4) is rarely used in practice because it leads to questionable behavior near free surfaces [5]. This happens because free boundary particles are underoccupied and, according to (4), they have . Therefore, they are equipped with large internal energy and, affected by its negative gradient, they start to vibrate. A remedy is provided in the following section, which leads to the possibility to use the symplectic form of WCSPH.
3 Treatment of free surfaces
3.1 Solving the equation for with an appropriate integration constant
The issue of unphysical behavior at free surfaces can be corrected by choosing an appropriate integration constant when solving the differential equation for ,
(5) 
where
Indeed, since WCSPH based on (1a) typically assumes zero internal forces at , it is actually (5) and not (4) that is equivalent to the standard WCSPH formulation, since (4) leads to nonzero density gradients near the boundary. Although schemes based on (1a) and (5) are completely equivalent at the ODE level, they become very much different once time discretization is considered. It appears that (1a) is currently used by most researchers, despite the advantages of Equation (5):

Inferring from positions using Equation (5) does not accumulate the density error.

The system of SPH evolution equations becomes symplectic (and not only Hamiltonian).

Therefore, symplectic structurepreserving integrators can be employed, which leads to simulations globally reversible in time.
However, even if we use the correct closed formula for density (5) and use a symplectic integrator, are facing another obstacle. For problems involving free surfaces, integration constant is typically much bigger for boundary particles and, therefore, whenever such particle submerges into interior of the fluid body, it generates a small void pocket around itself. This unwanted numerical artifact is removed in the following section.
3.2 Initial state correction (ISC)
The numerical artifact that particles that are originally at the boundary are treated differently than the rest of the particles, even if the submerge into the in the interior of the fluid, is rooted in the difference in the constants at the beginning. Therefore, w´e shall remove that artifact by correcting the initial state.
A solution is to correct the initial positions by a vector field
such thatThese nonlinear equations can be solved by a strategy based on the Newton’s method, details of which follow. Linearization with respect to leads to
On the left hand side, we identify the discrete divergence operator
(6a)  
see [5] for details. This linear problem for unknowns is underdetermined, so shall be restricted to in the form of a discrete gradient  
which can be also rewritten as  
(6b) 
Sparse linear system (6) for and has to be solved iteratively, updating and reevaluating according to (4), until with satisfactory precision. This initial state correction (ISC) procedure has been tested in our numerical examples, and the procedure converges quadratically except for problems with too much symmetry (for instance rectangle with particles arranged in square grid nearly always diverges)^{1}^{1}1For these cases, we recommend to initially disrupt particle positions by a random noise.. The whole procedure may be then summarized as follows:
Initial state correction (ISC) set randomly for all repeat until sufficient precision is acquired: for all recompute neighbor list find for all solve linear problem
This computation is expensive, but is performed only once per simulation. In place of the linear system with block matrix, one could also solve directly for :
However, the product leads to a matrix which is much more dense. A matrixfree iterative method could be used without evaluating the product directly, is a future possibility.
Let us note that after ISC, all constant fields will be correctly interpolated, at least initially. Weak compressibility helps to preserve this trait approximately in a way that does not accumulate error provided (4) is used. In this sense, ISC can be considered an alternative to (zeroth order) operator renormalization with the benefit that it does not break the symmetry of smoothing kernel. Unfortunately, this technique has two drawbacks. Firstly, it may slightly deform geometry in an unwanted way. Secondly, adapting it to a concrete setting may be difficult — for instance, when an anticlump term is included.
Both symplecticity of the system of equations and ISC improve the numerical solutions, but in order to obtain a globally reversible behavior we have to tackle a further problem caused at the walls of the simulation box, as in the following section.
3.3 Wall modelling
In WCSPH, walls can be modeled by a layer of dummy particles which behave as if they belonged to the fluid, except that their positions and velocities are fixed. In this approach, it is hoped that the forces preventing compression will deflect any incoming particle. This method, of course, violates balance of momentum, but this is plausible, since we can imagine that any momentum yielded by the fluid is grounded. More importantly though, dummy particles violate conservation of energy. To see this, consider a wave stopped by a wall. The fluid acts on the solid with equal and opposite force. Therefore, a wall particle is subject to acceleration, which should change its velocity after time step to
To preserve wall integrity, this velocity must be annulled, reducing the total energy by a small kinetic contribution. From this, we see that dummy particles are essentially dissipative.
One could argue that this is still physically reasonable because every fluidsolid collision creates a sound wave, leading to loss of energy. However, it is dubious whether such effect is correctly modeled by dummy particles. Moreover, it complicates efforts to make a strictly energyconservative scheme.
An alternative approach uses a fluidsolid force inspired by the LennardJones potential
by which a wall particle acts onto a fluid particle . Here,
and are numerical parameters [10]. Except for this interaction, wall particles are neither considered in the balance of mass (1b) nor in the calculation of density (5). Unlike the dummy particle approach, this solid wall modeling is conservative, with a modified energy
(7) 
This makes it potentially advantageous for problems involving inviscid flows.
In our case, when we seek a globally reversibleintime simulation, the conservativeness of the fluidwall interaction becomes important. However, even if we use the correct close formula for density (1a), the initial state correction, and the conservative fluidwall interaction, the symplectic Verlet scheme is not reversible globally in time because of the errors caused by the floatingpoint arithmetic. The following section contains the final ingredient necessary for that reversibility, the symplectic Verlet integrator with fixedpoint arithmetic.
4 Symplectic integrator and fixedpoint arithmetic
The symplectic character of evolution equations (1b) and (1c) is preserved in symplectic numerical schemes, for instance the classical Verlet scheme [11]. Consequence are for instance the longtime stability of the trajectories and conservation of of integrals of motion (energy, momentum, and angular momentum, among others). However, even application of the Verlet scheme in SPH does not lead to globally intime reversible simulations due to the floatingpoint errors, but a remedy is in the fixedpoint arithmetic, as we show below.
The usual discretization of the SPH equations using the Verlet scheme is
where is the acceleration composed of internal, gravitational, and wall forces. The scheme is of second order and symplectic, and therefore it conserves energy with error
that does not depend explicitly on time, provided that is sufficiently small [12]^{2}^{2}2Here, we rely on the fact that using the closed expression (5).. If we took as a separate state variable with its own evolution equation (1a) (usual WCSPH), symplecticity of the system of equations would be violated. The usual WCSPH equations are still Hamiltonian (generated by a Poisson bracket and a Hamiltonian), as checked by program [9], but they are not symplectic and they thus the (symplectic) Verlet scheme does not preserve their geometric structure, which causes energy to grow exponentially, unless some stabilization is used  see Figure 2.^{3}^{3}3In [8], it is suggested to correct this error by a midstep extrapolation of density. However, it should be noted this approach adds additional inaccuracy to the scheme.
Besides being symplectic, Verlet scheme is also time reversible, which means that after changing the sign of velocities, the scheme returns back in time. The scheme thus preserves the timereversal symmetry of the original system (1a), (1b), (1c). We can verify this by substituting
which produces the same set of equations
but in the reversed order. Unfortunately, the time reversibility fails in the floatingpoint arithmetic (FloPA), where addition is not associative. In particular, addition and subtraction of a value to a float does not recover . For instance
(8) 
Although this error is usually very small, it tends to accumulate in the presence of nonlinearities and destroys the reversibility in longer simulations.
There is, however, an easy solution suggested by [13] in the context of Nbody simulations of Solar System. It converts vectors , , and to the fixedpoint arithmetic (FixPA) just before they are used in the timestep evaluation. Since addition is associative in FixPA, this method allows for bitprecise time reversibility and longtime conservation of energy – see Figure 3. Note that the intermediate computations of the sums in (5) and (1b) can be still performed in FloPA and only then converted to FixPA. However, due to the nonassociativity of FloPA, summation of more than two elements in FloPA is orderdependent, and to ensure reversibility, the order of the summands must be chosen in a way that does not depend on the current state of the simulation (for instance, it is possible to perform summations always in the order of the particle indices).
density  1000 kg/  

spatial step  0.005 m  
num. sound speed  120 m/s  
gravitational constant  9.8 m/  
water column width  1m  
water column height  2m  
box width  4m  
box height  3m  
kernel support radius  
LennardJones radius  
LennardJones energy  
time step 
In summary, symplectic WCSPH with the correct closed formula for density, the initial state correction, the conservative wallfluid interaction, and the symplectic numerical Verlet scheme in the fixedpoint arithmetic finally lead to globallyintime reversible simulations. We have demonstrated the reversibility for instance on the dam break benchmark, and Appendix A contains the Gresho vortex benchmark. Although the simulation is symmetric, we can still observe the growth of entropy if we choose not see all the details of the simulation, which is the matter of the following section.
5 Emergence of the second law of thermodynamics
The second law of thermodynamics tells that the entropy of each isolated system grows until the system reaches the thermodynamic equilibrium. Where does this irreversible behavior come from when the underlying evolution equations for particle dynamics (here SPH) are reversible? In this section we illustrate that the emergence of the second law from completely reversible dynamics is caused by ignoring details of the dynamics, similarly as in [15].
The emergence of the second law is actually expected due to the result of Lanford [16] and following works [17], where a system of classical particles with a shortrange potential is shown to obey the Boltzmann equation. The collision term in the Boltzmann equation then causes the growth of the Boltzmann entropy. Figure 5 shows the growth of Boltzmann entropy in a reversible SPH breakingdam simulation.
(blue) in comparison with the equilibrium MaxwellBoltzmann distribution (orange) (
19) with fitted temperature . Right: Evolution of Boltzmann entropy (see Appendix B) in the reversible breakingdam simulation (blue). Despite the reversibility of the dynamics, the entropy increases and approaches equilibrium entropy corresponding to the temperature of the measured MaxwellBoltzmann distribution (blue line) from (21).But how to interpret the Lanford’s mathematical result and the observed growth of Boltzmann entropy of the SPH particles from the physical point of view? The SPH particles obey reversible Hamiltonian dynamics. Another possibility to describe their motion is the Liouville equation, which expresses reversible evolution of the Nparticle distribution function, . Indeed, Liouville equation follows from the particle mechanics, and, vice versa, if we set the distribution function to be the product of Dirac distributions in the positions and momenta of all particles, we recover Hamilton canonical equations. Liouville equation also conserves the Liouville entropy
(9) 
which thus remains constant. This is consistent with the reversibility of the underlying Hamiltonian particle mechanics.
Now we decide not to see all the positions and momenta of the individual particles, observing for instance only the probability distribution of momenta at given space and time. Such oneparticle distribution function
indeed does not contain the knowledge of positions and momenta of all particles and thus the entropy of the system expressed in terms of must be higher than the Liouville entropy (functional of ). The former is the Boltzmann entropy and it is obtained by maximization of the Liouville entropy subject to the constraint that the oneparticle distribution function is known [18]. Plugging the resulting back into the Liouville entropy, we obtain the Boltzmann entropy(10) 
see [18]. Because the Boltzmann entropy is the Liouville entropy evaluated at the point where it is maximal (keeping the knowledge of ), the Boltzmann entropy is higher than the Liouville entropy. Although the Liouville entropy remains constant in the dynamics, the Boltzmann entropy can grow in the dynamics following the Liouville equation.
The irreversibility of the dynamics of the oneparticle distribution function can be illustrated on the reversible breakingdam simulation with fixedpoint arithmetic. Section 4 contains details of the simulation, but let us also show the histogram of the particle velocity distribution in the middle of the simulation, just before the velocities are inverted. Figure 5 shows the measured histogram of the velocities of the SPH particles in comparison with the equilibrium MaxwellBoltzmann distribution function (19) obtained by fitting the effective temperature^{4}^{4}4Note that the effective temperature is not the physical temperature of the system because it is calculated from the overall motion of the macroscopic SPH particles. To see this, consider the extreme case of only one SPH particle with weight one kilogram, where the effective temperature obtained from the motion of that macroscopic particle does not coincide with the physical temperature of that particle.. Although all velocities were initially zero, their distribution approaches the equilibrium distribution. This approach towards equilibrium is reflected in the growth of Boltzmann entropy during the longtime reversible simulation.
In other words, when the particles are described by the Boltzmann equation, our knowledge is incomplete, since the positions and momenta of all the particles remain unknown, in contrast to the description by means of Hamilton canonical equations or Liouville equation. And this lack of knowledge makes the evolution of the oneparticle distribution function irreversible. The irreversibility is then explicitly expressed by the collision integral in the Boltzmann equation. The second law of thermodynamics thus emerges from completely reversible dynamics when our description is incomplete (not seeing all positions and momenta of the particles).
In summary, if we see all the positions and momenta in the SPH simulation, we can not see the second law of thermodynamics. Indeed, the simulation is reversible and the Liouville entropy remains constant. However, when we only focus on the oneparticle distribution function, we can see the growth of Boltzmann entropy and thus irreversible behavior. The second law of thermodynamics emerges when we our knowledge about the precise state of the system is incomplete [15].
Interestingly, since our discrete system is deterministic, reversible and can exist in only finitely many states, it is recurrent [19]. In other words, if the system evolves long enough, it comes back to the initial state exactly (provided that the particles are not allowed to escape to infinity). However, since the phase space of the system contains thousands of SPH particles, it has enormous amount of states, which makes it highly unlikely to observe the recurrence theorem in practice.
6 Conclusion
In this paper we have turned the usual weakly compressible smoothed particle hydrodynamics (WCSPH) to a symplectic form by finding the correct closed formula for the mass density. Because the closed formula for density leads to different treatment of particles that are initially near the boundary, the method of initial state correction (ICS) was introduced, due to which all particles are then treated in the same way. This leads to stable SPH simulations in the presence of free surfaces without any further stabilization.
In order to get simulations that preserve the energy, we use a conservative fluidwall interaction. A symplectic (Verlet) scheme, which is suitable for the symplectic WCSPH, implemented in the fixedpoint arithmetic then leads to globallyintime reversible SPH simulations. This is demonstrated on the dam break benchmark, where inversion of velocities at a later stage of the simulation eventually leads the system back to the initial state. The simulations are available in a new Julia package SmoothedParticles.jl [6].
Despite the reversibility of the simulations, we observe thermodynamic behavior when we do not use all the details of the simulation. The Boltzmann entropy, which depends only on the oneparticle distribution function and not on positions of individual particles, grows in the dambreak simulation and approaches the equilibrium value. In other words, we observe the emergence of the second law of thermodynamics from purely reversible dynamics, caused by reduction of our knowledge about the system.
In future, we would like to extend the WCSPH framework to nonisothermal fluids and solids while keeping the Hamiltonianity of the equations and numerical schemes.
References
 [1] R. Gingold and J. Monaghan, “Smoothed particle hydrodynamics: theory and application to nonspherical stars,” Mon. Not. R. Astron. Soc., vol. 181, no. 3, pp. 375–389, 1977.
 [2] M. Antuono, A. Colagrossi, S. Marrone, and D. Molteni, “Freesurface flows solved by means of sph schemes with numerical diffusive terms,” Computer Physics Communications, vol. 181, no. 3, pp. 532–549, 2010.
 [3] C. Zhang, X. Y. Hu, and N. A. Adams, “A weakly compressible sph method based on a lowdissipation riemann solver,” J. Comput. Phys., vol. 335, pp. 605–620, 2017.
 [4] Z.F. Meng, A.M. Zhang, P.P. Wang, and F.R. Ming, “A shockcapturing scheme with a novel limiter for compressible flows solved by smoothed particle hydrodynamics,” Computer Methods in Applied Mechanics and Engineering, vol. 386, p. 114082, 2021.
 [5] D. Violeau, Fluid Mechanics and the SPH Method: Theory and Applications. Oxford University Press, Oxford, UK, 2012.
 [6] O. Kincl and M. Pavelka, “SmoothedParticles.jl,” 2021, https://github.com/OndrejKincl/SmoothedParticles.jl
 [7] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, “Julia: A fresh approach to numerical computing,” SIAM Review, vol. 59, no. 1, pp. 65–98, 2017.
 [8] J. Monaghan, “Smoothed particle hydrodynamics,” Reports on Progress in Physics, vol. 68, p. 1703, 2005.
 [9] M. Kroeger and M. Huetter, “Automated symbolic calculations in nonequilibrium thermodynamics,” Comput. Phys. Commun., vol. 181, p. 2149–2157, 2010.
 [10] J. J. Monaghan, “Simulating free surface flows with sph,” Journal of computational physics, vol. 110, no. 2, pp. 399–406, 1994.
 [11] B. Leimkuhler, S. Reich, and C. U. Press, Simulating Hamiltonian Dynamics. Cambridge Monographs on Applie, Cambridge University Press, 2004.
 [12] J. Candy and W. Rozmus, “A symplectic integration algorithm for separable hamiltonian functions,” Journal of Computational Physics, vol. 92, no. 1, pp. 230 – 256, 1991.
 [13] H. Rein and D. Tamayo, “JANUS: a bitwise reversible integrator for Nbody dynamics,” Monthly Notices of the Royal Astronomical Society, vol. 473, pp. 3351–3357, 2017.
 [14] S. Koshizuka and Y. Oka, “Movingparticle semiimplicit method for fragmentation of incompressible fluid,” Nuclear science and engineering, vol. 123, no. 3, pp. 421–434, 1996.
 [15] M. Pavelka, V. Klika, and M. Grmela, “Generalization of the dynamical lackoffit reduction,” Journal of Statistical Physics, vol. 181, no. 1, pp. 19–52, 2020.
 [16] O. E. Lanford, “Time evolution of large classical systems,” in Dynamical systems, theory and applications (J. Moser, ed.), vol. 38 of Lecture Notes in Physics, pp. 1–111, Springer–Verlag, Berlin, 1975.
 [17] M. Pulvirenti, C. Saffirio, and S. Simonella, “On the validity of the boltzmann equation for short range potentials,” Reviews in Mathematical Physics, vol. 26, no. 02, p. 1450001, 2014.
 [18] M. Pavelka, V. Klika, and M. Grmela, Multiscale ThermoDynamics. Berlin: de Gruyter, 2018.
 [19] D. Wallace, “Recurrence theorems: a unified account,” Journal of Mathematical Physics, vol. 56, no. 2, p. 022105, 2015.
 [20] H. GrimmStrele, F. Kupka, and H. Muthsam, “Curvilinear grids for weno methods in astrophysical simulations,” Computer Physics Communications, vol. 185, no. 3, pp. 764–776, 2014.
Appendix A Gresho vortex benchmark
The Gresho vortex benchmark [20] prescribes a tangential velocity component of Eulerian fluid in polar coordinates as
The vortex is confined in a box with noslip wall (though different variants can be found in literature). Theoretically, in the absence of viscosity, the vortex should be stationary, but obtaining this in a numerical scheme is challenging. In our numerical experiment, we simulate the flow in timeinterval and measure the error using a discretized Sobolev norm
(11) 
The meaning of the factor is that the zero velocity field corresponds to error . In order to prevent formation of void space in the center of vortex, we had to include an anticlump force term in the form
no filter  passive filter  active filter  

square  28.16%  12.20%  53.07% 
hexagonal  27.84%  11.77%  55.79% 
Vogel  29.85%  13.31%  78.06% 
Vogel + ISC  28.19%  11.26%  55.25% 
density  1  

spatial step  1e2  
num. sound speed  20  
kernel support radius  
box size  1  
anticlump pressure  10  
time step  
filter frequency  30 
where is a smoothing kernel with support radius . Additionally, we consider noise reduction by the Shepard filter
where
This filter can be applied either in postprocessing (passive filter), or by setting every time steps (active filter). Results are summarized in Table 2 and Figure 6. Simulation parameters are listed in Table 3.
Appendix B The MaxwellBoltzmann entropy
In this section we derive the formula for Boltzmann entropy in terms of the MaxwellBoltzmann distribution and, subsequently, to find its equilibrium value.
b.1 The reduced Boltzmann entropy in terms of the MaxwellBoltzmann distribution
Let us start with the Boltzmann entropy in two dimensions expressed in terms of the oneparticle distribution function ,
(12) 
where the position is constrained to a box with volume and where
is the Planck constant. This formula can be obtained for instance by the principle of maximum entropy from the Liouville entropy, from where it also follows that the oneparticle distribution function is normalized to the number of particles, see
[18].Assuming that the distribution function depends only on the magnitude of the momentum (isotropic dependence), the Boltzmann entropy can be rewritten as
(13) 
where is the norm of the momentum. This dependence on the oneparticle distribution function has to be converted to a dependence on the MaxwellBoltzmann distribution function
, which tells the probability that the norm of velocity of a particle
is in the interval . This is done using the normalization,(14) 
which leads to the expression of the twodimensional MaxwellBoltzmann distribution function in terms of the oneparticle distribution,
(15) 
which is normalized to unity. Note that we assume that the oneparticle distribution is homogeneous in space (the total volume being ) and isotropic in momentum. These assumptions are valid in the thermodynamic equilibrium, but they approximately hold also before the equilibrium is reached, and for the purpose of showing that the Boltzmann entropy grow in our simulations, the approximation is satisfactory. Finally, the Boltzmann entropy in terms of the MaxwellBoltzmann distribution function becomes
(16) 
Since the number of particles is conserved in our simulation, we only evaluate the part of Boltzmann entropy that varies when changes,
(17) 
called reduced Boltzmann entropy. Note that the integral converges near the origin because the MaxwellBoltzmann distribution typically has linear growth there.
In our simulations, this reduced Boltzmann entropy is numerically integrated using the approximate MaxwellBoltzmann distribution function obtained from a histogram of the particle velocities at given time instant.
b.2 Equilibrium Boltzmann entropy
What is the final value of the reduced Boltzmann entropy when the system of particles reaches the thermodynamic equilibrium? This question can be answered in two steps. First, the equilibrium distribution function is calculated by the principle of maximum entropy (MaxEnt). Second, the equilibrium distribution function is plugged into the formula for the reduced Boltzmann entropy. Note, however, that the MaxEnt step depends on the energy (Hamiltonian) of the system and it depends on the complexity of the Hamiltonian whether the calculation can proceed purely analytically (without numerical solutions). Therefore, we use the simplest Hamiltonian approximating the true Hamiltonian of our SPH particles, consisting only of the kinetic energy of the particles. Maximization of the Boltzmann entropy subject to the constraints given by the total energy and total number of particles leads to the equilibrium oneparticle distribution function
(18) 
where and are the Lagrange multipliers corresponding to the two constraints (number of particles and total energy). From the normalization it follows that , while the other Lagrange multiplier can be interpreted as the inverse temperature, .
The equilibrium MaxwellBoltzmann distribution function then becomes
(19) 
As time proceeds in our simulations, the histogram of particle velocities approaches the equilibrium MaxwellBoltzmann distribution.
When this equilibrium distribution function is plugged back into the Boltzmann entropy, the reduced Boltzmann entropy becomes
(20) 
This equilibrium entropy, which depends on , can be calculated once the temperature is obtained by fitting the histogram of particle velocities to the equilibrium MaxwellBoltzmann distribution function (19). This makes sense, however, only at later stages in the simulations, when the histogram approaches the equilibrium distribution.
Instead of using the equilibrium temperature, which makes sense only in later stages of the simulations, we can express the equilibrium entropy in terms of energy, which can be measured anytime. The total kinetic energy of the twodimensional system of particles is equal to , which makes it possible to write the equilibrium entropy in terms of the energy,
(21) 
This value of equilibrium entropy is close to the Boltzmann entropy obtained directly from the approximated MaxwellBoltzmann distribution, while the entropy based on the equilibrium temperature is only approached at later stages of the simulations.
Acknowledgments
OK was supported by project No. START/SCI/053 of Charles University Research program. MP was supported by project No. UNCE/SCI/023 of Charles University Research program. OK and MP were also supported by the Czech Science Foundation (project no. 2022092S)