Work and Thermal Fluctuations in Crystal Indentation under Deterministic and Stochastic Thermostats: The Role of System–Bath Coupling

The Jarzynski equality (JE) was originally derived under the deterministic Hamiltonian formalism, and later, it was demonstrated that stochastic Langevin dynamics also lead to the JE. However, the JE has been verified mainly in small, low-dimensional systems described by Langevin dynamics. Although the two theoretical derivations apparently lead to the same expression, we illustrate that they describe fundamentally different experimental conditions. While the Hamiltonian framework assumes that the thermal bath producing the initial canonical equilibrium switches off for the duration of the work process, the Langevin bath effectively acts on the system. Moreover, the former considers an environment with which the system may interact, whereas the latter does not. In this study, we investigate the effect of the bath on the measurable quantity of the JE through molecular dynamics simulations of crystal nanoindentation employing deterministic and stochastic thermostats. Our analysis shows that the distributions of the kinetic energy and the mechanical work produced during the indentation processes are affected by the interaction between the system and the thermostat baths. As a result, the type of thermostatting has also a clear effect on the left-hand side of the JE, which enables the estimation of the free-energy difference characterizing the process.


Introduction
The Jarzynski equality (JE) is meant to link the statistics of non-equilibrium works to an equilibrium property: the free-energy difference between two equilibrium states of a given physical system [1]. In this sense, the JE complements the fluctuation-dissipation relationships that obtain non-equilibrium properties from equilibrium experiments. The present literature is vast and substantially includes diverse investigations on the validity of the JE. Most authors report verifications of the theory or provide reasons in support of the validity of the JE; such as those discussing optimal protocols for stochastic systems; see, e.g., Refs. [2][3][4][5]. However, investigations considering a reduced number of work measurements find violations of the JE or express some concern in connection with the Hamiltonian derivation; see, e.g., Refs. [6][7][8][9]. In particular, the study of a variable volume system [9] shows that the Hamiltonian of the JE is not universal. Moreover, Ref. [10] demonstrates that the JE is violated even in small systems complying with the Jarzynski theory due to the emergence of process-dependent irreversibilities at the nano-scale.
In this study, we investigate the effect of deterministic and stochastic thermostats on the measurable (or computable) quantity that appears in the JE. We consider that this assessment is required because there are two fundamentally different derivations of the JE that apparently lead to the same work-fluctuation expression, when in fact they refer to different types of experimental conditions and, hence, to distinct free-energy quantities. To where λ is a time-dependent parameter controlled by an external agent, H S is the energy of S, H E is that of E, and h int is the interaction energy of S with E. Initially, λ(0) = α, and S+E is in equilibrium with a thermal bath (B) at temperature T; hence, the statistics of its phases Γ are given by the canonical ensemble where Q S+E (λ) = e −βH(Γ;λ) dΓ , and β = (k B T) −1 characterizes the thermal bath. At time t = 0, the S+E is disconnected from B, and the external agent acts on λ for a time τ, which varies from its initial value α to its final value λ(τ) = ω. This fully deterministic process is repeated many times, with identical λ(t), but each time starting from a different initial condition Γ 0 = (x 0 , y 0 ) and different energy H(Γ 0 ; α) dictated by the initial canonical distribution f α . The fact that S+E and B do not exchange energy during the process may be justified by assuming that only a negligible amount of energy can be exchanged in the typically short process time τ. Then, the derivation of the JE-that proceeds through exact analytical calculations-suggests that the work distributions of practically any process taking λ from the given initial value α to the final value ω, during any process time τ > 0, can be used to compute the free-energy difference between the canonical state at temperature T with Hamiltonian H(Γ; λ(0) = α) and the canonical state at temperature T with Hamiltonian H(Γ; λ(τ) = ω). The calculations [1] can be summarized as follows.
Letting S t λ : M → M denote the phase space evolution operator so that Γ ∈ M evolves into S t λ Γ ∈ M at time t and letting x t be the coordinates and momenta of S at time t, the quantity is introduced and called work [1]. Indeed, when λ represents a position in space and the derivative of H S with respect to λ represents a force, W J corresponds to a mechanical work. Following Ref. [1], one obtains that and where is the difference of the free energies of S+E in the canonical equilibrium at temperature T with parameters ω and α, respectively. To obtain the free-energy difference essentially related to S, the following quantity is introduced which was proposed by Kirkwood [11] to treat subsystems of macroscopic dense fluids in thermodynamic equilibrium. The quantity H * S (x; λ) is the energy of S plus an average contribution coming from the interaction of S with E. Then, the system whose Hamiltonian is H * S , S * say, can be associated to the canonical ensemble and to the free energy i.e., the free energy of a hypothetical system with Hamiltonian H * S . This quantity is then linked to the solvated free energy of "S in E" [12].
It is found that which readily leads to the JE where · α is the canonical average with respect to the initial ensemble f α of S+E. As noted in Ref. [12], transformations of the free energy of a system strongly interacting with another one are not particularly interesting; hence, some kind of solvated free energy is preferred. There are two aspects of this theory that seem to remain unscrutinized in the present literature. First, the energy exchange between S+E and B may be small, but so is the free-energy variation computed by means of the JE. Moreover, E can be absent in the Hamiltonian derivation or, equivalently, h int may vanish. This is the case of the Langevin derivation, outlined in the following section. However, there is a second major difference between the two derivations: B acts on S during the Langevin process, whereas in the Hamiltonian derivation B is not considered.

The Langevin Derivation of the JE
The JE can also be derived under a stochastic framework [13], particularly using the Langevin equation [2,4] where X t is the vector of position coordinates of S, U is a time-dependent potential with t ∈ [0, τ], and W t is a Wiener process. Then, is a quantity related to the variations of the system's energy produced by changes in λ.
In the particular case in which λ represents a position in space, w t becomes a mechanical work. Additionally, when λ defines the position of the external agent exerting a force F on the system, w t describes the work carried out by the force (−F) acting on the agent.
Because there is no environment, the relevant free energy is that of S alone, which can be called intrinsic free energy. At a given λ, this is defined by where and Ω is the set of coordinates x. The quantity G(λ) is intrinsic in the sense that it solely refers to the properties of the system S. Thus, this framework is analogous to that adopted in the deterministic Hamiltonian derivation assuming weak coupling (h int ≈ 0). However, unlike the Hamiltonian derivation, the stochastic Langevin formulation allows S to interact with B assuming that the kinetic energy plays no explicit role in the work statistics.
To obtain the JE, one can start from the joint random variable (X t , w t ) and its probability density p t (x, a), where a ∈ R. Then, the quantity of interest with initial condition obeys Because the average of the exponential of −w t with respect to the initial distribution is expressed through by taking λ(0) = α and λ(τ) = ω, the JE is readily obtained Analogous to the Hamiltonian formulation, the JE is here derived irrespective of the form of λ(t), as long as λ(0) = α and λ(τ) = ω, for any process time τ. Note that Equation (20) is the JE for the intrinsic-not the solvated-free-energy difference of S. Thus, when E is absent (in the Hamiltonian formulation), both derivations refer to the same intrinsic quantity.

Comparison between the Two Derivations
The above summaries show that the Hamiltonian and Langevin theoretical derivations refer to two different classes of experiments. The former refers to systems that do not interact with a heat bath for the duration of the work process, whereas the latter to systems that do interact with a thermal bath. The described experiments differ also in assuming the presence or the absence of a third object, called environment (here, denoted as E). The conditions in which E acts on S can be related to protein stretching experiments [14], where S can be represented by the protein, E by the water in which the protein is immersed, and B can be the air of the laboratory in which the water pool is situated. The Langevin setting may also involve an environment E. However, the interaction of E with S is reduced due to the viscosity damping applied to S and not to an effective interaction energy (h int from Equation (1)). This is the case of the harmonic oscillators reviewed-along with other experimental settings-in Ref. [15].
In both derivations, the result is process independent. The protocol independence is less obvious under the Hamiltonian framework, where rapid transformations may not allow an efficient exchange of energy between S and B, and the number of particles as well as space and time scales strongly influence the system's response to external perturbations.
We consider that the role of system-bath coupling in the work statistics-which is the measurable quantity in the JE-requires further evaluation. Assuming a condition in which E is absent, both Hamiltonian and Langevin derivations of the JE lead to the intrinsic freeenergy variation and follow the same expression. The fact that nonequilibrium processes typically last only a short time supports that the exchanged energy may be negligible. However, this is a delicate assumption in at least two regards: (1) the values of the freeenergy difference are usually small as well, and (2) the process independence convoked by both approaches allows the process time to be long.
To investigate such questions, we consider the molecular dynamics model of crystal nanoindentation described in Ref. [10] using Nosé-Hoover and Langevin thermostats. Our analysis shows that the the coupling between S and B has a clear effect on the resulting distributions of the kinetic energy and the mechanical work. Consequently, the left-hand side of the JE is affected. Interestingly, the work statistics obtained in indentations protocols that produce fully reversible elastic deformations in the crystal lead to different distributions as a function of the thermostat.

Molecular Dynamics Setup: Crystal Nanoindentation
We carry out an extensive number of molecular dynamics (MD) simulations to investigate the properties of a minute (001)-oriented Ta crystal indented by a spherical nanoscopic tip. To run the simulations, we employ the open-source LAMMPS code [16]. The interatomic force field is modeled by means of the embedded-atom method potential developed by Ravelo et al. [17].
The Ta crystal has a cuboid shape of size 10.8 × 10.8 × 6 nm 3 and contains 40,293 particles. In the MD domain, the crystal's particles are sorted into two groups: the particles contained in the two lowermost atomic planes of the crystal create an effective floor that prevents the downward displacement of the crystal during indentation as their motion is restricted. The remaining particles constitute the hereafter called particle system, which includes a total number of N = 39,204 particles whose positions, q ≡ {q 1 ; q 2 ; . . . ; q N }, and velocities, v ≡ {v 1 ; v 2 ; . . . ; v N }, are, respectively, denoted by . . , N. The system's particles are free to move and interact with each other according to the prescribed ensemble properties and the interatomic potential. Periodic boundaries are applied to the lateral sides of the crystal.
The indenter is modeled by a spherical-shaped repulsive potential, of radius R = 3 nm and center C with coordinates q c = (x c , y c , z c ), where k is the indenter stiffness and Note that the indenter acts on particle i when δ i ≤ 0. Then, the repulsive force exerted on particle i, Here, k is set to 100 eV/Å 3 . Notice in Equation (22) that the direction of the repulsive force from the indenter is dictated by the (q i − q c ) vector. Figure 1a depicts the computational domain in our MD simulations, which contains the Ta crystal and the repulsive indenter. This computational approach has been largely employed in MD studies of indentation in metallic bodies (cf. Ref. [18] and references therein). We create an initial state in which the particles occupy the lattice positions of the Ta crystal and the velocities are taken from a normal distribution with 0 mean and a standard deviation scaled to produce a temperature T = 300 K. To generate an equilibrium canonical distribution at T, we run a preliminary 20 ps thermalization during which the particles follow NVT conditions with the Nosé-Hoover (NH) thermostat [19] controlling the system's temperature at T = 300 K. Thus, the positions and velocities of the particles at t = 0 (i.e., prior to indentation) are sampled from the canonical distribution produced during the NVT thermalization run.
The indentation run consists of a closed loading/unloading loop where the indenter moves vertically with constant velocity, v ind = 10 m/s. The motion of the vertical coordinate of the indenter center, z c , is described through where h(t) maps the penetration of the indenter during τ. Note that the maximum penetration, h max , is attained at t = τ/2, where z c (t = τ/2) = z 0 − h max and τ = 2h max /v ind ; see Figure 1b. We conveniently define z 0 such that the indenter tip and the top surface of the crystal are separated by a small vertical distance ∆ = 0.5 Å, so that the constituting particles are guaranteed to lie outside of the radius of action of Φ at t = 0 and t = τ; see Figure 1a. This also allows the particles to arrange initially into an unperturbed Ta bcc crystalline configuration during the thermalization run. The applied indentation load, P, is then defined as the sum of the vertical, repulsive force contribution, P = − ∑ N i=1 F iz , coming from the particles that satisfy δ i ≤ 0; see Equation (22). The computational timestep, dt, is set to 2 fs in all of the MD simulations.
In our indentation setting, the particles of the crystal constitute the system of interest S, the indenter-whose time dependent potential energy appears in the Hamiltonian of S-represents the external agent, and no environment is included, so h int = 0. The Hamiltonian of S is then given by H(q, p, q c ) = K(p) where q is the vector of position coordinates; p = mv is the vector of momenta of all particles, with m being the mass of a Ta atom and v being the set of their velocities; K(p) is the kinetic energy of the system (Section 2.2.2); and V(q) is the potential energy of the interatomic interactions. Note that, here, the time-dependent parameter of the Jarzynski theory is λ(t) = q c (t), i.e., the moving center of the indenter potential Equation (21).
To evaluate the effect of the coupling between the system and the bath on the mechanical and kinetic response of the system to the indentation processes, we adopt the following computational approaches.
(1) We perform MD indentation simulations using the (deterministic) NH thermostat with 3 NH chains [19] to implement a condition of constant number of particles N, volume V, and temperature (that represents the thermal bath at T = 300 K). To tune the coupling of the particles with the NH bath, we vary the NH thermostat parameter ω p that accounts for the frequency at which the particles are thermostatted; see the discussion given in Supplementary Section S1. With ω p = 100 dt, the energy exchange between the system and the NH bath is sensibly strong despite the short time of our indentation processes; see Supplementary Figure S1. (2) For ω p = 100,000 dt, the sluggishness of the heat flow between the system's particles and the NH bath describes similar conditions to those considered in the Jarzynski theory, which neglects the system-bath coupling. (3) By removing the thermostat-i.e., the ω p = ∞ limit-we obtain an adiabatic evolution of the system with unthermostatted particles. These conditions emulate those considered in the Hamiltonian derivation of the JE. (4) Lastly, we use a stochastic Langevin thermostat at T = 300 K, which acts on the system via a random force [20]. We impose a damping coefficient of γ L = 1 ps −1 that allows an efficient energy exchange between the system and the Langevin bath. For further details, see Supplementary Section S2. This approach reproduces the scheme adopted in the Langevin derivation of the JE.

Computation of Thermodynamic and Mechanical Properties
Classical MD provides a window into the microscopic dynamical behavior of the constituent particles of the system. Thus, MD simulations-unlike experiments-give access to all the necessary dynamical ingredients of particle systems, which enables the computation of equilibrium macroscopic properties by sampling from a statistical mechanical ensemble.
In this investigation, special attention is given to two specific quantities measured during the indentation runs (during which Φ effectively acts on the system): (i) the mechanical work exerted to the indenter by the particles and (ii) the total kinetic energy of the system's particles. In the following, we explain the post-simulation and on-the-fly algorithms that we employ to access to these quantities.

The Mechanical Work
The elementary mechanical work carried out on the indenter by the system, involves the infinitesimal displacements of the indenter, dq c = (0, 0, dz c ) and the opposite vertical forces, −F iz , where η is the step function, so η(δ i ) = 1 for δ i ≤ 0 and η(δ i ) = 0 for δ i > 0. (As discussed in Ref. [10], the elementary mechanical work carried out by the indenter on the system, dW I = ∑ N i=1 F i · dq i , differs not only in the sign from Equation (25) but also substantially. Thus, an external operator cannot obtain the work carried out on the system W I from external measurements of dW S , particularly in small systems made of classical particles). Note that the negative sign in Equation (25) comes from the force that particle i exerts on the indenter as derived from the action-reaction principle. Because −F iz ≥ 0, the term dz c dictates the sign of dW S . Additionally, note that dz c = −v ind dt during the loading stage and dz c = +v ind dt during unloading; see Equation (23).
We use an in-house post-simulation code to obtain dW S (q, q c ; t) by means of Equation (25). This equation is evaluated under evenly spaced infinitesimal time intervals of δt = 0.05 ps that corresponds to a net (infinitesimal) variation in the vertical indenter motion of 5 pm (=0.05 Å). Finally, the total mechanical work W S carried out during a time interval from t = 0 to t = τ is calculated as the sum of the computed elementary works from Equation (25), where l = τ/δt. In our MD indentation setup, W S (τ) = −W J (τ) [10].

The Total Kinetic Energy
We also obtain on-the-fly values of the total kinetic energy of the particle system, where K i is the instantaneous, translational kinetic energy of particle i.
To investigate the effect of the system-bath coupling on the system's thermodynamic properties during indentation, we assess the kinetic fluctuations of the system in terms of the statistical behavior of the quantity K tot . In our indentation simulations, K tot is evaluated every 0.05 ps (=25 dt). Since we perform swift nonequilibrium processes and our MD simulations refer to a relatively small number of particles (N = O(10 4 )), the consequences of the thermodynamic limit can be observed. In particular, note that K tot might not coincide with 3Nk B T/2, where T is the bath's temperature. When K tot is averaged over the ensemble of initial conditions, we refer to this quantity as K tot .

Results and Discussion
In this study, we present a statistical analysis of the work and the kinetic energy obtained during MD indentations performed in absence as well as in presence of thermal baths modeled by NH and Langevin thermostats. Attention is also given to the effect of the thermostat coupling on the left-hand side of the JE (Equation (11)), which is the measurable quantity in the Jarzynski theory. To this end, we consider two distinct indentation protocols that lead to different perturbations of the system. The load/unload indentation protocols are characterized in terms of maximum indenter penetration, h max , attained at t = τ/2. The dynamics of the particle system are defined according to the four distinct sampling methods described in Section 2.1. This gives us eight different indentation simulations. Each protocol is repeated over a large number of realizations (n = 1000). The individual realizations of the process corresponds to a different initial condition (at t = 0) drawn from the canonical distribution (2), produced during the NVT thermalization run during which the indenter potential Φ from Equation (24) is effectively 0; see Section 2.1.

Indentations with Elastic Deformations
The load/unload indentations with fixed maximum penetration h max = 0.1R + ∆ = 3.5 Å characterize the herein called elastic protocol; see Figure 1c. Our MD indentations with h max = 3.5 Å result in perturbations of the crystal that lead to elastic contacts between the indenter and the crystal's surface (Supplementary Figure S3), where the resulting P − h curves shown in Figure 2a follow a good agreement with the continuum elastic behaviour predicted by the Hertzian contact theory [21]. In this context, Figure 2a,b show that the unload stage approximately traces back the mechanical load path followed during loading, which suggests that the process is fully reversible. Additionally, notice that the load-unload curves are unaffected by the presence or absence of a thermostat. Interestingly, although the adiabatic and weakly thermostatted indentation are sampled from fundamentally different statistical ensembles, the corresponding P and K tot evolutions are practically identical; see the inset in Figure 2c. To confirm reversibility, Figure 2c further shows that the total kinetic energy of the particle system fluctuates around a constant K tot value irrespective of the thermostat coupling, even under unthermostatting conditions where there is strictly no coupling.  (see the table), where µ is the mean (in eV), σ 2 is the variance (in eV 2 ), and γ is the skewness of the K tot data.
Reversibility in the indenter-induced elastic perturbation is also evident in the work evolutions shown in Figure 3a,b. Along the time interval [0, τ], W S (t) gradually decreases from 0 to its minimum value at t = τ/2, and then, W S (t) increases during unloading, approximately matching the loading W S − h path; see Figure 3b. Note that the discontinuity in Figure 3a stems from the fact that the indenter's motion is inverted at t = τ/2 whereas the forces exerted by the indenter at t = (τ/2) − and at t = (τ/2) + remain identical. Figure 3c shows the W S histograms from n = 1000 realizations of the indentations with h max = 3.5 Å, which reveals fundamental details regarding the work fluctuations as a function of the thermostat coupling. We find that all W S distributions nearly adhere to normal distributions with values of skewness, γ, close to 0. In addition, the indentations with varying NH thermostat parameter ω p (namely with ω p = 100 dt and ω p = 100,000 dt) produce similar W S distributions with µ ≈ −0.44 eV and σ 2 ≈ 0.016 eV 2 . However, the W S data from the indentations with Langevin-thermostatted particles substantially differ from the other indentations, where the resulting W S distribution becomes considerably wider (with σ 2 = 0.116 eV 2 ) and shifts toward more negative values (with an average value of µ ≈ −2.4 eV). We attribute this to the damping induced by the Langevin thermostat, which hinders the (elastic) recovery to the initial state. This is evidenced in the W S − t/τ evolutions drawn in Figure 3b, which capture the gradual divergence of the W S values as t → τ with Langevin-thermostatted particles as compared with that with NH-thermostatted particles (with ω p = 100 dt). (Note that the difference between these W S (t = τ) values exists because it exceeds the statistical deviation of the W S data from the indentations with NHthermostatted particles. Additionally, note that a return to the initial state would require quasi-static transformations during which an efficient exchange of energy between the S and E occurs and W S (τ) = 0; see Ref. [10]. Clearly, these conditions are easier to obtain with Langevin baths as they tend to exchange energy with the system more efficiently than the NH thermostats that we employ).
In light of these results, we also find that the values of the measurable quantity in the JE, e −βW J , are affected by the thermostatting coupling; cf. the table in Figure 3c.
The plots of Figure 3d show the K tot histograms as a function of the thermostat coupling. We also find that these indentations produce normal K tot distributions with values of the skewness, γ, close to 0. The probability density function (PDF) describing such normally distributed K tot datasets can then be approximated by the general form of the Gaussian function, g(K tot , µ, σ) = (σ √ 2π) −1 exp[−(K tot − µ 2 )/2σ 2 ]; see the normalized histogram in Figure 3e. Contrarily to the W S distributions, the K tot histograms render similar values of the average (µ ≈ 1520 eV) and the variance (σ 2 ≈ 30-40 eV 2 ) regardless of the imposed thermostat coupling. Nonetheless, the indentation processes performed with a weak NH thermostat (ω p = 100,000 dt) and with unthermostatted particles statistically produce distributions marginally shifted toward larger K tot values as compared with those with thermostatted particles; cf. Figure 3d.

Indentations with Plastic Deformations
With deeper indenter penetrations than those produced by the elastic protocol, the mechanical response changes drastically. In contrast to the reversible elastic deformations discussed in the previous section, indentations with h max > 4.5 Å induce in the crystal non-reversible plastic deformations that persist over time. In general terms, crystal plasticity allows metals to sustain deformations beyond the elastic limit through the formation of non-reversible crystalline distortions. In the case of indented metallic crystals, plasticity manifests through the generation of crystalline defects under the indenter tip; see Supplementary Figure S4 and also Ref. [18]).
The indentations concerning a load/unload process with a maximum value of the indenter penetration of h max = 0.3R + ∆ = 9.5 Å characterize the herein called plastic protocol; see Figure 1c. In these indentations, the transition to a plastic stage is attained at penetrations of h ≈ h max /2. Hence, the imposed h max is well within the penetration range in which plastic features can be readily observed in the P − h curves.
The P-h curves from Figure 4a are characterized by an early elastic response followed by a marked load drop at the inception of plasticity. (Similar drops are observed in DNA pulling experiments [14]). With increasing penetration, further load drops result from the activation of additional plastic processes in the crystal. The P − h evolution then diverges from the elastic fit, thus manifesting the emergence of irreversibility in the system. Note that, unlike in the elastic protocol, the onset of plasticity leads to a marked increase in the instantaneous kinetic energy in the indentations with weakly NH-thermostatted (ω p = 100,000 dt) and unthermostatted particles; see Figure 4c. During unloading, the force vanishes at effective penetration values greater than 0 (or h f > ∆ in Figure 4a). In this regard, steep unloading curves are a fundamental manifestation of the formation of a plastic imprint during the indentation process, which remains in the particle system upon removal of the indenter tip from the crystal's surface. (For further details of such indentation responses in metals using nanometer-sized indenter tips, see Refs. [18,22]).
Irreversibility also becomes manifest in the W S time evolution obtained during the plastic protocol, where the load stage produces a greater absolute value of W S than during unloading (notice that W S takes a negative value over the [0, τ/2) time interval and a positive value over [τ/2, τ]); see Figure 5a,b. As a result, indentations with h max = 9.5 Å systematically lead to negative W S ; see Figure 5c. Interestingly, our analysis indicates that the W S histograms exhibit relatively similar moments irrespective of the thermostat coupling, with an expected W S value of ≈ −370 eV and a variance of σ 2 ≈ 3000 eV 2 . Somewhat unexpectedly, the W S distribution from the indentations with Langevin-thermostatted particles slightly differs from the deterministic approaches, as it only exhibits some left-hand skewness (γ ≈ 0.4). Our results from the plastic protocol indicate that the the left-hand side of the JE is uncomputable as e −βW J → 0; see the table in Figure 5c).
The resulting thermal fluctuations obtained in the plastic protocol are assessed through the K tot and K tot distributions shown in Figure 5d,e, respectively. Upon the occurrence of plastic deformations in the crystal, the K tot distributions from the indentations with strong system-bath coupling substantially differ in terms of the employed thermostat. Moreover, the indentations with a strong coupling show that the simulations using the Langevin thermostat lead not only to a wider normally distributed K tot histogram than that obtained with the NH thermostat (with σ 2 = 0.223 eV 2 and 0.045 eV 2 , respectively) but also to statistically larger values of K tot ; see Figure 5d. This essentially highlights the differences in performance of stochastic vs. deterministic thermostats. On the other hand, the indentations with the system weakly coupled to the thermal bath produce K tot bimodal distributions as a result of the thermal fluctuations obtained before and after the onset of plasticity in the crystal. In addition, when thermostatting is effectively inoperative in the plastic protocol, the thermal fluctuations during indentation become markedly wild, and thus, the K tot distributions shift toward larger values and become much wider, with a variance of σ 2 ≈ 120 eV 2 as compared with the indentations with thermostatted particles (σ 2 < 1 eV 2 ); see Figure 5e.   . (a,b): Evolution of dW S and W S in terms of t/τ. The plots concern a single realization using NH-thermostatted particles. (c): W S histograms from 1000 realizations. W S stands for W S (t = τ) obtained by means of Equation (26). (d,e): K tot and K tot histograms from 1000 realizations. The quantity K tot corresponds to the averaged value of K tot over τ obtained for each individual realization.

Conclusions
In this work, we study the effect that different deterministic and stochastic thermostats have on the mechanical and thermal properties of an indented Ta crystal. We perform MD simulations with unthermostatted and NH-and Langevin-thermostatted particles. The NH thermostats are characterized by the parameter ω p . With ω p = 100 dt, the NH bath effectively acts on the system while the NH thermostat with ω p = 100,000 dt allows only a limited exchange of energy between the system and the bath.
We present a systematic analysis of the work and kinetic fluctuations obtained in two distinct indentation protocols that produce reversible elastic and non-reversible plastic deformations in the crystal. Our main observations are summarized as follows: 1. In our indentation simulations, the system-bath coupling prescribed by the thermostats has a clear effect on the resulting work fluctuations. This is crucial when it comes to obtaining appropriate work statistics that enable free-energy difference calculations by means of the JE and related expressions. 2. In the MD indentations with unthermostatted and NH-thermostatted (with ω p = 100,000 dt) particles, the instantaneous kinetic energy of the system exhibits wild fluctuations when non-reversible plastic deformations are induced in the crystal. 3. The absence or presence of a stochastic thermostat in the dynamics of the particle system respectively represent the cases considered by the Hamiltonian and Langevin derivations of the JE. We find that the differences between the two approaches are substantial and bring about non-negligible effects in the calculation of the left-hand side of the JE. Such differences are clearly observable in the work distributions obtained under the fully reversible elastic protocol.