From proteins to grains: a journey through simple models.

i i “thesis” — 2018/9/7 — 8:45 — page i — #1 i i i i i i Departamento de F ´ ısica At´ omica Molecular y Nuclear Tesis Doctoral. Universidad de Sevilla From proteins to grains: a journey through simple models Doctorando: Carlos Alberto Plata Ramos Director: Antonio Prados Monta˜no Tutor: Diego G´omez Garc´ıa September 7, 2018 i i “thesis” — 2018/9/7 — 8:45 — page ii — #2 i i i i i i ii i i “thesis” — 2018/9/7 — 8:45 — page iii — #3 i i i i i i List of publications This thesis includes, at least partially, the research contained in the following works: •Carlos A. Plata, Fabio Cecconi, Mauro Chinappi, and Antonio Prados, Understanding the dependence on the pulling speed of the unfolding pathway of proteins, Journal of Statistical Mechanics P08003 (2015). •Alessandro Manacorda, Carlos A. Plata, Antonio Lasanta, Andrea Puglisi, and Antonio Prados, Lattice models for granular-like velocity fields: hydrodynamic description, Journal of Statistical Physics 164, 810 (2016). •Carlos A. Plata, Alessandro Manacorda, Antonio Lasanta, Andrea Puglisi, and Antonio Prados, Lattice models for granular-like velocity fields: finite-size effects, Journal of Statistical Mechanics 093203 (2016). •Carlos A. Plata and Antonio Prados, Global stability and theorem in lattice models with nonconservative interactions, Physical Review E 95, 052121 (2017). •Carlos A. Plata and Antonio Prados, Kovacs-like memory effect in athermal systems: linear response analysis, Entropy 19, 539 (2017). •Carlos A. Plata and Antonio Prados, Modelling the unfolding pathway of biomolecules: theoretical approach and experimental prospect. In Luis L. Bonilla, Efthimios Kaxiras, and Roderick Melnik (editors),Workshop on Coupled Mathematical Models for Physical and Biological Nanoscale Systems and Their Applications (Springer International Publishing, 2018). •Carlos A. Plata, Zackary N. Scholl, Piotr E. Marszalek, Antonio Prados, Relevance of the speed and direction of pulling in simple modular proteins, Journal of Chemical Theory and Computation 14, 2910 (2018). Other works that are not included in the thesis are: •Antonio Prados and Carlos A. Plata,Comment on “Critique and correction of the currently accepted solution of the infinite spherical well in quantum mechanics” by Huang Young-Sea and Thomann Hans-Rudolph, Europhysics Letters 116, 60011 (2017). iii i i “thesis” — 2018/9/7 — 8:45 — page iv — #4 i i i i i i iv i i “thesis” — 2018/9/7 — 8:45 — page v — #5 i i i i i i Table of contents 1 Introduction 1 1.1 Biophysics.................................... 2 1.1.1 Single-molecule experiments . . . . . . . . . . . . . . . . . . . . . . 4 1.1.2 Theoretical developments . . . . . . . . . . . . . . . . . . . . . . . 5 1.1.3 Unfolding pathway and its pulling dependence . . . . . . . . . . . 9 1.1.4 SummaryofpartI........................... 11 1.2 Granulargases ................................. 13 1.2.1 Granular hydrodynamics . . . . . . . . . . . . . . . . . . . . . . . . 15 1.2.2 Irreversibility: H-theorem....................... 17 1.2.3 Memory effects: Kovacs experiment . . . . . . . . . . . . . . . . . 21 1.2.4 Summary of part II . . . . . . . . . . . . . . . . . . . . . . . . . . 25 I Predicting the unfolding pathway of modular systems with toy models 27 2 The basics of modelling modular systems 29 2.1 The model fundamentals . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.2 Pulling the system: perturbative solution . . . . . . . . . . . . . . . . . . 32 2.2.1 Asymmetryterm............................ 36 2.2.2 Kineticterm .............................. 37 2.3 The critical velocities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 2.4 Comparing our theory with simulations . . . . . . . . . . . . . . . . . . . 40 2.4.1 Free energies of the units: shape and physical parameters . . . . . 40 2.4.2 Numericalresults............................ 42 2.5 Moving closer to the experiment . . . . . . . . . . . . . . . . . . . . . . . 50 2.5.1 Finitestiffness ............................. 50 2.5.2 Location of the elastic reaction . . . . . . . . . . . . . . . . . . . . 51 2.5.3 Units with different contour lengths . . . . . . . . . . . . . . . . . 52 2.6 Experimental prospect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 v i i “thesis” — 2018/9/7 — 8:45 — page vi — #6 i i i i i i vi Table of contents 3 Testing the model with molecular dynamics 57 3.1 Thecandidateprotein ............................. 58 3.2 All-atom molecular dynamics simulation . . . . . . . . . . . . . . . . . . . 58 3.2.1 Pulling trajectories . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.2.2 Unfolding criterion . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.3 Pulling the CC construct . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.3.1 C-pulling ................................ 63 3.3.2 N-pulling ................................ 65 3.4 Independence of the initial conformations in the pulling stage . . . . . . . 66 II Understanding granular matter with simple lattice models 71 4 Continuum limit, fluctuating hydrodynamics and finite size effects 73 4.1 Model and previous results . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.1.1 Evolution equations . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.1.2 Physical interpretation . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.1.3 Average fields in the continuum limit . . . . . . . . . . . . . . . . . 76 4.1.4 Fluctuating hydrodynamics . . . . . . . . . . . . . . . . . . . . . . 77 4.2 Dynamics of the one-particle distribution function . . . . . . . . . . . . . 78 4.3 The one-particle distribution function for some physical states . . . . . . . 80 4.3.1 The Homogeneous Cooling State . . . . . . . . . . . . . . . . . . . 81 4.3.2 The Uniform Shear Flow steady state . . . . . . . . . . . . . . . . 82 4.4 Study of the fluctuations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 4.4.1 Velocity fluctuations . . . . . . . . . . . . . . . . . . . . . . . . . . 86 4.4.2 Effect of velocity fluctuations on the total energy . . . . . . . . . . 87 4.5 Finite size effects: exact solution of the HCS . . . . . . . . . . . . . . . . 88 4.5.1 Eigenvalues for large systems . . . . . . . . . . . . . . . . . . . . . 92 4.5.2 Longtimelimit............................. 94 5 Global stability, derivation of the H-theorem 97 5.1 Basicsofthemodel............................... 98 5.1.1 The stochastic forcing . . . . . . . . . . . . . . . . . . . . . . . . . 98 5.1.2 Nonequilibrium steady states . . . . . . . . . . . . . . . . . . . . . 99 5.2 Proof of global stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 5.3 Why cannot HBbe the “good” Lyapunov functional? . . . . . . . . . . . 102 5.4 H-theorem for some specific NESS . . . . . . . . . . . . . . . . . . . . . . 104 5.4.1 USF state: simulations . . . . . . . . . . . . . . . . . . . . . . . . . 107 5.4.2 Numerical results in the uniformly heated system . . . . . . . . . . 109 6 Memory effects in athermal systems 113 6.1 Linear theory for the Kovacs memory effect . . . . . . . . . . . . . . . . . 113 6.1.1 Master equation approach . . . . . . . . . . . . . . . . . . . . . . . 113 6.1.2 Macroscopic equations approach . . . . . . . . . . . . . . . . . . . 117 6.2 Kovacs hump in the model of granular media . . . . . . . . . . . . . . . . 118 6.2.1 Kineticapproach............................ 118 6.2.2 Kovacs Hump in Linear Response . . . . . . . . . . . . . . . . . . . 121 i i “thesis” — 2018/9/7 — 8:45 — page vii — #7 i i i i i i Table of contents vii 6.2.3 Nonlinear Kovacs Hump . . . . . . . . . . . . . . . . . . . . . . . . 122 6.3 Numericalresults................................ 123 6.3.1 Validation of the first Sonine approximation . . . . . . . . . . . . . 123 6.3.2 KovacsHump.............................. 123 LinearResponse ............................ 123 NonlinearRegime ........................... 125 6.3.3 Time evolution of the H-Functional ................. 126 7 Conclusions 129 Bibliography 133 A Stability threshold 143 B Fluctuating expression for the dissipation 145 C Simulation method 149 D Lees-Edwards boundary conditions 151 E Derivation of the expression for dH/dt in a general driven state 153 F Completion of the proof of global stability 155 i i “thesis” — 2018/9/7 — 8:45 — page viii — #8 i i i i i i viii Table of contents i i “thesis” — 2018/9/7 — 8:45 — page 1 — #9 i i i i i i Chapter 1 Introduction How...? Why...? Curiosity is a natural instinct common to the whole mankind. Questioning is the very first step in any intellectual process, either science in general or physics in particular make no exception. Statistical mechanics was born to answer a question: how is the macroscopic world we see related to their microscopic components? Bernoulli, Maxwell, Boltzmann, Gibbs... all of them helped to establish the foundations of statistical mechanics. Nevertheless, the answer is not complete nowadays; fortunately for statistical physicists, there is still a lot of work to do. On the one hand, the scope is getting broader. Right now, we attempt to understand molecular biophysics, ecology, social sciences, and much more with the mathematical tools of statistical mechanics. On the other hand, the “right” theoretical framework for nonequilibrium statistical mechanics, in contrast to its equilibrium counterpart, is still under active development. This thesis is devoted to the analysis, through the lens of statistical mechanics, of two simple models motivated within two quite different fields: biophysics and kinetic theory of granular gases. The use of simple models for understanding, reproducing and predicting nature is a cornerstone in physics. The goal of this kind of modeling is to catch the essence of a complex system with the minimal, simplest, possible ingredients. The advantage of this approach is twofold. First, simplicity enables a (more) rigorous mathematical treatment, leaving the number of necessary approximations to a minimum. Second, the low number of ingredients allows us to isolate the features of a system that are responsible for the emergence of a certain behavior. Our work is divided into two parts, corresponding to the two aforementioned models. As stated above, we study them with the usual tools of statistical mechanics. That means that, depending on our level of description, our starting point is either Langevin-type equations, for the description of fluctuating physical quantities, or Fokker-Planck/master equations, for the description of probability density functions. On the one hand, in part I, we put forward and study a elasticity model for modular proteins capable of predicting the unfolding pathway of these macromolecules. On the other hand, we analyze a lattice model mimicking the main features of shear modes in granular gases in part II. 1 i i “thesis” — 2018/9/7 — 8:45 — page 8 — #16 i i i i i i 8 Chapter 1. Introduction Figure 1.4: Sketch of the WLC model. The protein chain tends to align with the pulling force F. strong pulling the extension of the system approaches its contour length Lcand the force diverges as either (L−Lc)−1for the FJC model, or (L−Lc)−2for the WLC model. Both the FJC and WLC models have been used for fitting real experiments with biomolecules, giving reasonably good results [29–31]. Specifically, the following WLC fit [29] FP kBT=L Lc +1 4Lc Lc−L2 −1 4,(1.6) is asymptotically valid along all the length range, and it is intensively employed in the literature. However, the aforementioned paradigmatic models do not take into account the internal structure of the chain. In fact, the internal structure of the biomolecule is responsible for its different states, folded or unfolded for instance. A coarse-grained modeling usually involves considering each unit within a macromolecule as a two-state system, which can be in either a folded or an unfolded state. To account for this, some models [32–34] consider a WLC with several possible values of the contour lengths: each branch of the force extension curve is fitted by a WLC model with a different value of contour length. Transitions between folded and unfolded states typically follow the development of Kramers theory [35, 36] carried out by Bell [37] and Evans [38]. The Bell-Evans expression provides the transition rates between states, given the applied force and the details of the free energy barrier. Quite recently, some more theoretical models [39–41], closer to the approach that will be followed in this thesis, have been proposed to analyze the elasticity of modular biomolecules. These models successfully explain the sawtooth pattern observed in the experiments: interestingly, an equilibrium-statistical-mechanics theory is sufficient to understand their emergence. In a nutshell, the models start from a free energy where each module gives an additive contribution thereto, the individual contributions being double well functions of the corresponding unit’s extension. By maximizing the probability dis- i i “thesis” — 2018/9/7 — 8:45 — page 9 — #17 i i i i i i Section 1.1. Biophysics 9 tribution function within the right statistical ensemble (force-control or length-control), or equivalently minimizing the corresponding free energy, the equilibrium force-extension curve of the system is obtained. As a consequence of the different metastable configurations (folded or unfolded) for each module, the force-extension curve presents different branches, see top left panel of figure 1.5. When dynamics is incorporated, hysteresis processes might appear. In figure 1.5, some force-extension curves of a system of 8 units for different velocities are shown. The hysteretic behavior strongly depends on the pulling rate employed. More specifically, an interesting interplay between the pulling velocity and the temperature is observed. At a certain value of the temperature, the system mainly sweeps the equilibrium forceextension curve if the pulling speed is slow enough: this is the quasistatic regime of pulling (bottom left). However, for higher velocities, hysteresis comes about (top right). This phenomenon is accentuated as the temperature is lowered, since thermally activated transitions are not possible in “cold” systems. This allows the system to sweep the whole metastable branches (bottom right), a regime that is usually said to correspond to the “maximum hysteresis path” [42] (also adiabatic pulling [41]). In fact, the maximum hysteresis path regime described above stems from the interplay between the pulling velocity and the temperature: the pulling velocity must be slow enough to allow the system to move over the (metastable) equilibrium branches but fast enough to prevent the system from going from the folded to the unfolded basin (or vice versa) by thermal activation. This regime will be of crucial relevance in the development of the theory presented in part I of this thesis. 1.1.3 Unfolding pathway and its pulling dependence The unfolding pathway is, roughly, the order and the way in which the structural blocks of a macromolecule unfold. The force-extension curve obtained in single-molecule experiment characterizes the elasto-mechanical behavior of the macromolecule and provides basic and essential information about the unfolding pathway. Some studies show that the pulling velocity plays a relevant role in determining the unfolding pathway, see for example [12,17,43,44]. Different unfolding pathways have been observed depending on (i) pulling direction, that is, which of the ends (C-terminus or N-terminus) the molecule is pulled from and (ii) the pulling speed [12,17,43,44]. Intuitively, it has been claimed that it is the inhomogeneity in the distribution of the force across the protein, for high enough pulling speeds, that causes the unfolding pathway to change [17]. Nevertheless, to the best of our knowledge, a theory that explains this crossover is still lacking. Our interest in this problem was triggered by the dissenting unfolding pathway observed in the maltose binding protein (MBP) in experiments [16] and simulations for high pulling velocity [17]. This molecule unfolds in four steps, which led Bertz and Rief to identify four internal substructures they call “unfoldons”, see figure 1.6. More specifically, AFM experiments allowed them to characterize these four unfoldons in the MBP, labeling them as M1, M2, M3 and M4, with M1 being the closest to the C-terminus and M3 the closest to the N-terminus [16]. Pulling the molecule at a typical speed of 10−9nm/ps, Bertz and Rief found a well-defined unfolding pathway: the weakest unfoldon (M1), that is, the one characterized by the lowest opening force, unraveled first. Thereafter, the i i “thesis” — 2018/9/7 — 8:45 — page 10 — #18 i i i i i i 10 Chapter 1. Introduction Figure 1.5: (Top left) Equilibrium force rips in the force-extension curve for a system with 8 units (lengthcontrol). Different colors are used to represent the stability: solid black for stable parts, dotted red for metastable parts and black arrows for force rips. In a quasistatic pulling process, the system follows the solid black curve with a series of first-order transitions in the force. (Top right) Hysteresis cycle for a system composed of 8 units, at a relatively high pulling speed and moderate temperature. Two traces are plotted: solid blue for unfolding and dashed red for refolding. (Bottom left) Same plot as in the top right panel, but for a smaller pulling velocity. Aside from thermal fluctuations, the system almost sweeps the equilibrium curve: the pulling is basically quasistatic (Bottom right) The same plot as in the bottom left panel, but for an almost vanishing temperature. Thermal fluctuations are so small that the system approaches the T= 0 behavior (adiabatic pulling): therein, the branches are swept up to the end of the metastability region. Taken from [41]. i i “thesis” — 2018/9/7 — 8:45 — page 11 — #19 i i i i i i Section 1.1. Biophysics 11 Figure 1.6: Structure of the Maltose Binding Protein (PDB ID: 4MBP). The color code identifies the four unfoldons: M1 blue; M2 green; M3 red; M4 gold. (Left) Crystallographic structure. (Right) Topological diagram. Taken from [17]. remainder of the unfoldons opened sequentially, from weakest to strongest (M4). Later, Guardiani et al. had a more detailed look into the unfolding of the MBP [17]. They showed, by means of a combination of Go model simulations and steered molecular dynamics, that the unfolding pathway is more complex and seems to depend on both the velocity and direction of pulling. C-pulling simulations of the Go model always showed a pathway compatible with Bertz and Rief’s experiment, see top panel of figure 1.7. However, N-pulling simulations of the Go model displayed a different behavior: for small velocities, again Bertz and Riefs pathway was found, but for high enough pulling speed it was M3, the closest to the N-terminus, that opened first, as shown in the bottom panel of figure 1.7. Steered molecular dynamics simulations at a pulling speed of 5 ·10−3nm/ps gave results that were consistent with those from the Go model at high velocities, showing the different pathways (M1 vs. M3 opening first) depending on the pulled terminus (C vs. N). In this context, Guardiani et al. introduced in [17] a toy model to qualitatively explain the observed pathways, which is the starting point of part I of this thesis. This simple model is akin to those employed in [39–41] to investigate the force-extension curves of modular proteins, their main difference stemming by the incorporation—in the simplest way—of the spatial structure of the chain into Guardiani et al.’s model. Numerical simulations of the latter presented a phenomenology that was compatible with both the Go model and the steered molecular dynamics results [17]. One of our aims is to develop a theoretical framework for this simple model, in order to get a deeper understanding of the dependence of the unfolding pathway on the pulling velocity and direction. 1.1.4 Summary of part I Despite the large number of models developed to unravel the nature of biomolecules, not everything is neat. In part I, we attempt both to better understand and to predict the unfolding pathway of modular biomolecules. Our approach to this problem follows the philosophy presented above, we put forward i i “thesis” — 2018/9/7 — 8:45 — page 12 — #20 i i i i i i 12 Chapter 1. Introduction Figure 1.7: Evolution of the gyration radius for unfoldons M1, M2, M3, and M4 as functions of the end- to-end distance (˚ A) at different pulling speeds. The curves show averages over 50 runs of simulations of the Go model, for C-pulling (above) and N-pulling (bottom). For C-pulling, it is always M1, the weakest unfoldon, that opens first (solid blue line). Notwithstanding, for N-pulling, the first unfoldon that unravels depends on the pulling velocity: for slow pulling (v= 0.001˚ Aτ−1) it is still M1, but for fast enough pulling (v= 0.1˚ Aτ−1) it is M3, the closest to the pulled end as depicted in figure 1.6 (dot-dashed red line). Taken from [17]. i i “thesis” — 2018/9/7 — 8:45 — page 13 — #21 i i i i i i Section 1.2. Granular gases 13 a simple model with the minimal ingredients to grasp the essence of the system. We do so in chapter 2. Therein, we introduce the model, which portray the conformation of the protein into a 1d chain with different units. Each unit contribute to the global free energy with a function that only depends on its own extension with a double-well shape. The perturbative solution of the dynamical system leads to our predicting of the first module that unfolds. Numerical results are in excellent agreement with our theoretical predictions. Additionally, we carry out some modifications of the model, in order to get closer to real experiments. The results are unchanged at the lowest (leading) order, thus proving the robustness of our analysis. Finally, we test the theoretical scheme developed in chapter 2 in a simple biomolecule, which comprises two coiled-coil structures, in chapter 3. This analysis is carried out by means of steered molecular dynamics simulations of the coiled-coil construct. First, characterizing the molecule is required: in particular, we need to introduce a criterion for considering it unfolded or folded. Such a criterion allows us to make a systematic comparison with the theoretical framework. A thorough statistical study of the simulations output provides a significant test of theory and validates the usefulness of the approach. Some technical details that are omitted in the main body of this part of the thesis are given in appendix A. 1.2 Granular gases A granular material is made of macroscopic particles that are called grains [45,46]. They can be found almost everywhere, for example dust, sand, seeds, pills, iceberg groups or asteroid populations are all instances of granular matter. Improving our knowledge of granular matters has a clear technical, industrial, and even economic interest. To support this statement, it suffices to turn our thoughts to transport and storage industry, agriculture or construction. In figure 1.8, some typical examples of granular material are shown. All granular materials share some typical properties. First, grains are solid and macroscopic, that is, their dynamics is governed by classical mechanics laws and they fill a space that is excluded to the other grains. Second, their interactions are nonconservative in the following sense: when two grains collide, some of their energy is “lost” into internal degrees of freedom, mainly due to deformation, as heat. Third, the characteristic energy of a grain is much greater than the thermal energy. Therefore, the temperature of the medium in which the grains are immersed is largely irrelevant: the typical energy to lift a grain by its own diameter is several orders of magnitude larger than the thermal energy [48]. Indeed, in granular fluids one can define a relevant granular temperature from kinetic theory, linked to velocity fluctuations. This granular temperature has nothing to do with the “conventional” temperature, which as already stated plays no role, but is related to the energy injection mechanism that is needed to keep grains moving. The result of enclosing a granular material and shaking it rapidly is a fluidized granular material: the granular fluid that we have just referred to above. The amount of available space and the intensity of the shaking determine the regime of fluidization [49]. When interactions are dominated by two-particle instantaneous (hard-core-like) collisions, one usually speaks of a granular gas. This regime is typically achieved in experiments when i i “thesis” — 2018/9/7 — 8:45 — page 14 — #22 i i i i i i 14 Chapter 1. Introduction Figure 1.8: Examples of granular matter: balls, gravel, lentils and sesame seeds. Taken from [47]. the packing fraction is of the order of ∼1% or less, and the peak acceleration is many times the gravity acceleration. The gas regime has played a crucial role in the development of granular kinetic theory: in the dilute limit, one may retrace the classical molecular kinetic theory after having relaxed the constraint of energy conservation [50]. Granular collisions are, in fact, inelastic: this occurs because each grain is approximated as a rigid body and the collisional internal dynamics is replaced by an effective energy loss, usually characterized by a normal restitution coefficient α∈(0,1]. This is the smooth hard-particle model: particles only interact when at contact, then the component of the relative velocity along the direction joining their centers is reversed and shrunk by a factor α, whereas the other components remain unchanged. Therefore, α= 1 corresponds to elastic collisions, whereas α→0 describes the completely inelastic limit. Most of granular kinetic theory rests upon many variants of the basic model of smooth hard particles. Important variants include roughness and rotation [51–55] as well as the consideration of velocity-dependent inelasticity [56]. Notwithstanding, the simple model of smooth inelastic hard spheres suffices to explain the basic phenomenology of granular gases. In this context, the Boltzmann equation for inelastic hard spheres constitutes the foundation of many investigations in the realm of granular phenomena, with both numerical and analytical approaches [57]. Granular fluids exhibiting separation between fast microscopic scales and slow macroscopic ones have led to several procedures to build a granular hydrodynamics [58, 59]. Note that the scale separation hypothesis is less clear in the granular case, as compared to molecular gases. First, one has the spontaneous tendency of granular gases to develop i i “thesis” — 2018/9/7 — 8:45 — page 15 — #23 i i i i i i Section 1.2. Granular gases 15 strong inhomogeneities even at the scale of a few mean free paths. Second, a granular system is typically of “small” size: it is usually constituted by a few thousands of grains [60,61]. The latter limitation cannot be easily relaxed even in theoretical studies: stability analyses have shown that spatially homogeneous states are unstable for too large sizes in some typical states [62]. The intrinsic “small” size of granular gases makes it essential to address another task: an adequate and consistent description of fluctuations, which are always important in a small system [63–65]. The number of particles N, ranging from 102to 104, is large enough to make it possible to apply the methods of statistical mechanics, but definitely much smaller than the Avogadro number. Interestingly, this is also the case for biomolecules, as stated before in section 1.1.2; the special relevance of fluctuations is a point that links the two parts of this thesis. Unfortunately, there is no general theory currently available for mesoscopic fluctuations out-of-equilibrium. Notwithstanding, important steps in the deduction of a consistent fluctuating hydrodynamics for inelastic hard spheres have been recently taken in the context of kinetic theory [66]. In this regard, the quite broad framework of Macroscopic Fluctuation Theory [67] cannot be employed. In its current state of development, Macroscopic Fluctuation Theory does not include macroscopic equations with advection terms and momentum conservation, such as those in the “granular” Navier-Stokes equations. In the last decades, lattice models have proved to be a flexible tool to identify the essential steps in a rigorous approach to the hydrodynamic limit, both at the average [68,69] and fluctuating [67] levels of description. Fluctuating hydrodynamics in linear and nonlinear lattice diffusive models have been recently investigated, both in the conservative [70–73] and in the nonconservative cases for the energy field [74–78]. Later, a lattice model, which in some simplified way mimics the velocity field of a granular gas, has been put forward to incorporate momentum conservation [79]. This model is the central pillar of our investigations in part II of this thesis. 1.2.1 Granular hydrodynamics Evolution equations for the “slow” fields, in space rand time t, density n(r, t), velocity u(r, t), and granular temperature T(r, t) constitute the full granular hydrodynamics. These equations can be derived from the Boltzmann equation for inelastic hard spheres, through a Chapman-Enskog procedure closed at the Navier-Stokes order [46, 59]. For generic dimension d, they are given by ∂tn+∇·(nu)=0,(1.7a) ∂tu+u·∇u+ (nmkB)−1∇·P= 0,(1.7b) ∂tT+u·∇T+2 dnkB [P: (∇u) + ∇·q] + ζT = 0,(1.7c) in which mis the mass of the particles. The energy dissipation rate is ζ(r, t)>0, being ζ→0 in the elastic limit, while the pressure tensor P(r, t) and the heat flow q(r, t) read Pij =pδij −η∇iuj+∇jui−2 dδij∇·u,q=−κ∇T−µ∇n, (1.8) i i “thesis” — 2018/9/7 — 8:45 — page 16 — #24 i i i i i i 16 Chapter 1. Introduction Additionally, the bulk pressure p, the shear viscossity η, the heat-temperature conductivity κand the heat-density conductivity µare given by constitutive relations [46]. We briefly recall that the transport coefficients depend on the hydrodynamic fields. Specifically, for hard-spheres one has that η∼√T,κ∼√T,µ∼n−1T3/2and ζ∼n√T[46]. This thesis is not an exhaustive investigation of granular hydrodynamics and its rich catalogue of possible stationary and nonstationary regimes [49,50]. Our aim is to study the model introduced in part II, validating its use for the simplified investigation of some peculiar states in the granular realm. With this intention, we highlight three essential aspects that give some contact points of the aforementioned model with actual granular fluids. First, the existence of a spatially homogeneous nonstationary solution, that is, the “Homogeneous Cooling State” (HCS). Second, the instability of such a state with respect to perturbations with long enough wavelength (small enough wavenumber). And third, the existence of the Uniform Shear Flow (USF) stationary state, in which the energy loss due to collisions is balanced—on average—by the heating brought about by the velocity difference, that is, the shear, imposed between the system boundaries [80]. All such aspects stem from a key difference with respect to the hydrodynamics of molecular fluids, which is the presence of the energy sink term ζT in (1.7). When spatial homogeneity is assumed along with periodic boundary conditions, and initial conditions n(r, t = 0) = n,u(r, t = 0) = 0, and T(r, t = 0) = T(0), (1.7) are reduced to ˙ T(t) = −ζ(t)T(t).(1.9) Since ζ(t)∝T(t)1/2for hard-spheres, (1.9) leads to the well known Haff’s law [81] THCS(t) = T(0) 1 + ζ(0)t 2−2 .(1.10) A different collisional model that is often used to simplify the kinetic theory approach is the so-called gas of pseudo-Maxwell molecules [82]: its peculiarity is that ζ(t) = ζ(0) is constant and therefore Haff’s law simplifies to an exponential decay, THCS(t) = T(0) exp[−ζ(0)t].(1.11) This spatially homogeneous solution, with monotonically decreasing temperature, is generally called “Homogeneous Cooling State” (HCS). Of course, it can be predicted at the more fundamental and general level of the Boltzmann [83] or even Liouville [84] equations. The HCS is not stable if the system size exceeds some critical value. Spatial perturbations of the velocity and density fields are amplified when the system is large enough [62, 85]. A linear stability analysis shows that the fastest amplification occurs for shear modes, which correspond to a transverse perturbation of the velocity field. For instance, a nonzero ycomponent of umodulated along the xdirection, that is, uy(x, t). The critical wavelength Lcseparating the stable from the unstable regime depends upon the restitution coefficient as L2 c∝(1 −α2)−1.(1.12) Velocity perturbations are not really amplified, because the amplitude of their fluctuations (temperature) always decay: the instability is observed only when the rescaled velocity i i “thesis” — 2018/9/7 — 8:45 — page 17 — #25 i i i i i i Section 1.2. Granular gases 17 field u(r, t)/pTHCS(t) is considered. Perturbations in the other fields (density, longitudinal velocity and temperature), the evolution of which are coupled with the velocity field, are also amplified, but with a slower rate and for longer wavelengths. There is a range of system sizes such that the only linearly unstable mode is the shear mode [50]. This entails that the velocity field is incompressible and density does not evolve from its initial uniform configuration. Such a regime may be observed for a certain amount of time, longer and longer as the elastic limit is approached. In two dimensions, (1.7) is obeyed with constant density and, for instance, ux= 0 whereas the hydrodynamic fields uyand Tonly depend on x. In this situation, we have that (1.7) reduces to ∂tuy(x, t)=(nm)−1∂x[η∂xuy(x, t)],(1.13a) ∂tT(x, t) = 1 nkB η[∂xuy(x, t)]2+1 n∂x[κ∂xT(x, t)] −ζT. (1.13b) In section 4.1.3 we will see that our lattice model is well described, in the continuum limit, by completely analogous equations. It is interesting to put in evidence a particular stationary solution of the system (1.13). Seeking time-independent solutions thereof, one finds ∂x[η∂xu(s) y(x)] = 0,η kB [∂xu(s) y(x)]2=−∂x[κ∂xT(s)(x)] + nζT(s)(x).(1.14) The general situation is that both the average velocity and temperature profiles are inhomogeneous: this is the so-called Couette flow state, which also exists in molecular fluids. Nevertheless, in granular fluids, there appears a new steady state in which the temperature is homogeneous throughout the system, T(s)(x) = T, and the average velocity has a constant gradient, ∂xu=a: this is the Uniform Shear Flow (USF) state, characterised by the equations ∂2 xu(s) y(x) = 0, η[∂xu(s) y(x)]2=nkBζT (s).(1.15) Such a steady state is peculiar of granular gases where the viscous heating term is locally compensated by the energy sink term. In granular gases of hard spheres, it has been proven that the USF state is linearly stable for perturbations in the direction of the shear [86]. 1.2.2 Irreversibility: H-theorem In thermodynamics and statistical mechanics, proving the global stability of the equilibrium state—or, in general, of the relevant stationary state—usually involves the introduction of a suitable Lyapunov functional [87]. A Lyapunov functional of the probability distribution function (PDF) has the following three properties: (i) It is bounded from below. (ii) It monotonically decreases with time. (iii) Its time derivative vanishes only when the PDF is the equilibrium one. i i “thesis” — 2018/9/7 — 8:45 — page 24 — #32 i i i i i i 24 Chapter 1. Introduction variables should be taken into account, Twill not remain constant but separate from its steady (initial) value and have either a maximum or a minimum. In molecular systems, there always appears a maximum in the Kovacs hump. This does not have to be the case for the granular temperature, because the granular gas is an intrinsically dissipative, out-of-equilibrium system. Defining the shifted time variable τ¼ζ0ffiffiffiffi Ts pðt−twÞ=2, we have to solve Eqs. (6) with the initial conditions βðτ¼0Þ¼1and A2ðτ¼0Þ¼aini 2=as 2, where aini 2is the value of the excess kurtosis in the final state of the waiting time window. Since as 2is small (jas 2j≤0.07) across the whole range of restitution coefficients, while βand A2are of the order of unity, we expand both βand A2in powers of as 2to obtain an approximate solution of Eqs. (6) [46], a2ðτÞ−as 2∼ðaini 2−as 2Þe−4Bτ;(9a) βðτÞ−1∼3ðaini 2−as 2Þ 16ð4B−3Þðe−3τ−e−4BτÞ:(9b) The relaxation of the excess kurtosis to its steady value is exponential, while that of the rescaled temperature βis the sum of two exponentials with different relaxation times. The sign of β−1is the same as that of as 2because (i) 4B>3and (ii) ðaini 2−as 2Þand as 2have the same sign as a function of the restitution coefficient for the arbitrary “cooling”(ξ0>ξ>ξ1) protocol in Fig. 1. In fact, Eq. (6b) predicts that dA2=dτis initially positive and, thus, ja2j>jas 2jin the whole waiting time window [46]. In addition, the steady excess kurtosis as 2changes sign at αc¼1=ffiffiffi 2 p≃0.707:as 2>0for α<αcwhile as 2<0for α>αc[47]. Thus, for small inelasticity (α>αc), β−1<0, and βhas a minimum, while the granular temperature T¼Ts=β2has a maximum. This behavior is completely similar to that of glassy systems, so we may speak of a normal Kovacs hump in the weakly dissipative case. On the contrary, for high inelasticity, α<αc, β−1>0, and βdisplays a maximum, which corresponds to a minimum of T; an anomalous Kovacs hump appears. In Fig. 2,theabovetheoreticalpredictionfortheKovacs hump is tested against numerical computations. The latter are obtained by means of direct Monte Carlo simulations [48] of the Boltzmann-Fokker-Planck equation (1).Twovaluesof FIG. 1 (color online). Top: Sketch of the drive time dependence. Bottom: Ensuing temperature evolution. At t¼0, the gas is at temperature Tsðξ0Þ, in the nonequilibrium steady state corresponding to a value of the driving ξ0. At t¼0, the drive is suddenly decreased to ξ1≪ξ0, which is kept for a waiting time tw. At t¼tw, the granular temperature is measured, and the driving is cranked up to a new value ξ, such that TsðξÞ¼TðtwÞ. The question mark is for the two possible scenarios: a positive hump with a Tmaximum (normal behavior, solid line) or a negative anomalous hump (dot-dashed line). At long times, T reaches its steady value TsðξÞ. 0.0 0.5 1.0 1.5 2.0 0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 1 0.0 0.5 1.0 1.5 2.0 0.00015 0.00010 0.00005 0.00000 1 FIG. 2 (color online). Plot of the Kovacs hump for α¼0.3(top) and α¼0.8(bottom). Monte Carlo simulation curves (points) for a system of 104hard disks (d¼2) averaged over 105(top) and 1.5×106trajectories (bottom). They are compared to the theoretical curve (9b): The dashed line corresponds to the predicted values of as 2,aini 2, and B, while the solid line is obtained by taking these three parameters from the simulation (see, e.g., Fig. 3, from which Bis directly measured). The sign of β−1changes from the highly inelastic (top) to the weakly inelastic (bottom) case. Note that a maximum of βcorresponds to a minimum of T¼Ts=β2(and vice versa) so that the Kovacs hump is anomalous in the highly inelastic region. PRL 112, 198001 (2014) PHYSICAL REVIEW LETTERS week ending 16 MAY 2014 198001-3 Figure 1.11: Emergence of anomalous Kovacs response in a uniformly heated granular gas. Specifically, β=pTs/T, where Tsis the steady value of the temperature, and τmeasures time in the number of collisions per particle. Therefore, a maximum in the temperature corresponds to a minimum in βand vice versa. (Top) Highly inelastic case (restitution coefficient α= 0.3), for which the anomalous Kovacs response is clearly observed. (Bottom) Inelasticity is decreased (α= 0.8, closer to the elastic case α= 1) and the Kovacs response becomes normal. The crossover from normal to anomalous response takes place at α= 1/√2; a detailed discussion can be found in [110], from which the figure is taken. i i “thesis” — 2018/9/7 — 8:45 — page 25 — #33 i i i i i i Section 1.2. Granular gases 25 1.2.4 Summary of part II As advanced above, part II is devoted to the thorough analysis of a lattice model that attempts to catch the essential phenomenology of the shear modes of granular fluids. In each chapter, we explore different aspects of the model. We start by introducing the model in chapter 4. Therein, we focus on the continuous hydrodynamic-like limit that, to the lowest order, leads to equations completely analogous to those in (1.13). In that limit, we study different relevant physical states. Specifically, we analyze the Homogeneous Cooling State and the Uniform Shear Flow state through its one-particle velocity distribution. We also go beyond the aforementioned lowest order analysis in two ways: (i) looking into the behavior of the fluctuating fields and (ii) solving exactly the Homogeneous Cooling State on the lattice. In chapter 5, we turn our attention to the stability of the NESS of this model system. Not only have we proven the global stability of a quite general family of states, but also clarified the inadequacy of HB, given by (1.16), as a Lyapunov functional. We finish the chapter with a rigorous proof of an H-theorem. To the best of our knowledge, our result constitutes the first proof, even in simple models, of an H-theorem within the context of systems with nonconservative interactions. Finally, chapter 6 is devoted to the analysis of Kovacs-like memory effects. We develop a general theoretical framework for the linear response analysis in athermal systems, starting from either the master equation for the probability distribution function (mesoscopic description) or the evolution equations for the macroscopic moments (macroscopic description). Our results are particularized for a variant of our lattice model of granular gas, and they show an excellent agreement with simulations. Although we test the theory in our specific model, it is worth noting the quite broad range of physical systems that our developed theoretical framework can be directly applied to. Appendices B-F deal with some technicalities that are skipped within the main text of this part of the thesis. i i “thesis” — 2018/9/7 — 8:45 — page 26 — #34 i i i i i i 26 Chapter 1. Introduction i i “thesis” — 2018/9/7 — 8:45 — page 27 — #35 i i i i i i Part I Predicting the unfolding pathway of modular systems with toy models 27 i i “thesis” — 2018/9/7 — 8:45 — page 28 — #36 i i i i i i i i “thesis” — 2018/9/7 — 8:45 — page 29 — #37 i i i i i i Chapter 2 The basics of modelling modular systems Over the course of this first part of the thesis, we study in depth a simple model of elasticity for biomolecules. This model, introduced by Guardiani et al. [17] for simulation purposes, lacked a thorough theoretical analysis. Herein, we specifically look for a theory capable of predicting the unfolding pathway in modular biomolecule submitted to mechanical pulling. We expect that this theory could also explain the unfolding pathway observed in experiments [16] and simulations [17] of the maltose binding protein. As described in the introduction chapter, the unfolding pathway depends on the pulling velocity: at very low pulling rates, it is the weakest unit that unfolds first, while at higher rates the first unit to unravel is the pulled one. This chapter is dedicated to the development of the aforementioned theory and its plan is detailed below. We start introducing the basics of the model in section 2.1. The system dynamical response to mechanical pulling is obtained by means of a perturbative approach in section 2.2. Section 2.3 is devoted to the obtention of the solution of the dynamical equations, which allows us to study the emergence of a set of critical velocities at which the unfolding pathway changes. Section 2.4 is dedicated to the comparison of the theoretical predictions with numerical results of the model equations. We seek more realistic variants of the model in section 2.5. Finally, we propose a possible experiment for testing our theory in section 2.6. 2.1 The model fundamentals Let us consider a certain modular biomolecule comprising Nmodules. The paradigmatic example is a polyprotein composed of N, possibly different, modules. Notwithstanding, we may also be considering a protein domain with Nunfoldons. From now on, we will refer to these modules or unfoldons, indistinctly, as modules or units. When the system is submitted to an external force F, the simplest description is to portray it as a onedimensional chain in the direction of the force. We denote the end-to-end extension of 29 i i “thesis” — 2018/9/7 — 8:45 — page 30 — #38 i i i i i i 30 Chapter 2. The basics of modelling modular systems the i-th unit by xi. In a real AFM experiment, the molecule is attached as a whole to the AFM device and stretched. Following Guardiani et al. [17], we model this system with a sequence of nonlinear bonds, as in figure 2.1. Therein, the endpoints of the i-th unit are denoted by qi−1and qi, so that its extension xiis xi=qi−qi−1, i = 1, . . . , N. (2.1) In this basic model, we consider that the left end of the first unit is fixed, that is, q0= 0 for all times. We assume that the inertia terms can be neglected and the evolution of the system follows the coupled overdamped Langevin equations γ˙qi=−∂ ∂qi A(q0, . . . , qN) + ηi,(2.2) in which ηiare Gaussian white noise terms. They verify hηi(t)i= 0,hηi(t)ηj(t0)i= 2γkBTδijδ(t−t0),(2.3) with γbeing the friction coefficient of each unit (the same for all), Tthe temperature of the fluid in which the protein is immersed, and kBthe Boltzmann constant. The global free energy function of the system is A(q0, . . . , qN) = N X i=1 ai(qi−qi−1) + ap(qN).(2.4) In (2.4), ap(qN) is the contribution to the free energy introduced by the force control or length control device, see below, while ai(xi) is the contribution to Astemming from the i-th unit, which is only function of its own extension xi. The total length of the system is given by N X i=1 xi=qN.(2.5) In force control experiments, the applied force Fis a given function of time, whereas in length control experiments the device (portrayed by the spring in figure 2.1) tries to keep the total length qNequal to the desired value L, also a certain function of time. The corresponding contributions to the free energy are ap(qN) = −FqN,force control,(2.6a) ap(qN) = 1 2kc(qN−L)2,length control,(2.6b) in which kcstands for the stiffness of the length control device. The length is perfectly controlled in the limit kc→ ∞, when qN=Lfor all times. For the sake of a common notation, we have not used different letters for Helmhotz or Gibbs free energies. It has to be understood that, on the one hand, introducing (2.6a) in (2.4), we obtain a Gibbs free i i “thesis” — 2018/9/7 — 8:45 — page 31 — #39 i i i i i i Section 2.1. The model fundamentals 31 Figure 2.1: Sketch of the basic model for a protein with four units. Each unit is depicted by a rectangle with nonharmonic free energy ai(xi). The beads mark the coordinates qiof their endpoints, so that the i-th unit’s extension xi=qi−qi−1(by definition, q0= 0). Finally, the spring represents the device attached to the pulled end q4, which controls either the applied force (force control) or the end-to-end distance q4(length control). The device contribution to the free energy is ap(q4), see (2.4) and (2.6). energy. On the other hand, if we use (2.6b) in (2.4), and the limit kc→ ∞ (taking into account that ap→0 since the force kc(L−qn) goes to a constant Lagrange multiplier), the result is a Helmholzt free energy whose minimum, restrained to the total length constraint, gives the equilibrium configuration of the system. An apparently similar system, briefly discussed in section 1.1.2, in which each module of the chain follows the Langevin equation γ˙xi=−∂A/∂xi+ηi, has been analyzed in the literature [40,41]. In this approach, the modules are completely independent in force control experiments, the global free energy is the sum of individual ones and, because of this, Langevin equations completely neglect the spatial structure of the chain. While this simplifying assumption poses no problem for the characterization of the force-extension curves in [41], it is not suited for the investigation of the unfolding pathway. In this context, the spatial structure plays an essential role. The spatial structure of biomolecules can be described in quite a realistic way by using the framework proposed by Hummer and Szabo several years ago [115], but our simplified picture in figure 2.1 makes an analytical approach feasible. Now, we look into the unfolding pathway of this system. As the evolution equations are stochastic, this pathway may vary from one trajectory of the dynamics to another. Nevertheless, in many experiments [16,17, 43] a quite well-defined pathway is observed, which suggests that thermal fluctuations do not play an important part in its determination. Physically, this means that the free energy barrier separating the unfolded and folded conformations at coexistence—that is, at the critical force, see below—is expected to be much larger than the typical energy kBTfor thermal fluctuations. Therefore, we expect the thermal noise terms in our Langevin equations to be negligible and, consequently, they will be dropped in the remainder of our theoretical approach. Of course, if the unfolding barrier for a given biomolecule were only a few kBTs, the thermal noise terms in the Langevin equations could not be neglected and our theoretical approach would have to be changed. In order to undertake a theoretical analysis of the stretching dynamics, we introduce one further simplification of the problem. We consider that the device controlling the length is perfectly stiff, thus the total length qN=Ldoes not fluctuate. We expect this assumption to have little impact on the unfolding pathway: otherwise, the latter would be more a property of the length control device than of the chain. In fact, we show in section i i “thesis” — 2018/9/7 — 8:45 — page 32 — #40 i i i i i i 32 Chapter 2. The basics of modelling modular systems 2.4.2 that the unfolding order is not affected by this simplification. For perfect length control, the mathematical problem is identical to that of the force control situation, but now the force Fis an unknown (Lagrange multiplier) that must be calculated at the end by imposing the constraint qN=Pixi=L. Therefore, the extensions xi’s obey the deterministic equations γ˙x1=−a0 1(x1) + a0 2(x2),(2.7a) γ˙xi=−2a0 i(xi) + a0 i+1(xi+1) + a0 i−1(xi−1),1< i < N, (2.7b) γ˙xN=−2a0 N(xN) + a0 N−1(xN−1) + F, (2.7c) F=γvp+a0 N(xN).(2.7d) We have introduced the pulling speed vp≡˙ L, (2.8) which is usually time independent. We assume that ai(xi) allows for bistability in a certain range of the external force F, in the sense that ai(xi)−Fxiis a double-well potential with two minima, see figure 2.2. Therefore, in that force range, each unit may be either folded, if xiis in the well corresponding to the minimum with the smallest extension, or unfolded, when xibelongs to the well with the largest extension. If the length is kept constant (vp= 0), there is an equilibrium solution of (2.7), a0 1(xst 1) = a0 2(xst 2) = ··· =a0 N(xst N) = Fst,(2.9) and Fst is calculated with the constraint Pixst i=L. This solution is stable as long as a00 i(xst i)>0 for all i. If all the units are identical, ai(x) = a(x), the metastability regions of each module— the range of forces for which the equation a0 i(x) = Fhas several solutions—coincide. Therefore, as briefly introduced in section 1.1.2, we obtain stationary branches corresponding to Junfolded units and N−Jfolded units that have been analyzed in detail in [39,41]. If all the modules are not identical, the metastability regions do not perfectly overlap since the units are not equally strong: the weakest one is that for which the equation a0 i(x) = Fceases to have multiple solutions at a smaller force value. It is important to note that if we change all the forces a0 i(xi) to ˜a0 i(xi) = a0 i(xi)−F0 and Fto ϕ=F−F0, we have the same system (2.7) but with ˜a0 iand ϕinstead of a0 iand F, respectively. Then, we may use the free energies for any common value of the force F0and interpret the Lagrange multiplier as the excess force from this value to be applied to the system. A similar result is also found if the length is controlled by using a device with a finite value of the stiffness kc. A constant force only shifts the equilibrium point of a harmonic oscillator: (qN−L) must be substituted by (qN−L−F0/kc). 2.2 Pulling the system: perturbative solution We write the i-th unit’s free energy as ai(x) = a(x) + ξ δai(x),(2.10) i i “thesis” — 2018/9/7 — 8:45 — page 33 — #41 i i i i i i Section 2.2. Pulling the system: perturbative solution 33 in which a(x) is the “main” part, common to all the units, and ξδai(x) represents the separation from this main contribution. If all the units are perfectly identical, ai(x) = a(x) for all ior, equivalently, δai(x) = 0. In principle, in an actual experiment, the splitting of the free energy in (2.10) can be done if the free energy aiof each unit is known: we may define the common part as the “average” free energy over all the units, a(x)≡a(x)≡N−1PN i=1 ai(x), and ξδai(x)≡ai(x)−a(x). From a physical point of view, the dimensionless parameter ξ > 0 measures the importance of the heterogeneity in the free energies. Our theory could be applied to a situation in which the free energy deviations δaiwere stochastic and followed a certain probability distribution, for instance to represent the slight differences among very similar units, as done in [41] to analyze the force-extension curves. Also the forces a0 i(x) in the evolution equations are split as a0 i(x) = a0(x) + ξ δfi(x), δfi(x)≡δa0 i(x).(2.11) As already noted above, we can use the free energies for any common value of the force F0, and interpret Fas the extra applied force from this value. In what follows, we consider the main part a(x) with two, equally deep, minima corresponding to the folded (F) and unfolded (U) configurations. Therefore, our “origin of force” F0corresponds to the critical force for the main, common, contribution ato the units’ free energies. Figure 2.2 presents a qualitative picture of the free energy and its derivative. The two minima correspond to lengths `Fand `U, with `F< `U. Also the point `bat which a00(`b) = 0 is marked. It is the condition a00(`b) = 0 that essentially determines the stability threshold, as it provides the limit force Fb=a0(`b)>0 at which the folded basin ceases to exist for the “main” potential. In the deterministic approximation considered here, thermal fluctuations are neglected and, for F < Fb, the folded unit cannot jump over the free energy barrier hindering its unfolding: it has to wait until, at F=Fb, the only possible extension is that of the unfolded basin. Of course, neglecting thermal noise restricts in some way the range of applicability of our results, see section 2.3 for a more detailed discussion and also the numerical section 2.4.2. Keeping the above discussion in mind, now we analyze the limit of stability of the different units. The asymmetry correction δfishifts the threshold force for the different units and the extension xi,b at which the i-th unit loses its stability is obtained by solving the equation a00 i(xi,b) = a00(xi,b) + ξδf0 i(xi,b) = 0. Linearizing in both the displacement xi,b −`band ξone gets a00(`b) + a000(`b)(xi,b −`b) + ξδf0 i(`b) = 0.(2.12) Noting that a00(`b) = 0, we obtain xi,b =`b−ξδf0 i(`b) a000(`b).(2.13) The corresponding force is Fi,b ≡a0 i(xi,b) = Fb+ξ δfi(`b),(2.14) in which we have consistently dropped terms of the order of ξ2. Then, units with δfi(`b)< 0 (δfi(`b)>0) are weaker (stronger) than average. See appendix A for more details. i i “thesis” — 2018/9/7 — 8:45 — page 40 — #48 i i i i i i 40 Chapter 2. The basics of modelling modular systems computed in a more direct way,1 vend c= max jvj(N).(2.34) Consistency of the theory requires that v(k+1) c> v(k) c. This can be proved right away. The consistency condition implies that δfαk+2 (`b)>δfαk+1 (`b)(νk+2 −νk)−δfαk(`b)(νk+2 −νk+1) νk+1 −νk ,(2.35) in which νk=αk(αk−1). Due to (2.33), αk+1 minimizes vαk(j). Therefore, in particular, vαk(αk+1)< vαk(αk+2), which is readily shown to be equivalent to (2.35) and proves the inequality. We have a trivial case for α1=N, when the pulled unit is precisely the weakest and it is always the first to unfold for any pulling speed. The simplest nontrivial case appears when all the modules has the same free-energy with the exception of the weakest, and α16=N, (2.29), (2.31) and (2.34) reduce to γv(1) c ξ=γvend c ξ=2N[δfN(`b)−δfα1(`b)] N(N−1) −α1(α1−1).(2.36) Note that the situation is quite simple, since there exist a single critical velocity vc= v(1) c=vend c. For vp< vcthe weakest module unfolds first whereas for vp> vcthe last one unfolds first. If more units have different free energies, the situation may be more complex, as shown in the previous paragraph. There appear intermediate critical velocities, which define pulling speed windows where neither the weakest unit nor the last one is the first to unfold. In order to obtain these regions, we need to recursively evaluate (2.33). 2.4 Comparing our theory with simulations 2.4.1 Free energies of the units: shape and physical parameters Different shapes for the double-well potentials have been considered in the literature. They can be classified in, mainly, two different classes: simple Landau-like quartic potentials which are employed to understand the basic mechanisms underlying the observed behaviors [17,39,41], and more complex realistic potentials, when the aim is obtaining a more detailed, closer to quantitative, description of the experiments [40,41,121,122]. For the simplest modeling of a double well, we use a quartic potential aq(x) aq(x) = ε 4h(x−σ)2−α2i2.(2.37) The physical meaning of the parameters are straightforward: εscales the shape of the potential, x=σgives the position of the maximum of the barrier, whereas the minima 1In our original publication [120], this equation had a typo, specifically an extra “min” in the subindex. i i “thesis” — 2018/9/7 — 8:45 — page 41 — #49 i i i i i i Section 2.4. Comparing our theory with simulations 41 are at x=σ±α. This simple dependence allows us to compute analytically the borders of the metastability region, specifically one can obtain the stability threshold `b=σ−α √3,(2.38) and the forces within the metastability region fulfill |F|< Fb=a0(`b) = 2√3 9εα3.(2.39) Below, we use a nondimensional version of this potential. It is done by defining a length scale [x] and a energy scale ε[x]4, aq(x) = 1 4h(x−σ)2−α2i2(2.40) where we have not used any special notation for the nondimensional variables in order not to clutter our formulae. The more realistic proposal we use in this work has been put forward by Berkovich et al. some years ago [121,122]. Therein, the free energy of a module is represented by the sum of a Morse potential, which mimics the enthalpic minimum of the folded state, and a WLC term [29], which accounts for the entropic contribution to the elasticity of the unfolded state. Specifically, Berkovich’s free energy aB(x) is written as aB(x) = U01−e−2bx−Rc Rc2−1+kBT 4PLc 1 1−x Lc−1−x Lc +2x2 L2 c!.(2.41) This shape has shown to be useful for some pulling experiments with actual proteins as titin I27 or ubiquitin [121, 122]. Therein, each parameter has a neat physical interpretation. First, in the WLC part, we recall (see section 1.1.2) that we have: (i) the contour length Lc, which is the maximum length for the totally extended protein, and (ii) the persistence length P, which measures the characteristic length over which the chain is flexible. Both of them, Lcand P, can be thought in terms of the number of amino acids in the chain. Note that the derivative of the WLC term leads precisely to the force-extension curve given by (1.6), with the identification L↔x. Second, for the Morse contribution, we have: (iii) Rc, which gives the location of the enthalpic minimum and (iv) U0and b, which measure the depth and the width (in a nontrivial form) of the folded basin. The explicit function giving the stability threshold `bin terms of the parameters in Berkovich’s potential cannot be obtained. However, we can always estimate `bnumerically, solving a00 B(`b) = 0 for a specific set of parameters. Once more, we use below nondimensional variables. In order to do so, we define a force scale [F] and take Lcas the length unit. Accordingly, dimensionless variables are introduced with the definitions µ=U0/(Lc[F]), β= 2bLc/Rc,ρ=Rc/Lc,A= kBTLc/(4PU0). Thus, the corresponding dimensionless potential reads aB(x) = µh1−e−β(x−ρ)i2−1 + A1 1−x−1−x+ 2x2.(2.42) i i “thesis” — 2018/9/7 — 8:45 — page 42 — #50 i i i i i i 42 Chapter 2. The basics of modelling modular systems 2.4.2 Numerical results Here, we check the agreement between our theory and the numerical integration of the evolution equations. First, we discuss the validity of the simplifications introduced in the development of the theory, namely (i) negligible thermal noise and (ii) perfect length control. Second, we look into the critical pulling speed, showing that there appears such a critical speed in the simulations and comparing this numerical value with our theory. We consider a system composed of N= 4 units, such as the maltose binding protein [17], see section 1.1.3. First, each unit is assumed to be characterized by a quartic bistable free energy. In reduced variables, the free energies have the form ai(x) = ia(x), with i= 1 for i6= 1, 1<1, where a(x) = aq(x) as given by (2.40), with σ= 0 and α= 3. Here, the value of σis different from the one in [17] (σ= 8). Its only effect is a shift of the origin of the extensions, our choice implies that a positive (negative) sign of the extension corresponds to an unfolded (folded) configuration. The value of the friction coefficient is, also in reduced variables, γ= 1. We use these dimensionless reduced variables to make it easier to compare our results with those in [17]. As stated above, the function (2.40) is one of the simplest, but reasonable, choices to describe the free energy of different unit of the same modular biomolecule. Using the notation introduced in (2.11), we have δfi(x)=0, i 6= 1, δf1(x) = −ξa0 q(x),(2.43) with ξ= 1 −1. Equations (2.13) and (2.14) give us the limits of stability up to first order in the asymmetry ξ, xi,b =`b,∀i, Fi,b =Fb, i 6= 1, F1,b = (1 −ξ1)Fb.(2.44) For this simple example, (2.44) is exact. The weakest unit is the first one, because F1,b is the minimum value of the force at the limit of stability. For the values of the parameters we are using, `b=−α/√3 = −1.73 and Fb=a0(`b) = 2√3α3/9 = 10.4. Since we are writing the free energies for a common given value of the force, all the units have their two minima equally deep at the same force. This assumption is made to keep things simple: the main ingredient for having an unfolding pathway that depends on the pulling speed is to have different values of the forces Fi,b at the stability threshold for the different units. In the case we are considering, the weakest unit is the first one, while the others share the same free energy. This means that we have the simplest scenario for the critical velocity in our theoretical approach: either the weakest (for vp< vc) or the last (for vp> vc) unit that opens first, as discussed at the end of the section 2.2. Here, (2.36) for α1= 1 and N= 4 reduces to γvc ξ=2 3Fb.(2.45) To start with, we consider the relevance of the noise terms in (2.2). In figure 2.3, we plot the integration of the Langevin equations together with the deterministic approximation [95] for a particular case: the first unit’s free energy corresponds to 1= 0.8 (ξ= 0.2), the stiffness of the length control device is kc= 5, the temperature is T= 1, and the pulling speed is vp= 0.38. For these values of the parameters, taken from [17], the critical velocity in (2.45) is vc= 1.4; thus we are considering a subcritical velocity, i i “thesis” — 2018/9/7 — 8:45 — page 43 — #51 i i i i i i Section 2.4. Comparing our theory with simulations 43 Figure 2.3: Evolution of the units’ extensions as a function of the system length qN. The pulling speed is vp= 0.38 and the length control device has a stiffness kc= 5. The symbols correspond to a typical realization of the Langevin process (2.2) with T= 1, whereas the lines correspond to the deterministic (zero noise) approximation. vp< vc. Thermal fluctuations are small, which entail that the same unfolding pathway is observed in the deterministic and the majority of the stochastic trajectories. Let us consider in more detail the relevance of thermal noise: from a physical point of view, it may be inferred by looking at the height of the free energy barrier at the critical force in terms of the thermal energy kBT. For the values of the parameters we are using, this barrier is around 20kBT, which explains why thermal noise is basically negligible in figure 2.3. If the temperature is decreased from T= 1 to T= 0.25, the barrier is so high, around 80 times the thermal energy, that essentially all the stochastic trajectories coincide with the deterministic one. On the other hand, if the temperature is increased to T= 4, the barrier in only a few kBT, and we expect that the deterministic approximation ceases to be valid. In order to further clarify the role played by the temperature, we present figure 2.4. Both panels display bar graphs with the frequencies with which each unit unfolds first in the stochastic trajectories. Specifically, the statistics shown has been obtained with 1000 trajectories of the Langevin equations (2.2) with perfect length control, and several different values of the temperature. In the top panel, a subcritical velocity vp= 0.38 < vc is considered, so that the weakest (first) unit is expected to unfold first. In the bottom panel, the numerical data for a supercritical velocity vp= 2 > vcare shown, for which the pulled (fourth) unit would unfold first. The effect of thermal noise is quite similar in both cases. For the lowest temperature T= 0.25, the frequency of the deterministic pathway is close to unity and, for the temperature in figure 2.3, T= 1, its frequency is i i “thesis” — 2018/9/7 — 8:45 — page 44 — #52 i i i i i i 44 Chapter 2. The basics of modelling modular systems Figure 2.4: First-unfolding frequency for each of the units. These frequencies have been obtained by integrating the Langevin equations with perfect length control, for different values of the temperature. (Top) Numerical frequencies obtained in 1000 trajectories, for a subcritical pulling speed vp= 0.38 < vc. (Bottom) The same as in the top panel, but for a supercritical pulling speed vp= 2 > vc. As the temperature decreases, the frequency of the deterministic unfolding pathway approaches unity in both cases. i i “thesis” — 2018/9/7 — 8:45 — page 45 — #53 i i i i i i Section 2.4. Comparing our theory with simulations 45 still very large, clearly larger than any of the others. On the other hand, for the highest temperature, T= 4, thermal noise is no longer negligible. In the following, we restrict the analysis to the physically relevant case in which the deterministic approximation gives a good description of the first unfolding event. In figure 2.5 (top panel), we look into the same pulling experiment as before, but now we compare the deterministic evolution of the extensions for two finite values of the stiffness to the kc→ ∞ limit. Consistently with our expectations, the unfolding pathway is not affected by this simplification. For the smaller values of kc, the length is not perfectly controlled, but the length control improves as kcincreases, as seen in the bottom panel. Notwithstanding, the curves in the top panel, which correspond to different values of kc, are almost perfectly superimposed when plotted as a function of the real length of the system qN(but not of the desired length L). This means that the real length qNis a good reaction coordinate, as already stated in section 2.2. We have integrated the deterministic approximation (2.7) (zero noise) of the Langevin equations for different values of the pulling speed, and extracted from them the numerical value of the critical velocity as a function of the asymmetry ξ= 1−1. In order to obtain this numerical prediction, we initially set vpequal to the theoretical critical velocity given by (2.45). Then, we recursively shift it by a small amount δvp, such that δvp/vc= 0.0001, until the pathway changes. We compare the values so obtained to the theoretical expression (2.45), in figure 2.6. We find an excellent agreement for ξ.0.1, for ξ > 0.1 there appear some quantitative discrepancies. They stem from two points: (i) the perturbative expansion used for obtaining (2.36) from (2.27) and (ii) the intrinsically approximate character of (2.27), since `bgives rigorously the limit of stability only for the static case vp= 0. Therefore, we have looked for the solution of (2.27) in the numerical integration of the deterministic equations. This is the dashed line in figure 2.6, which substantially improves the agreement between theory and numerics because we have eliminated the deviations arising from point (i) above. In fact, for the case we have studied in the previous figures, which corresponds to a not so small asymmetry ξ= 0.2, the improved theory gives an almost perfect prediction for the critical velocity. Note that the new numerical estimate of vcby solving (2.27) is always below that given by (2.36). This can be easily understood: since the value `cat which x1and x4intersect is lower than `b, the pulling velocity needed for the crossing is lower than (2.36). Now we consider the more realistic Berkovich’s potential, given by (2.41) and (2.42). We take the values of the parameters from [41,123], namely P= 0.4 nm, Lc= 30 nm, Rc= 4 nm, b = 2, U0= 100 pN nm,(2.46) and T= 300 K. Defining the force scale [F] = 100 pN, the values of the nondimensional parameters in (2.42) are µ= 0.0333, β= 30, ρ= 0.133 and A= 0.776. In dimensionless variables, Fb= 0.527 (52.7 pN) and `b= 0.157 (4.70 nm). The relevant time scale is set by the friction coefficient γ, [t] = γLc/[F]. In turn, γis given by the Einstein relation D=kBT/γ, where Dis the diffusion coefficient for tethered proteins in solution. We consider a typical value D= 1500 nm2/s, also taken from [121], so that γ= 0.0028 pN nm−1s. We consider a system of 4 units, again with all the units but the first being identical. Then, ai(x) = aB(x), i6= 1, and the first unit is the weakest because a1(x) = (1−ξ)aB(x). The situation is then similar to the one we have already analyzed with the quartic potential i i “thesis” — 2018/9/7 — 8:45 — page 46 — #54 i i i i i i 46 Chapter 2. The basics of modelling modular systems Figure 2.5: (Top) Evolution of the units’ extensions as a function of the system length qN. The symbols correspond to the integration of the deterministic equations, for kc= 5 (filled) and kc= 50 (empty), whereas the line corresponds to the limit of perfect length control, kc→ ∞. The pulling speed is the same as in figure 2.3, that is, vp= 0.38. (Bottom) Comparison between the desired and actual lengths, Land qN. Clearly, the length control improves as kcincreases. i i “thesis” — 2018/9/7 — 8:45 — page 47 — #55 i i i i i i Section 2.4. Comparing our theory with simulations 47 Figure 2.6: Phase diagram for the unfolding pathway in the pulling velocity-asymmetry plane for the quartic potential. Two well-defined regions are observed. These regions are separated by the curve critical velocity vcvs. asymmetry ξof the first unit. The numerical values for vc(circles) are compared to the theoretical expression (2.45) (solid line). The dashed line corresponds to the alternative approach discussed in the text, which improves the agreement with the numerical results for ξ > 0.1. Error bars have been omitted because they are smaller than point size. (2.40), but there is a difference that should be noted: here, aB(x) is the free energy at zero force, whereas for the quartic potential aq(x) was the free energy at the critical force F0for which the folded and unfolded minima were equally deep. Then, the force here must not be interpreted as the extra force from F0, but as the whole force applied to the polyprotein. On the basis of our theory, we expect the simplest situation with only one critical velocity vc, below (above) which the weakest unit (the pulled unit) unfolds first. This is also indeed the case in the numerical simulations, and we compare the theoretical and numerical critical velocities in figure 2.7. A very good agreement is found again, up to values of the asymmetry ξof the order 0.1−0.2. The above discussion shows that the validity of the theory presented here is not restricted to simple potentials like the quartic one; on the contrary, it can be confidently applied to situations in which the units are described by realistic potentials. For the typical parameters we are using, the theoretical critical velocity vcfor the Berkovich potential equals 1270nm/s for an asymmetry ξ= 0.1. The latter can be regarded as a conservative estimate of the largest asymmetries for which our theory gives an almost perfect account of the unfolding pathway. Interestingly, this pulling speed corresponds to the upper range of velocities employed in AFM experiments, for instance see table I of [9]. Therefore, testing our theory in real AFM experiments with modular proteins should be achievable, see also section 2.6. Finally, we consider a more complex situation, in which more than one unit is different i i “thesis” — 2018/9/7 — 8:45 — page 48 — #56 i i i i i i 48 Chapter 2. The basics of modelling modular systems Figure 2.7: Phase diagram for the unfolding pathway in the pulling velocity-asymmetry plane for the Berkovich potential (2.42). As in figure 2.6, there appear two well-defined regions, separated by the curve critical velocity vcvs. asymmetry ξ. Numerical values for vc(circles) compare very well with our theoretical prediction (2.36) (solid line). Again, the dashed line corresponds to the alternative approach discussed in the text, which once more improves the agreement theory-simulation for the larger values of ξ. from the rest. Therefore, there may exist more than one critical velocity as discussed in section 2.3. To be concrete, we have considered a system with 4 units in which a2(x) = a3(x) = aq(x), a1(x) = (1 −ξ)aq(x), as before, but the pulled unit free energy is changed to a4(x) = (1 + 3ξ/2)aq(x). In this situation, we have two different critical velocities: for very low pulling speeds, the weakest unit is the first to unfold, but there appears a velocity window inside which neither the weakest nor the pulled unit is the first to unfold. This stems from the fact that the first and the third unit reach simultaneously the limit of stability for a velocity v1(3) = 4ξγ−1Fb/3 that is smaller than the velocity v1(4) = 5ξγ−1Fb/3 for which the first and the last would do so. The physical reason behind this is the pulled unit’s threshold force being larger enough than that of the third one. We recall that vi(j) is the velocity for which the i-th and the j-th unit reach simultaneously their limits of stability. Afterwards, the third unit and the fourth attain the limit of stability in unison for a velocity v3(4) = 2ξγ−1Fb, and the following picture emerges from our theory. Using the notation introduced in section 2.3, we define two critical velocities, γv(1) c ξ=4Fb 3,γv(2) c ξ= 2Fb,(2.47) such that: (i) for vp< v(1) c, it is the weakest unit that unfolds first, (ii) for v(1) c< vp< v(2) c, i i “thesis” — 2018/9/7 — 8:45 — page 49 — #57 i i i i i i Section 2.4. Comparing our theory with simulations 49 Figure 2.8: Phase diagram for the unfolding pathway in the pulling velocity-asymmetry plane for the more complex with two critical velocities. Now, we have three well-defined regions, separated by two curves, respectively, critical velocities v(1) c(blue) and v(2) c(black) vs. asymmetry ξ. Note that our theoretical approach is able to reproduce the existence of the three different pulling regimes. Nevertheless, the discrepancies between theoretical and numerical values for the critical velocities are larger than in figures 2.6 and 2.7. it is the third unit that unfolds first, and (iii) for vp> v(2) c, the first unit to unfold is the pulled one. We check the more complex scenario described in the previous paragraph in figure 2.8. On the one hand, our theory correctly predicts the existence of the three pulling regimes described above. On the other hand, even for very small asymmetries, there appear some noticeable discrepancies between theory and simulation. The validity of the perturbative expansion for obtaining the critical velocities, expressed by condition (2.27), is strongly supported by the accurateness of the theoretical prediction for the simplest case with only one critical velocity, see figures 2.6 and 2.7. Then, we believe that this discrepancy stems from the intrinsically approximate character of the condition a00 i= 0 for determining the stability threshold for finite pulling velocity vp6= 0. Thus, improving the present theory should involve the derivation of a more accurate condition for the stability threshold in this case. This refinement to our theoretical framework, which probably makes a multiple scale analysis necessary for lengths close to the condition a00 = 0, is an open question that deserves further investigation. i i “thesis” — 2018/9/7 — 8:45 — page 56 — #64 i i i i i i 56 Chapter 2. The basics of modelling modular systems Figure 2.12: (Top) Force-extension curve corresponding to the pulling experiment in figure 2.11. Specifically, we consider two pulling velocities vp= 10 nm/s (subcritical, solid) and vp= 22 nm/s (supercritical, dashed). (Bottom) Zoom of the region of interest, which clearly shows the shift between the peaks stemming from the increased contour length of the mutant unit. i i “thesis” — 2018/9/7 — 8:45 — page 57 — #65 i i i i i i Chapter 3 Testing the model with molecular dynamics The theoretical approach developed in the previous chapter is clearly a drastic simplification of reality. In fact, the final result of our theory is deterministic, in the sense that the unfolding pathway is a definite one, the randomness coming from thermal fluctuations being effectively “suppressed” by the fast enough pulling velocity. In reality, as also discussed in chapter 2, the unfolding pathway does have some degree of stochasticity, stemming from the interactions between the molecule under study and the fluid where it is immersed, which are encoded in the Gaussian white noises of our Langevin description. Testing the theory is a mandatory step of scientific method. The task of performing real pulling experiments with a modular protein that matches all the requirements of our theoretical framework is not simple at all. Therefore, in this chapter, we focus on steered molecular dynamic (SMD) of engineered systems. Since the first investigations made by Grubm¨uller [126] and Schulten [127], these computational techniques have shown to be of crucial relevance in the current development of biophysics [128–130]. To test our theoretical predictions, we consider a particularly simple structure composed of coiled coils. This kind of structures is common in nature, which makes it extremely useful as a model system [131–133]. The organization of the chapter is as follows. Section 3.1 is devoted to stress the relevance of the hypothesis assumed in the theory developed in the previous chapter. Also, we discuss its applicability to specific molecules in specific ranges of velocities. Therein, we introduce the construct we work with throughout this chapter, which comprises two consecutive coiled-coil structures. The simulation procedure along with the method employed for the data analysis are presented in section 3.2. Specifically, we explain how we acquire the initial configurations from which the construct is pulled. In section 3.3, the results of our simulations are put forward, and they are compared with our theoretical prediction. A key point in our numerical analysis is our assumption about the independence of the initial conditions chosen for the pulling stage. Then, in section 3.4, we prove this data acquisition to be actually uncorrelated, within the accuracy of our statistics. 57 i i “thesis” — 2018/9/7 — 8:45 — page 58 — #66 i i i i i i 58 Chapter 3. Testing the model with molecular dynamics 3.1 The candidate protein We expect our theoretical approach to hold for some molecules within a specific range of pulling velocities. One of the obvious requirements our protein candidate must meet is a negligible interaction between repeats, since we have assumed no nearest-neighbor interaction terms in the global free energy. Regarding the range of velocities, if pulling is very slow and quasistatic, the first unfolding event occurs at the length value for which the free energy minima over the branch with all the units folded and over the branch with only one unit unfolded are equally deep [39]. Therein, the jump between branches occurs by thermal activation over the free energy barrier separating them. As a consequence, the completely folded branch is only partially swept, as marked by the dashed vertical (red) lines in figure 3.1. In contrast, there is a range of fast pulling velocities that do not give the system enough time to be thermally activated over the barrier, but are slow enough to allow it to sweep completely the part of the branches that corresponds to metastable equilibrium states. In this case, the jump between branches comes about at the limit of metastability, only when the folded minimum disappears [41], as depicted in the top panel of figure 3.1. This is marked by the solid vertical (blue) lines in the bottom panel of figure 3.1. This range of adiabatic velocities we are interested in, as defined in chapter 2, are also said to lead to the “maximum hysteresis path” [42], regime briefly introduced at the end of section 1.1.2. We have designed a simple homopolyprotein, which we employ below to test whether it fits our theoretical description. We have extracted the structure of an antiparallel coiled-coil motif (CC) from the archeal box C/D sRNP core protein (Protein Data Bank entry 1NT2), which comprises 67 residues and whose N-terminus and C-terminus are, respectively, arginine (ARG) and isoleucine (ILE) [134]. This structure has been proven to be useful as a mechanical folding probe [124]. We use this CC as the building blocks of the molecule: our system is simply a concatenation of two CC motives connected by a linker, which is composed of two consecutive pairs of alternated residues of glycine and serine. We expect this linker not to introduce any significant interaction between the two domains. The initial conformation of the constructed model structure and orientation of the two CC repeats is shown in figure 3.2. The end- to-end vector points from the N-terminus to the C-terminus, aligned with the x-axis, whereas both axial directions of the two CC structures are located as parallel as possible to the z-axis. According to the theoretical framework we have developed along this work, since we have two identical units, we expect that if we pull from one end of the designed molecule, the first repeat to unfold will be precisely the closest to the moving end. We perform SMD simulations to analyze the degree of agreement of the obtained numerical results with the theory. 3.2 All-atom molecular dynamics simulation Our molecular dynamics simulations start from the initial conformation shown in figure 3.2. First, we add hydrogen atoms using VMD Automatic PSF Builder [135]. Then, we create a water box of size 300˚ A×70˚ A×120˚ A, in the x-, y-, and z-axes, respectively, which is long enough in the direction of pulling (the x-axis) to contain the unfolded protein. i i “thesis” — 2018/9/7 — 8:45 — page 59 — #67 i i i i i i Section 3.2. All-atom molecular dynamics simulation 59 Figure 3.1: (Top) Schematic free energy landscape of a single repeat for three different forces F1< F2< F3. The system starts in the folded (F) state, which is the absolute minimum for F1. If thermal noise is negligible, the repeat remains in the folded state for F2even when the unfolded state is more stable. For F3, the F state disappears, and the repeat finally unfolds. (Bottom) Qualitative picture of the stability branches in a modular system with two units. The blue line shows the unfolding pathway followed in the limit of the so-called “maximum hysteresis path”, when the pulling speed is high enough to make the system sweep the whole branches, including their metastable parts. In this limit, the jumps between consecutive branches take place by the mechanism shown in the top panel because the system does not have enough time to jump over the barrier separating the folded and unfolded states. In other words, the “fast enough” pulling speed effectively suppresses thermal fluctuations. Conversely, in the quasistatic limit the transition from folded to unfolded takes place at the lengths (dashed red lines) at which the branch with one more unfolded unit becomes more stable, that is, when its free energy becomes smaller. For the quasistatic case, the system has always time to find its way through the barrier. i i “thesis” — 2018/9/7 — 8:45 — page 60 — #68 i i i i i i 60 Chapter 3. Testing the model with molecular dynamics Figure 3.2: Initial conformation of the homopolyprotein comprising two CCs in the SMD simulations. The pulling direction is aligned with the x-axis, whereas the axial directions of the CCs are aligned with the z-axis. Also, NaCl is introduced in the system replacing water molecules until the concentration reaches 150 mM/L and the charge is neutralized. Finally, simulations are performed using NAMD2 2.10 [136], with two different stages: (i) the “equilibration stage” at 310 K and (ii) the “pulling stage” with velocity v0= 1.4·10−2nm/ps and stiffness 4860 pN/nm. We have also considered faster pulling velocities, namely 2v0and 5v0, as detailed below. The molecule’s behavior for these faster pulling velocities has been investigated in order to elucidate whether or not the unfolding pathway becomes more deterministic as the pulling speed is increased. These typical pulling speeds in steered molecular dynamics simulations are higher by several orders of magnitude than the experimental ones. However, they are necessary to investigate this kind of system with the available computer power. In addition, this high velocity range is especially relevant for our present purposes, since we are interested in exploring the maximum hysteresis path limit. Note that the considered value for the stiffness of the elastic reaction is also two orders of magnitude higher than the typical ones in AFM experiments, and thus closer to the perfect length control situation assumed in the theory developed in chapter 2. 3.2.1 Pulling trajectories A notable number of pulling trajectories NTare needed in order to obtain a meaningful statistical analysis of the unfolding pathway. The final configuration of the molecule in the equilibration stage is taken as the initial condition for the pulling stage. In order to generate different initial conditions for pulling, we have considered one “long” trajectory in the equilibration stage and collected the molecule configurations at several different times i i “thesis” — 2018/9/7 — 8:45 — page 61 — #69 i i i i i i Section 3.2. All-atom molecular dynamics simulation 61 tk, with tk+1 −tk>0.1 ns, as the initial conditions for the different pulling trajectories k= 1, . . . , NT. We have checked the “independence” of the trajectories obtained from these initial conditions, in the sense that the unfolding pathway from two consecutive initial conditions, corresponding to kand k+ 1, are not correlated. Details are given in section 3.4. The duration of the pulling stage ∆tpis chosen to allow the molecule to unravel. The size of a single CC motif in its axial direction is around 5 nm, specifically 4.82 nm between the two Cαs most separated in the axial direction. Therefore, a motif can be considered as completely unfolded when its end-to-end distance, measured between their Cαs in the terminal residues ARG and ILE, exceeds 10 nm. For the “base” velocity v0= 1.4·10−2nm/ps, we have chosen ∆tp= 1.6 ns, so that the total length increment is v0∆tp= 22.4 nm. For the faster pulling velocities, 2v0and 5v0, we have decreased ∆tp accordingly. In SMD simulations, the length of each repeat can be measured as a function of time and thus we may introduce the basic (and the simplest) classification of trajectories by labeling them as “good” (G) if the CC motif closest to the pulled end unfolds sooner than the furthest, and “bad” (B) otherwise. Clearly, it is G-trajectories that agree with the prediction of the theory developed in chapter 2, when particularized for identical units. The above basic classification of SMD trajectories as G or B can be refined, so as to have a more accurate description of the unfolding pathway, as done below in section 3.2.2. In our theoretical framework, the unfolding is sequential: when the first unit unfolds, the second one has not reached its limit of stability and thus remains in the folded state. However, in the simulations, this perfectly sequential unfolding is not always found: when the first unit unfolds, the second one can be partially unfolded, see also section 3.3. 3.2.2 Unfolding criterion Following the discussion in the previous paragraph and for the sake of accuracy, we define a criterion for distinguishing between different subtypes of trajectories. Therefore, we incorporate a quantitative measurement of the degree of unfolding for each repeat; first in a geometric way, and second taking into account the fraction of native contacts [137]. As already said above, the size of the CC motif in its axial direction is approximately 5 nm and thus we consider a motif to be completely unfolded when its end-to-end distance is greater than 10 nm. Also, we define when a motif is partially unfolded in our simulations by introducing an “unfolding threshold”, that is, a length below which we consider the unit to be still folded. Specifically, we take this unfolding threshold to be 7 nm, which corresponds to an opening angle of 90◦in a rigid rods picture, as depicted in figure 3.3. The above “geometric” choice for the unfolding threshold has some degree of arbitrariness, especially in relation to the length (or angle between rigid-rods) chosen for the unfolding threshold. In order to give a physical basis for this choice, we look into the fraction of native contacts [137] as a function of the total length for a single CC motif in figure 3.4. For this purpose, we have performed a SMD simulation in which a single CC motif, which is initially folded, is pulled at a speed of 3.75·10−3nm/ps from either its C-terminus (blue) or its N-terminus (red). It can be observed how the number of native contacts decreases along the trajectory, with a well-defined plateau arising between (approximately) 7 and 10 nm. Clearly, the borders of this plateau demarcate the region where the main i i “thesis” — 2018/9/7 — 8:45 — page 62 — #70 i i i i i i 62 Chapter 3. Testing the model with molecular dynamics Figure 3.3: Simple unfolding criterion for the CC under study. This criterion is based on a rigid-rod picture. Depending on the value of its end-to-end distance, we consider the molecule as (i) folded if it is shorter than 7 nm, (ii) partially unfolded when it is between 7 and 10 nm, and (iii) completely unfolded if it is longer than 10 nm. Table 3.1: Definition of the different trajectory types in a SMD C-pulling simulation. Types I and IV are the closest trajectories to a deterministic pathway, agreeing and disagreeing, respectively, with the prediction of our model. Type Subtype First repeat that unfolds State of the other repeat Length of the “folded” repeat Good (G) I C-terminus folded <7 nm Good (G) II C-terminus partially unfolded >7 nm Bad (B) III N-terminus partially unfolded >7 nm Bad (B) IV N-terminus folded <7 nm bonds that keep the double-stranded CC folded are broken. Therefore, the above-defined geometric thresholds for considering the molecule folded/partially unfolded/completely unfolded agree with the corresponding limits in the fraction-of-native-contacts picture. Consistently with the above-described criteria, let us index the different subtypes of trajectories by (I, II, III, IV), attending to their degree of agreement with the theoretical prediction. G-trajectories are split into I and II subtypes: when the pulled motif is the first that unfolds, the other unit can be either still folded (type I) or partially unfolded (type II). Similarly, B-trajectories are divided into III and IV subtypes: when the nonpulled motif unfolds first, the pulled unit can be either partially unfolded (type III) or still folded (type IV). The explicit distinction between the different cases, for a simulation in which the molecules is pulled from its C-terminus, is shown in table 3.1. Obviously, if the molecule is pulled from its N-terminus, the same classification of trajectories applies but with the role of the repeats C and N reversed. i i “thesis” — 2018/9/7 — 8:45 — page 63 — #71 i i i i i i Section 3.3. Pulling the CC construct 63 Figure 3.4: Percentage of native contacts as a function of the total length. Both the data from C-pulling (blue) and N-pulling (red), at a speed 3.75 ·10−3nm/ps are plotted. The borders of the plateau agree quite well with the thresholds up to (from where) we consider the molecule to be folded (unfolded) in the geometrical picture, marked with vertical dashed lines. 3.3 Pulling the CC construct In this section, we present the SMD results corresponding to the pulling of the CC construction described above. We consider the two possible experiments, pulling from the C-terminus (C-pulling) and from the N-terminus (N-pulling), separately. 3.3.1 C-pulling First, we pull the molecule from its C-terminus at the base velocity v0= 1.4·10−2nm/ps. We plot the evolution of the distance between the end terminals of each repeat in this C-pulling experiment in figure 3.5. The red line stands for the pulled repeat (C-terminus) whereas we plot in blue the length of the other repeat (N-terminus). It can be seen how the pulled repeat clearly unfolds first in type I. Although from different categories, types II, III and IV seem to share a common feature. In the initial part of the trajectory, it is the pulled repeat the fastest to lengthen but its unfolding comes to a standstill before being completed, and the second repeat takes advantage of this impasse to increase its extension. Due to thermal fluctuations, we do not expect to obtain a perfect agreement with our theory, but a preponderance of the deterministic (type I) trajectories. As discussed before, a “fast” pulling velocity is necessary to be in the “maximum hysteresis path” limit, in which our theory is expected to hold. Therefore, in addition to the base velocity v0= 1.4·10−2nm/ps, we have carried out simulations at velocity 2v0. Specifically, we have done 31 trajectories for each velocity and collected their statistical information in table 3.2, which is completely compatible with our theoretical expectation. There is already a preponderance of type I trajectories at the base velocity v0, at which they represent i i “thesis” — 2018/9/7 — 8:45 — page 64 — #72 i i i i i i 64 Chapter 3. Testing the model with molecular dynamics Figure 3.5: Representative plots of the different types of trajectories for SMD C-pulling simulations. Each panel corresponds to a given type, as labeled. Both repeats’ extensions are plotted: N-repeat (blue) and C-repeat (red). Our model predicts that the pulled repeat (C-terminus) is the first that unfolds. Videos corresponding to each of the trajectory types can be found in the supporting information to [138]. i i “thesis” — 2018/9/7 — 8:45 — page 65 — #73 i i i i i i Section 3.3. Pulling the CC construct 65 Table 3.2: Statistical analysis of the output of 31 runs of SMD C-pulling simulations of the two CCs construct. At the base velocity v0, there is already a significant preponderance of type I trajectories (45%) as compared to type IV ones (16%). Type I trajectories clearly prevail as the pulling velocity is increased to 2v0, at which its frequency is boosted up to 67%. In addition, no type IV trajectories are found at this higher velocity. This behavior clearly agrees with our theoretical prediction. Type Occurrence at v0(%) Occurrence at 2v0(%) I 45 68 II 16 19 III 23 13 IV 16 0 Good = I + II 61 87 Bad = III + IV 39 13 almost half of the total number of trajectories, against a reduced fraction, only 16%, of type IV trajectories. Moreover, the prevalence of type I trajectories increases with the pulling velocity, as expected in our theoretical framework: at velocity 2v0, more than two thirds of the SMD runs are of type I and there are no type IV trajectories. We can get additional insight by calculating the frequencies of G (sum of I+II) and B (sum of III+IV) trajectories. This simplification of the types of events shows, even more clearly, the preponderance of the deterministic path as the velocity is increased: the frequency of G-trajectories shows a clear increase from 61% at velocity v0to 87% at velocity 2v0(the frequency of B-trajectories decreases accordingly from 39% to 13%). This is the reason why we have not considered even higher velocities for C-pulling. The above results show that increasing the pulling speed effectively diminishes the relevance of thermal activation effects and makes the unfolding more “deterministic”, in the sense of increasing the prevalence of type I trajectories. It should also be stressed that detailed analysis of the rest of type trajectories highlights a branching from type I due to an impasse of the length of pulled repeat, as already stated at the beginning of this section. Wrapping things up, our theory seems to predict the unfolding mechanism displayed in the SMD of this CC homopolyprotein—at least, when it is pulled from the C-end. 3.3.2 N-pulling Of course, our one-dimensional theory is completely left-right symmetric, since the free energy only depends on the extensions. Therefore, if we perform the same kind of SMD simulations but pulling from the N-terminus, we expect similar results (within statistical errors for the limited number of trajectories). Thus, the unfolding should start, preponderantly, from the unit closer to the N-terminus. Nevertheless, we show below that the situation is more complex, which we understand as a signature of anisotropy in the considered molecule. Table 3.3 presents the statistics of the different types of trajectories for SMD simulations in which we pull from the N-terminus. In terms of the same base velocity v0= 1.4·10−2nm/ps, we have conducted numerical experiments at velocities v0, 2v0and 5v0. Specifically, the statistics have been obtained again from 31 runs for each pulling velocity. The situation is much more complex than for C-pulling, since for both v0and i i “thesis” — 2018/9/7 — 8:45 — page 72 — #80 i i i i i i i i “thesis” — 2018/9/7 — 8:45 — page 73 — #81 i i i i i i Chapter 4 Continuum limit, fluctuating hydrodynamics and finite size effects Throughout this part of the thesis, we mainly investigate a lattice model for granular media inspired by two previous models on the lattice, specifically those in [76,141]. The model differs from these two previous proposals in a few crucial aspects. First, in [141], the velocity field evolved under the enforcement of the so-called kinematic constraint, which is disregarded here. Second, in [76], only the energy field was considered, and therefore momentum conservation was absent. Lattice models have been of paramount importance for understanding rigorously the conditions needed to have a hydrodynamic description, both at the average [68, 69] and fluctuating [67] levels. Recently, fluctuating hydrodynamics has been employed to derive the large deviation function in the context of energy-conserving [70–73] and even in energy-dissipating [75–78] models. Both in the conservative and nonconservative case, momentum conservation has not been taken into account. This shortcoming may be relevant, since it is known that momentum conservation is linked to the appearance of long-ranged correlations in out-of-equilibrium systems [142,143]. More specifically, spatial long-range velocity correlations in the homogeneous cooling state can be partly explained by fluctuating hydrodynamics, but require a more refined treatment to be fully investigated. This chapter is dedicated to introduce our model, from the basics to more involved issues. Our approach intends to contain the essential ingredients to investigate granular fluids, but reducing the complexity of the mathematical framework. In this way, we expect to get a transparent interpretation of the physical results. In particular, we aim to elucidate the “perturbative” nature of the continuum limit and calculate the corrections thereto [144]. Such corrections give interesting information about the structure—in space and time—of the correlated granular fluctuations and reveal new phenomena, which are peculiar of inelastic collisions. The organization of this chapter is the following. In section 4.1, we briefly revise the main aspects of the model, focusing on its continuum, hydrodynamic-like, limit. A 73 i i “thesis” — 2018/9/7 — 8:45 — page 74 — #82 i i i i i i 74 Chapter 4. Continuum limit, fluctuating hydrodynamics and finite size effects more detailed description of the model can be found in [79, 145, 146]. The evolution of the one-particle distribution function and some physically relevant stationary states are thoroughly analyzed, respectively, in sections 4.2 and 4.3. Section 4.4 is devoted to the study of the fluctuations of the system, including a discussion of their relation with some instabilities that appear in the model. Finally, we present an exact solution for the twoparticle correlation function of the system in section 4.5, which is valid for arbitrary size and provides further insight into the instabilities. 4.1 Model and previous results 4.1.1 Evolution equations Let us consider a 1d lattice with Nsites. First, we define the dynamics in discrete time. After the p-th step of the dynamics, the particle at the l-th site has a velocity vl,p ∈ R, and the configuration for the system at time pis denoted as vp≡ {v1,p, ..., vN,p}. One individual trajectory of the stochastic process is built in the following way: the configuration of the system changes from time pto time p+ 1 because a pair of nearest neighbors (l, l + 1) is chosen at random and collides inelastically, that is, vp+1 =ˆ blvp where the operator ˆ bltransforms the pre-collisional velocities (vl,p, vl+1,p) into the postcollisional ones (vl,p+1, vl+1,p+1) and leaves all other sites unaltered. The post-collisional velocities are given by vl,p+1 =vl,p −1 + α 2∆l,p,(4.1a) vl+1,p+1 =vl+1,p +1 + α 2∆l,p,(4.1b) where ∆l,p =vl,p −vl+1,p and the normal restitution coefficient is α∈(0,1]. Note that, by its own definition, this stochastic process is Markovian. In the following, we use a notation such that the evolution operator ˆ blacts naturally on observables, for example, ˆ blvl,p =vl,p+1. Momentum is always conserved, (ˆ bl−1)(vl,p +vl+1,p) = 0,(4.2) whereas energy, if α6= 1, is not, (ˆ bl−1)(v2 l,p +v2 l+1,p)=(α2−1)∆2 l,p/2<0.(4.3) The collision rule (4.1), which corresponds to the simplest one used in granular fluids [46], is valid for bulk sites. It must be complemented with suitable evolution equations for the sites next to the system boundaries, which depend on the physical situation at hand. See below for details. The evolution equation for the velocities can be cast in the form vl,p+1 −vl,p =−jl,p +jl−1,p, jl,p =1 + α 2∆l,pδyp,l,(4.4) which is nothing but a discrete continuity equation. Therein, jl,p is the momentum current from site lto site l+ 1 at time p,δyp,l is Kronecker’s delta, and ypis a homogeneously i i “thesis” — 2018/9/7 — 8:45 — page 75 — #83 i i i i i i Section 4.1. Model and previous results 75 distributed random integer in [1, L], where Lis the number of possible colliding pairs. For periodic boundary conditions, L=N, whereas for thermostatted boundaries L=N+ 1. We have only kinetic energy, Kp=PN l=1 el,p at time p, where el,p =v2 l,p. By squaring (4.4), the evolution equation for the energy at site lreads el,p+1 −el,p =−Jl,p +Jl−1,p +dl,p.(4.5) There are two contributions to the evolution of the energy: (i) the “flux” term −Jl,p + Jl−1,p, and (ii) a sink term dl,p stemming from the inelasticity of collisions. The energy current Jl,p from site lto site l+ 1 and energy dissipation dl,p at site lare Jl,p = (vl,p +vl+1,p)jl,p, dl,p =α2−1 4δyp,l∆2 l,p +δyp,l−1∆2 l−1,p<0,(4.6) respectively. The above stochastic dynamics generates the trajectories corresponding to the Markov process described by the following master equation in continuous time [145,146], ∂τPN(v, τ|v0, τ0) = ω L X l=1 |∆l|β"PN(ˆ b−1 lv, τ|v0, τ0) αβ+1 −PN(v, τ|v0, τ0)#,(4.7) where PN(v, τ|v0, τ0) is the conditional probability density of finding the system in state vat time τprovided it was in state v0at time τ0. Above, ωis a constant with dimensions of frequency that determines the time scale. Moreover, the operator ˆ b−1 lis the inverse of ˆ bl, that is, ˆ b−1 lchanges the post-collisional velocities into the pre-collisional ones for the colliding pair (l, l+1). At the p-th step of each dynamical trajectory, the continuous time τis incremented by δτp=−Ωp(L)−1ln z, Ωp(L) = ω L X l=1 |vl,p −vl+1,p|β,(4.8) in which zis a stochastic variable homogeneously distributed in the interval (0,1) and β≥0 is a parameter that affects the collision rate. On the one hand, the collision rate becomes independent of the relative velocity for β= 0, similarly to the case of pseudo- Maxwell molecules [141,147]. On the other hand, β= 1 and β= 2 are analogous to the hard core [83] and “very hard-core” [82, 148] collisions, respectively. From now on, we focus on the Maxwell case β= 0, although the case with β6= 0 is considered in chapter 6. The initial condition for (4.7) is clearly PN(v, τ0|v0, τ0) = δ(v−v0). Moreover, the one-time probability distribution PN(v, τ) verifies the same equation but with an arbitrary (normalized) initial condition PN(v,0). 4.1.2 Physical interpretation Literally taken, there is no mass transport in the model, particles are at fixed positions on the lattice and they only exchange momentum and kinetic energy. As discussed in i i “thesis” — 2018/9/7 — 8:45 — page 76 — #84 i i i i i i 76 Chapter 4. Continuum limit, fluctuating hydrodynamics and finite size effects Figure 4.1: Sketch of the lattice model of granular media. Sites are characterized by a velocity transverse to the lattice and exchange momentum and energy through nearest-neighbor collisions. section 1.2.1, this can be a valid assumption in an incompressible regime which is expected when the velocity field is divergence free, for instance during the first stage of the development of the shear instability, or in the so-called Uniform Shear Flow. We are also disregarding the so-called kinematic constraint, at difference with the approach in [141,149]: a colliding pair is chosen regardless of the sign of its relative velocity. Our model has a nice physical interpretation: the dynamics occurs inside an elongated 2d or 3d channel, the lattice sites represent positions on the long axis, while the transverse (shorter) directions are ignored; the velocity of the particles do not represent their motion along the lattice axis but rather along a perpendicular one, as depicted in figure 4.1. Let us easily imagine that the (hidden) component along the lattice axis is of the order of the perpendicular component, but in random direction. On the one hand, this justifies our disregarding of the kinematic constraint, while on the other, the collision rate may still be considered proportional to some power βof the relative velocity. A fair confirmation of this interpretation comes from the average hydrodynamics equations derived below. As anticipated in section 1.2.1, they replicate the transport equations (1.13) for granular gases in d > 1, but restricted to the shear (transverse) velocity field. 4.1.3 Average fields in the continuum limit In the large system size limit as L→ ∞, a continuum limit may be introduced by considering that the average velocity ul,p =hvl,piand energy El,p =hv2 l,piare smooth functions of space and time. Of course, the local temperature Tl,p =El,p −u2 l,p is also assumed to be smooth. Specifically, we introduce hydrodynamic continuous space and time variables, respectively,1 x=l−1 L, t =ωτ L2.(4.9) For β= 0, the balance equations for the average velocity u(x, t) and energy E(x, t) = u2(x, t) + T(x, t) (4.10) 1Note that the definition of the continuous space variable below is slightly different from that given in our original publications [145, 150]. This choice has no consequences for the equations derived in the continuum limit for the bulk sites, but allows us to simplify the derivations of some specific boundary conditions. i i “thesis” — 2018/9/7 — 8:45 — page 77 — #85 i i i i i i Section 4.1. Model and previous results 77 read ∂tu(x, t) = −∂xjav(x, t),(4.11a) ∂tE(x, t) = −∂xJav(x, t) + dav(x, t).(4.11b) Therein, the average momentum and energy currents, jav(x, t) and Jav(x, t), respectively, are given by jav(x, t) = −∂xu(x, t), Jav(x, t) = −∂xE(x, t),(4.12) and the dissipation field dav(x, t) is dav(x, t) = −νT, (4.13) with ν= (1 −α2)L2.(4.14) In order to write the average dissipation we have made use of the molecular chaos assumption or Stosszahlansatz. Specifically, we have assumed that hvl,pvl±1,pi=ul,pul±1,p +O(L−1).(4.15) In (4.14), we have introduced the macroscopic dissipation coefficient ν, which is the relevant parameter in the hydrodynamic space and time scales [145,146]. Note that if νis of order unity, the collisions have to be quasi-elastic 1 −α21. It is straightforward to combine (4.11), (4.12) and (4.13) to write closed equations for the hydrodynamic fields: average velocity and temperature, ∂tu(x, t) = ∂xxu(x, t),(4.16a) ∂tT(x, t) = −νT(x, t) + ∂xxT(x, t) + 2 [∂xu(x, t)]2.(4.16b) These equations must be solved submitted to suitable boundary conditions, which depend on the physical state under scrutiny. Note the resemblance between our hydrodynamic equations and those in (1.13). 4.1.4 Fluctuating hydrodynamics The balance equations (4.11) may also be written at the fluctuating level of description, by considering that v(x, t) and e(x, t) are fluctuating quantities, whose averages are u(x, t) and E(x, t). In this way, fluctuating balance equations are written for both v(x, t) and e(x, t), which are the continuum limit versions of the microscopic balance equations (4.4) and (4.5), namely ∂tv(x, t) = −∂xj(x, t), j(x, t) = −∂xv(x, t) + ξ(j)(x, t),(4.17a) ∂te(x, t) = −∂xJ(x, t) + d(x, t), J(x, t) = −∂xe(x, t) + ξ(J)(x, t).(4.17b) In the equations above, (j, J) are the fluctuating currents for momentum and energy, and (ξ(j), ξ(J)) are their corresponding noises. These noises have been shown to be Gaussian i i “thesis” — 2018/9/7 — 8:45 — page 78 — #86 i i i i i i 78 Chapter 4. Continuum limit, fluctuating hydrodynamics and finite size effects and white [145,146]. The amplitudes of their correlations hξ(γ)ξ(γ0)ican be cast in matrix form, hξ(γ)(x, t)ξ(γ0)(x0, t0)i=L−1Ξ(γγ0)δ(x−x0)δ(t−t0),(4.18) where (1,2) for γor γ0correspond to (j, J). These amplitudes have been computed within the local equilibrium approximation in [145], with the result Ξ= 2T(x, t)1 2u(x, t) 2u(x, t) 2[T(x, t)+2u2(x, t)].(4.19) The average velocity u(x, t) and the temperature T(x, t) must be calculated in the state corresponding to the physical situation of interest. Finally, the dissipation field d(x, t) is given by d(x, t) = −νθ(x, t) = −νe(x, t)−v2 R(x, t),(4.20) where v2 Ris the regular part of v2, defined as v2 R(x, t)≡v2(x, t)−L−1θ(x, t) lim ∆x→0δ(∆x).(4.21) This regular part of the velocity field has the property hv2 R(x, t)i=u2(x, t), as shown in appendix B. Equation (4.20) tells us that the fluctuations of the dissipation field are enslaved to those of the fluctuating temperature field θ(x, t). This is so because the dissipation noise ξ(d)is subdominant as compared to the current noises, since it scales as L−3instead of as L−1, as proven in [145]. 4.2 Dynamics of the one-particle distribution function Here, we apply the usual procedure of kinetic theory and map the master equation into a BBGKY hierarchy. In particular, we focus on the evolution equation for the one-particle distribution function at site land at time τ, which we denote by P1(v;l, τ). By definition, P1(v;l, τ) = ZdvPN(v, τ)δ(vl−v).(4.22) It is easy to show that none of the terms in the sum (4.7) contribute to the time evolution of P1except those corresponding to l−1 and l, because the collisions involving the pairs (l−1, l) and (l, l + 1) are the only ones which change the velocity at site l. Therefore, ∂τP1(v;l, τ) = ω× (Z+∞ −∞ dvl−1|∆l−1|β"P2(ˆ b−1 l−1{vl−1, v};l−1, l, τ) αβ+1 −P2(vl−1, v;l−1, l, τ)# +Z+∞ −∞ dvl+1|∆l|β"P2(ˆ b−1 l{v, vl+1};l, l + 1, τ) αβ+1 −P2(v, vl+1;l, l + 1, τ)#), (4.23) i i “thesis” — 2018/9/7 — 8:45 — page 79 — #87 i i i i i i Section 4.2. Dynamics of the one-particle distribution function 79 where P2(v, v0;l, l+1, τ) is the two-particle probability distribution for finding the particles at the l-th and (l+ 1)-th sites with velocities vand v0, respectively. For the special case β= 0, the evolution equation for P1can be further simplified, because the terms on the rhs of (4.23) coming from the loss (negative) terms of the master equation can be integrated. We get ∂τP1(v;l, τ) = ω−2P1(v;l, τ) + 1 αZ+∞ −∞ dvl−1P2(ˆ b−1 l−1{vl−1, v};l−1, l, τ) +1 αZ+∞ −∞ dvl+1P2(ˆ b−1 l{v, vl+1};l, l + 1, τ).(4.24) The equation for P1, either (4.23) for a generic βor (4.24) for β= 0, could be converted to a closed equation for P1by introducing the molecular chaos assumption, which in our present context means that P2(v, v0;l, l + 1, τ) = P1(v;l, τ)P1(v0;l+ 1, τ) + O(L−1).(4.25) By neglecting the O(L−1) terms in (4.25), we obtain a pseudo-Boltzmann or kinetic equation for P1, which determines the evolution of the one-time and one-particle averages under the assumption of O(L−1) correlations. Note that this “smallness” of two-particle correlations do not prevent them from being long-ranged. In the continuum (hydrodynamic-like) limit defined in (4.9), the evolution equation for the one-particle distribution function (4.23) can be written in a simpler form. The main idea is the quasi-elasticity of the microscopic dynamics, stemming from (4.14), i.e. α= 1 −L−2ν/2 + O(L−4). Then, P2ˆ b−1 l−1{vl−1, v};l−1, l, τ=P2v−ν 4L2∆l−1, vl−1+ν 4L2∆l−1;l−1, l, τ =h1 + ν 4L2∆l−1(∂vl−1−∂v)iP2(v, vl−1;l−1, l, τ) + O(L−4).(4.26) Moreover, we make use of the molecular chaos assumption (4.25) and identify P1(v;l, τ) = P1(v;x= (l−1)/L, t =ωτ/L2),(4.27) that is, we consider P1to be a smooth function of the hydrodynamic space and time variables xand t. Under the hypotheses outlined above, a lengthy but straightforward calculation gives for arbitrary βthat ∂tP1(v;x, t) = ∂xZ+∞ −∞ dv0|v0−v|β[P1(v0;x, t)∂xP1(v;x, t) −P1(v;x, t)∂xP1(v0;x, t)] −ν 2∂vZ+∞ −∞ dv0(v0−v)|v0−v|βP1(v0;x, t)P1(v;x, t).(4.28) The most important corrections to this equation emanate from finite size effects that give rise to nonzero correlations, i.e. the corrections to molecular chaos hypothesis that are i i “thesis” — 2018/9/7 — 8:45 — page 80 — #88 i i i i i i 80 Chapter 4. Continuum limit, fluctuating hydrodynamics and finite size effects expected to be of the order of L−1, see section 4.5 for details. This equation must be supplemented with suitable boundary conditions depending on the physical state under scrutiny. Note that the divergence structure of the rhs of (4.28) stems from the fact that P1(v;x, t) is a locally conserved quantity. If we took moments in (4.28), we would obtain the hydrodynamic equations for a generic value of β. It is clear, from the structure of this equation, that these hydrodynamic equations would be not closed and constitutive relations for the momentum and energy currents and the dissipation fields would be needed. We also recall that in our model there is no mass transport and therefore the “field density” n(x, t) is uniform and constant in time, n(x, t) = Z+∞ −∞ dv P1(v;x, t)=1,(4.29) which is consistent with the result ∂tn(x, t) = 0 obtained by marginalising vin (4.28). Again, for the case β= 0, the equation for P1in the continuum limit can be simplified, ∂tP1(v;x, t) = ∂xxP1(v;x, t) + ν 2∂v[(v−u(x, t))P1(v;x, t)] (4.30) The above equation is not linear for P1, since the average momentum is a functional thereof, u(x, t) = R+∞ −∞ dv vP1(v;x, t). Note that, consistently, by taking moments in (4.30), the average hydrodynamic equations (4.16) are reobtained. Moreover, the evolution equations for all the moments are closed (at any order) under the molecular chaos assumption, without further knowledge of the one-particle distribution function P1. Note that an appealing physical picture for P1(v;x, t) arises in the continuum limit. In fact, P1(v;x, t)dv dx = N X l=1 P1(v;l, t)dv ∆xΘ(L−1l−x)Θ(x+dx −L−1l) =L−1 N X l=1 P1(v;l, t)dv Θ(L−1l−x)Θ(x+dx −L−1l),(4.31) in which Θ(x) is Heaviside step function. The product of Heaviside functions selects the range of l’s corresponding to the interval (x, x +dx). Thus, P1(v;x, t)dv dx can be interpreted as the fraction of the total number of particles with velocities in the interval (v, v +dv) and positions in the interval (x, x +dx), which makes it neater the connection with the usual kinetic approach. 4.3 The one-particle distribution function for some physical states In this section, we restrict ourselves to the pseudo-Maxwell case β= 0. Therein, we analyze two physically relevant states that are typical of dissipative systems such as granular fluids. Specifically, we investigate the Homogeneous Cooling State (HCS) and the Uniform Shear Flow (USF) state. Although the model also admits stationary solutions i i “thesis” — 2018/9/7 — 8:45 — page 81 — #89 i i i i i i Section 4.3. The one-particle distribution function for some physical states 81 of Couette Flows for appropiate boundary conditions [145,146], they are not treated here. The theoretical results below are compared to numerical results, which have been obtained by means of Monte Carlo simulations described in appendix C. Specifically, we always plot the scaled one-particle distribution defined as ϕ(c;x, t) = pT(x, t)P1(v;x, t), c =v−u(x, t) pT(x, t),(4.32) in order to avoid visualizing the much sharper distributions P1that arise in some situations as a consequence of the cooling. 4.3.1 The Homogeneous Cooling State Let us consider a system with periodic boundary conditions. Moreover, the initial condition is “thermal”, that is, the random variables vl,0are Gaussian with zero average and unit variance. This means that we are choosing the uniform initial temperature to be the temperature unit. From this initial condition, the system evolves into the state known as the Homogeneous Cooling State (HCS), in which the system remains homogeneous and the temperature decays in time. At the (average) hydrodynamic level, one has [145,146] u(x, t)=0, THCS(x, t) = T(t= 0)e−νt.(4.33) Since the collision frequency of our model is velocity-independent, we are dealing with pseudo-Maxwell molecules and thus the temperature decays exponentially in time. This is the expected behavior, which replaces the typical algebraic decay for hard particles, known as Haff’s law [81]. As stated in 1.2.1, the HCS is known to be unstable: it breaks down in too large or too inelastic systems [62, 85]. In our model and in the continuum limit, this condition (studied in [79]) reads ν > νc, or, equivalently, L > Lc, where νc= 8π2, Lc= 2π√21−α2−1/2.(4.34) Note that the scaling of the critical length Lcis the same that that introduced in (1.12). When ν < νc(or L < Lc) there is no unstable mode. This instability mechanisms is completely analogous to the one found in granular gases for the shear mode [151]. Note that the amplification appears in the rescaled velocity ˜u(x, t) = u(x, t)/vth(t), being vth(t) = pTHCS(t), and not in the velocity u(x, t). Interestingly, the one-particle distribution function can be exactly calculated in the HCS. Taking into account the homogeneity of the state and the vanishing of the average velocity u(x, t), (4.30) for P1(v;t) simplifies to ∂tP1(v;t) = ν 2∂v[vP1(v;t)] .(4.35) This equation can be integrated right away to give P1(v;t) = eνt/2P1(veνt/2;t= 0).(4.36) i i “thesis” — 2018/9/7 — 8:45 — page 88 — #96 i i i i i i 88 Chapter 4. Continuum limit, fluctuating hydrodynamics and finite size effects which has to be integrated with the initial condition ˜ E(0) = 1. We have omitted the ν-dependence of ψHCS in order not to clutter our formulae. Therefore, up to order of L−1, we have ˜ E(t) = 1 + δ˜ E(t), δ ˜ E(t) = ψHCS ν t L.(4.61) which is expected to be valid as long as νψHCSt/L 1. Equation (4.61) can be also obtained by solving perturbatively the evolution equations for the temperature and twoparticle correlation function [146,150]. There is a critical dissipation value νψsuch that ψHCS vanishes, that is, νψ=νc/4=2π2, ψHCS(νψ)=0,(4.62) and the finite-size correction in (4.61) changes sign. Therefore, at this point we find a change in the time-derivative of δ˜ E(t). For large system sizes, the energy decays faster (slower) than the Haff’s law for ν < νψ(ν > νψ) because ψHCS <0 (ψHCS >0). Comparisons between these theoretical results and Monte Carlo simulations give an excellent agreement [146,150]. 4.5 Finite size effects: exact solution of the HCS In this section, we further analyze the velocity correlations. The average equation for the granular temperature (or the energy) in the HCS is closed only when the correlation hvlvl+1iis neglected, since it is expected to be of the order of L−1. In other words, the evolution equation for the temperature is closed in the molecular chaos approximation. Interestingly, for the case of Maxwell molecules we are considering in this work, we can account for the effect of the correlations in an exact way, thus going beyond molecular chaos. We assume that the system is in a spatial-translation-invariant state, such as the HCS. We define the set of spatial correlations of the velocity at time τas Ck(τ) = hvj(τ)vj+k(τ)i.(4.63) Here, krepresents the distance between the involved sites in the correlation. Note that the average temperature at any site jis given by C0, T(τ)≡C0(τ) = hv2 j(τ)i.(4.64) As a consequence of momentum conservation, in the center of mass frame we have the “sum rule” C0(τ)+2 L−1 2 X k=1 Ck(τ) = 0,∀τ, (4.65) where we have considered an odd L. i i “thesis” — 2018/9/7 — 8:45 — page 89 — #97 i i i i i i Section 4.5. Finite size effects: exact solution of the HCS 89 The evolution equation of the correlations is readily obtained from the master equation (4.7),2 ω−1∂τC0= (α2−1)(C0−C1),(4.66a) ω−1∂τC1=1−α2 2(C0−C1) + (1 + α)(C2−C1),(4.66b) ω−1∂τCk= (1 + α)(Ck+1 +Ck−1−2Ck),2≤k≤(L−1)/2,(4.66c) CL+1 2=CL−1 2,∀τ. (4.66d) Above, we have omitted the τ-dependence of the correlations to keep our notation simple, and written the evolution equations for odd L, because the “upper” boundary condition (for the maximum value of k) is simpler to write. For even L, the boundary condition would read CL 2+1 =CL 2−1. Of course, our choice of Las odd is irrelevant in the large system size limit L1. The hierarchy (4.66) can be exactly solved by reducing it to the eigenvalue problem of a certain matrix. We carry out this approach to the problem also for odd L. The problem for an even number particles may be solved by following an utterly similar strategy, but the boundary conditions are a little more involved to write. We do not present here these calculations because they do not provide any additional physical insight. First, it is useful to introduce a change of variables in order to make the matrix symmetric. Specifically, we define c0=C0, ck=√2Ck,1≤k≤(L−1)/2.(4.67) Second, we rewrite the hierarchy (4.66) as ω−1(1 + α)−1∂τc0=−(1 −α)c0+1−α √2c1,(4.68a) ω−1(1 + α)−1∂τc1=1−α √2c0−3−α 2c1+c2,(4.68b) ω−1(1 + α)−1∂τck=ck−1−2ck+ck+1,2≤k≤(L−3)/2,(4.68c) ω−1(1 + α)−1∂τcL−1 2=cL−3 2−cL−1 2,(4.68d) in which we have extracted the common factor (1 + α) on the rhs of (4.66) and made use of (4.66d) to write (4.68d) for cL−1 2. Now, we can solve the system above by a standard eigenvector method, that is, we seek solutions of the form ck=eλ(1+α)ωτ φk.(4.69) We denote the eigenvalues by λand its corresponding eigenvector by φ,φkis thus the 2Equation (4.66b) had a typo in our original publication [150]. Specifically, the prefactor of (C0−C1) was (1 −α2) instead of (1 −α2)/2, as pointed out by the authors of [160]. i i “thesis” — 2018/9/7 — 8:45 — page 90 — #98 i i i i i i 90 Chapter 4. Continuum limit, fluctuating hydrodynamics and finite size effects k-th component thereof. In this way, we reach the system λφ0=−(1 −α)φ0+(1 −α) √2φ1,(4.70a) λφ1=(1 −α) √2φ0−(3 −α) 2φ1+φ2,(4.70b) λφk=φk−1−2φk+φk+1,2≤k≤(L−3)/2,(4.70c) λφL−1 2=φL−3 2−φL−1 2.(4.70d) Equations (4.70) are a system of second-order difference equations for φkwith contant coefficients, in which (4.70c) is the general equation and (4.70b) and (4.70d) are their boundary conditions. On top of that, (4.70a) acts as an extra condition that ensures momentum conservation, as shown below (see also [146,150]). The general solution of (4.70c) is of the form φk>0=rk[117], which substituted into (4.70c) has two solutions (r1, r2) that verify r1r2= 1,(4.71a) r1+r2= 2 + λ. (4.71b) We introduce a new variable q∈[0, π] such that r1=eiq and r2=e−iq, as suggested by (4.71a). Note that |r1|=|r2|= 1, if one of the roots were larger than one it would lead to correlations increasing with k, which is physically absurd. Moreover, from a purely mathematical point of view, restricting ourselves to |r1|=|r2|= 1 leads to a complete set of eigenvectors. From (4.71b), we obtain λ(q) = 2(cos q−1),(4.72) and the corresponding eigenvector is given by φk>0(q) = A eikq +B e−ikq,(4.73a) φ0(q) = 1−α √2 (2 cos q−1−α)(A eiq +B e−iq).(4.73b) The boundary conditions (4.70b) and (4.70d) determine the constants Aand B, and also the allowed values of the “index” q. The determinant of the linear system for Aand Bmust be zero, which is equivalent to impose that qmust be a zero of the function g(q) = 2 sin L+ 3 2q−(5 + 3α) sin L+ 1 2q+ (5 + 7α) sin L−1 2q −(3 + 5α) sin L−3 2q+ (1 + α) sin L−5 2q.(4.74) This function has (L+1)/2 different zeros in the half-open interval [0, π), which we denote by qn:q0= 0, qnis the n-th nonvanishing zero of g(q), n= 1,...,(L−1)/2. In addition q=πalso makes g(q) vanish, but it does not correspond to an eigenvalue because the associated eigenvector would be identically zero, as shown below. In this way, we find (L+ 1)/2 eigenvalues λn= 2(cos qn−1),(4.75) i i “thesis” — 2018/9/7 — 8:45 — page 91 — #99 i i i i i i Section 4.5. Finite size effects: exact solution of the HCS 91 Figure 4.4: Plot of the function g(q) defined in (4.74). The system size is L= 11. Its zeros qndetermine the eigenvalues, as given by (4.75). The first zero is always q0= 0, and there are (L−1)/2 additional zeros qi6= 0, i= 1,...,(L−1)/2. The inset shows a zoom of the small qregion, to facilitate the identification of q1and q2. the corresponding eigenvectors of which give a complete set for our problem. In figure 4.4, we plot the function g(q) for L= 11, which has six zeros in the interval [0, π). The proportionality relation obtained between Aand Bmakes it possible to write the eigenvector corresponding to the eigenvalue λnup to a normalization constant Nn, φ(n) k>0=Nncos L 2−kqn,(4.76a) φ(n) 0=(1 −α)Nn √2 (2 cos qn−1−α)cos L 2−1qn.(4.76b) The above expressions clearly show that there is no eigenvector for q=π, since all its components are zero (recall that Lis odd). The constant Nnis chosen to obtain a orthonormal set of eigenvectors, in the sense that L−1 2 X k=0 φ(n) kφ(n0) k=δnn0.(4.77) We do not give the explicit expression for Nnbecause it is quite involved and unnecessary for our purposes. The eigenvector corresponding to q0= 0 is particularly simple, (4.76) implies that φ(0) 0=N0 √2, φ(0) k>0=N0,N0=r2 L.(4.78) i i “thesis” — 2018/9/7 — 8:45 — page 92 — #100 i i i i i i 92 Chapter 4. Continuum limit, fluctuating hydrodynamics and finite size effects Then, the orthogonality relation of φ(0) and φ(n)(n6= 0) makes it possible to write a “sum rule” for the components of the latter eigenvectors, φ(n) 0+√2 L−1 2 X k=1 φ(n) k= 0, n > 0.(4.79) This sum rule is connected with (4.65), which stemmed from momentum conservation. It also allows us to write φ(n) 0in a more convenient form for some calculations, φ(n) 0=−Nn √2csc qn 2sin L−1 2qn,(4.80) which does not depend explicitly on α. Finally, we have all the ingredients to build the general solution of (4.68) as the sum ck= L−1 2 X n=1 aneλn(1+α)ωτ φ(n) k,(4.81) where anis given in terms of the initial conditions by an= L−1 2 X k=0 φ(n) kck(0).(4.82) The sum in (4.81) starts from n= 1 because a0= 0, since a0=N0  c0(0) √2+ L−1 2 X k=1 ck(0) =N0 √2 C0(0) + 2 L−1 2 X k=1 Ck(0) = 0.(4.83) We have made use of momentum conservation, as expressed by the sum rule (4.65), to obtain the last equality. 4.5.1 Eigenvalues for large systems Here, we would like to derive an approximate expression for the eigenvalue spectrum in the large system size limit L1. Therefore, we consider that the microscopic dynamics is quasi-elastic by introducing the macroscopic dissipation coefficient ν, (1 −α2)L2=ν, as in (4.13). The eigenvalues are given by the zeros of the function g(q) in (4.74), and we expand this function for q1 by introducing the scaling Q=qL, with the result tan Q 2ν 2Q2L−2−Q4L−4+1 2Q5L−5= 0.(4.84) We are assuming that Qis of the order of unity and have neglected terms of the order of L−6. i i “thesis” — 2018/9/7 — 8:45 — page 93 — #101 i i i i i i Section 4.5. Finite size effects: exact solution of the HCS 93 Figure 4.5: Plot of the difference Q1−Q(0) 1as a function of the system size L. Specifically, we present the plots for ν=π2,ν=νψ= 2π2and ν= 4π2. For each of these values, there are two curves: the theoretical curve Q(1) 1/L (lines) and the numerical estimate of Q1(symbols). See the legend for details. The finite-size correction is especially small for ν=νψ. In order to obtain an analytical approximation for the eigenvalues, we propose an expansion of Qn=qnLin powers of L−1,Qn=Q(0) n+Q(1) nL−1+O(L−2). To the lowest order, we obtain Q(0) 1=rν 2,(4.85a) Q(0) n= 2(n−1)π, n = 2,...,(L−1)/2.(4.85b) Moreover, the finite size corrections are Q(1) 1=ν 8 tan 1 2pν 2,(4.86a) Q(1) n=16(n−1)3π3 8(n−1)2π2−ν, n = 2,...,(L−1)/2.(4.86b) Note that Q(1) 1vanishes at ν=νψ= 2π2whereas it diverges at ν=νc= 8π2. The former property is connected to the change of sign in the finite-size correction to the cooling rate of the HCS, whereas the latter gives rise to the instability of the HCS, as discussed in sections 4.4. In figure 4.5, we check the above expansion for the zeros of the function g(q). Specifically, we do so for the first zero q1: the numerical estimation of q1is compared with the expansion in (4.85) and (4.86) by plotting Q1−Q(0) 1as a function of the system size L. It is observed that this difference tends to zero as the system size increases, for all the i i “thesis” — 2018/9/7 — 8:45 — page 94 — #102 i i i i i i 94 Chapter 4. Continuum limit, fluctuating hydrodynamics and finite size effects considered values of ν. The finite size correction is especially small for ν=νψ= 2π2, for which the theoretical correction Q(1) 1vanishes. Therefore, finite size corrections are as small as possible for this case, which makes it particularly adequate to investigate the hydrodynamic description, as done in section 4.3 and also in [79,145,146]. We want to emphasize that the instability of the HCS is reobtained here as a crossing between the first two nonzero eigenvalues: for ν=νc= 8π2, we have that Q(0) 1=Q(0) 2. On the one hand, for ν < νc, the largest nonvanishing eigenvalue is λ1(λ1<0) and dominates the long-time dynamics of the system: the energy C0and all the correlations Ckdecay with exp[λ1ω(1 + α)τ] = exp(νr HCSt), see below. On the other hand, for ν > νc, the dominant term is the one corresponding to Q2≃2πand the long time behavior of the system becomes independent of ν. The large system size limit of the eigenvalues is then λn=−Q(0) n 2 L2"1 + L−12Q(1) n Q(0) n +O(L−2)#.(4.87) Moreover, the exponent in (4.81) controlling the time dependence of the contribution for each mode is λn(1 + α)ωτ ∼ −2Q(0) n 2"1 + L−12Q(1) n Q(0) n +O(L−2)#t, (4.88) which shows the relevance of the hydrodynamic scale t, defined in (4.9), in the large system size limit. 4.5.2 Long time limit Equation (4.81) gives the general time evolution for the velocity correlations. Here, we show that these correlations tend to their HCS values in the long time limit, provided that ν < νc, that is, we are below the instability. Let us consider the scaled correlations ˜ Ck ˜ Ck(τ) = Ck(τ) C0(τ)=ck(τ) √2c0(τ),(4.89) that is,, we scale the correlations with the energy C06= 0. For long enough times, the only relevant contribution to (4.81) stems from the maximum (minimum in absolute value) eigenvalue λ1. Thus, the time dependence for all the correlations Ck(or ck) are the same and, consequently, the quotient in (4.89) becomes time-independent for long enough times. Making use of (4.80) and (4.81), ˜ Ck=φ(1) k √2φ(1) 0 =−sin q1 2csc L−1 2q1cos L 2−kq1,(4.90) which is nothing but the discrete version of the continuum solution obtained in [146,150]. i i “thesis” — 2018/9/7 — 8:45 — page 95 — #103 i i i i i i Section 4.5. Finite size effects: exact solution of the HCS 95 We can also derive the rate at which the energy and all the correlations are decaying in the long time limit. Particularizing (4.88) for n= 1, we have that λ1(1 + α)ωτ ∼ −νt 1−L−1ψHCS +O(L−2)=−νr HCSt, (4.91) where νr HCS is the “renormalized” by fluctuations cooling rate introduced in [146, 150], by applying a multiple scale analysis to find finite size corrections to the hydrodynamic description. Thus, the energy is given by C0(t) = T(t= 0) exp(−νr HCSt) (4.92) and the correlations Ck, with k > 0, follow from (4.90). i i “thesis” — 2018/9/7 — 8:45 — page 96 — #104 i i i i i i 96 Chapter 4. Continuum limit, fluctuating hydrodynamics and finite size effects i i “thesis” — 2018/9/7 — 8:45 — page 97 — #105 i i i i i i Chapter 5 Global stability, derivation of the H-theorem Here, we aim at investigating the global stability and the possibly associated H-theorem for a generalization of the model introduced in the previous chapter. At difference with the approach in [96,97], our analysis is not restricted to spatially homogeneous situations: we consider the whole space and velocity dependence of the one-particle probability distribution function. More specifically, a general energy injection mechanism is introduced, in which the system may be driven through both the boundaries and the bulk. We show that, under quite general conditions, the steady state is globally stable. Independently of the initial preparation, the system always ends up in the steady state. Interestingly, it is not necessary to have an H-theorem to prove this: it suffices to show that His decreasing in the long-time limit, not for all times. In this regard, we find a situation that is similar to the proof of the tendency towards the equilibrium curve in systems whose dynamics is governed by master equations with time-dependent transition rates [161–167]. Our proof of global stability also makes it possible to show the inadequacy of Boltzmanns HB, defined in (1.16), as a candidate for a Lyapunov functional in inelastic systems. Not only is this done for the simplified models considered throughout this thesis, but for a general collision term that does not conserve energy in collisions. Therefore, this result also applies to the inelastic Boltzmann or Enskog equations, employed for granular fluids. Specifically, the main idea is the possibility of reversing the sign of dHB/dt by a suitable choice of the initial PDF. Thus, dHB/dt cannot have a definite sign. In this sense, our result can be understood as a generalization of that in [99], within the first Sonine approximation of the inelastic Boltzmann equation, to an arbitrary collision kernel with nonconservative interactions. Having proved global stability by showing that His a nonincreasing functional for long times, a natural question arises. Is Ha Lyapunov function, that is, a nonincreasing functional for all times? There does not seem to be a unique proof, valid for any driving mechanism, even within the framework of our simplified model. Notwithstanding, we 97 i i “thesis” — 2018/9/7 — 8:45 — page 104 — #112 i i i i i i 104 Chapter 5. Global stability, derivation of the H-theorem Figure 5.1: Time evolution of the functional HB. Specifically, we consider two different initial conditions in a uniformly heated system. Both simulations start from a Gaussian distribution with a homogeneous temperature, which slightly differs from the stationary one: (i) T(t= 0) = 1.1Ts(blue circles) and (ii) T(t=0)=0.9Ts(red triangles). Both functionals are symmetric with respect to its stationary value in agreement with the prediction of the linear approximation. Consistently, the mean value of both curves (dashed line) remains approximately constant throughout. We have considered a system with parameters ν= 20 and ξ= 50, size N= 330, and averaged over 3000 trajectories. and, in general, ln P(s) 1is not a sum of constants of motion. In fact, again the sign of dHB/dt|lin is reversed when ∆P1→ −∆P1, similarly to the homogeneous case. In figure 5.1, we have numerically checked this prediction for the sheared system, with the resulting evolution of HBbeing completely similar to that for the uniformly heated case in figure 5.1. 5.4 H-theorem for some specific NESS Here we prove that the functional H[f] is monotonically decreasing for all times in some specific physical situations. Our proof applies both to the sheared and the uniformly heated systems described in section 5.1.2. To be as general as possible, we consider a system that is both heated and sheared: a6= 0 and ξ6= 0. In this situation, the boundary conditions for the PDF are given by (4.43), which lead to the Lees-Edwards boundary conditions (4.40) for the averages u(x, t) and T(x, t). The steady solution of the hydrodynamic equations is us(x) = ax−1 2, Ts=2a2+ξ ν.(5.29) On the one hand, the average velocity has a linear profile, similarly to the situation in i i “thesis” — 2018/9/7 — 8:45 — page 105 — #113 i i i i i i Section 5.4. H-theorem for some specific NESS 105 Figure 5.2: Time evolution of the functional HB. Analogously to figure 5.1, we consider two different initial conditions, but now for a sheared system. Both simulations start from a Gaussian distribution with the stationary average velocity profile but a homogeneous temperature, which slightly differs from the stationary one: (i) T(t=0)=1.1Ts(blue circles) and (ii) T(t= 0) = 0.9Ts(red triangles). We have considered a system with parameters ν= 20 and a= 5, size N= 660, and averaged over 6000 trajectories. the USF state. On the other hand, the temperature remains homogeneous but its steady value has two contributions, one coming from the shear and the other from the stochastic thermostat. The viscous heating 2(∂xu)2and uniform heating ξterms cancel the cooling term −νT for all x. The stationary solution of the kinetic equation is the Gaussian distribution (4.45) corresponding to the hydrodynamic fields in (5.29). Of course, the USF state and NESS of the uniformly heated system in section 5.1.2 can be easily recovered as particular cases of (4.45): for (a6= 0, ξ = 0) and (a= 0, ξ > 0), respectively. Then, we turn now to the question of the existence of an H-theorem, that is, the existence of a nonequilibrium entropy ensuring the monotonic approach of the one-particle PDF to the steady state. Our starting point is the following expansion of the one-particle PDF in Hermite polynomials, P1(v;x, t) = 1 p2πT(x, t)exp −[v−u(x, t)]2 2T(x, t)(1 + ∞ X n=3 γn(x, t)Hn "v−u(x, t) pT(x, t)#), (5.30) which is known as the Gram-Charlier series [173–176]. Therein, u(x, t) and T(x, t) are the (exact) average velocity and temperature stemming from the hydrodynamic equations for the considered distribution. Note that the argument of the Hermite polynomials is the scaled velocity cdefined in (4.32). The above expansion is suggested by the Gaussian shape of the stationary PDF in (4.45). Making use of the definition (4.32) of the scaled i i “thesis” — 2018/9/7 — 8:45 — page 106 — #114 i i i i i i 106 Chapter 5. Global stability, derivation of the H-theorem one-particle distribution function ϕ(c;x, t) and the orthogonality relation of the Hermite polynomials [177], it is readily obtained that ϕ(c;x, t) = 1 √2πexp −c2 2"1 + ∞ X n=3 γn(x, t)Hn(c)#,(5.31) with γn(x, t) = 1 n!Zdc Hn(c)ϕ(c;x, t).(5.32) Also, we could have written γnas a combination of moments of the distribution. Some comments on the Gram-Charlier expansion are pertinent. First, note that n≥3 in the sum: γ1=γ2= 0 because the zero-th order Gaussian contribution exactly gives the first two moments hvi(x, t) = u(x, t) and hv2i(x, t) = u2(x, t) + T(x, t). Second, if P1(v;x, t) were symmetric with respect to v=u, that is, h(v−u)2n+1i= 0 for all n∈N, only even values of nwould be present in the sum and one recovers the usual expansion in Sonine-Laguerre polynomials of kinetic theory. Finally, it is worth stressing that the series (5.30) converges for functions such that the tails of ϕ(c;x, t) approach zero faster than e−c2/4for c→ ±∞ [176,178,179]. Now we substitute the Gaussian stationary solution (4.45) and the Gram-Charlier series (5.30) into the three contributions to dH/dt, given by (5.16), (5.17) and (5.18). For the inelastic term, it is readily obtained that dH dt inel =ν 2TsZdx (u−us)2.(5.33) For the diffusive and noise terms, the key ideas are a change of the integration variable from vto c= (v−u)/√Tand the use of the recursion relations and the orthogonality property of the Hermite polynomials [177]. Working along these guidelines, we arrive at dH dt diff =−Zdx T u0 T−u0 s Ts2 −u02 s T2 sZdx (u−us)2 −1 √2πZdx dc e−c2/2 1 + P∞ n=3 γnHn(c)(T0 2TH2(c) + ∞ X n=3 γ0 nHn(c) − ∞ X n=3 γnu0 √TnHn−1(c) + ∞ X n=3 γnT0 2T[Hn+2(c) + nHn(c)])2 ,(5.34) dH dt noise =−ξ 2Zdx T 1 T−1 Ts2 −ξ 2T2 sZdx (u−us)2 −ξ 2√2πZdx dc e−c2/2 1 + P∞ n=3 γnHn(c) 1 T"∞ X n=3 γnnHn−1(c)#2 ,(5.35) where the prime stands for spatial derivative. In (5.33), (5.34) and (5.35) there are several terms multiplying Rdx (u−us)2: they cancel out when we take into account the i i “thesis” — 2018/9/7 — 8:45 — page 107 — #115 i i i i i i Section 5.4. H-theorem for some specific NESS 107 equation for the (spatially homogeneous) stationary temperature. Therefore, the sum of the remaining terms leads right to dH dt =A(t) + B(t),with both A(t), B(t)≤0,(5.36) being A(t) = −Zdx T "u0 T−u0 s Ts2 +ξ 21 T−1 Ts2#,(5.37) and B(t) the sum of the second and third line in (5.34) with the second line in (5.35), that is, B(t) = −1 √2πZdx dc e−c2/2 1 + P∞ n=3 γnHn(c)   ξ 2T"∞ X n=3 γnnHn−1(c)#2 +T0 2TH2(c) + ∞ X n=3 γ0 nHn(c)− ∞ X n=3 γnu0 √TnHn−1(c) + ∞ X n=3 γnT0 2T{Hn+2(c) + nHn(c)}#2   , (5.38) In conclusion, dH/dt ≤0 for all times and we have shown that the H-theorem holds for the sheared and heated system. Rigorously, our proof holds for those PDFs such that the above Hermite expansion converges. Note that the proof remains valid for the approach to any NESS, whose PDF is a Gaussian with a homogeneous temperature, independently of the corresponding boundary conditions. In section 5.2, we have already demonstrated that dH/dt only vanishes for P1(v;x, ∞) = P(s) 1(v;x), but the same result can be rederived here in a different way. By imposing that both A(t) and B(t) vanish in the long time limit and making use of the hydrodynamic equations for the averages, it can be shown that u(x, ∞) = us(x), T(x, ∞) = Tsand γn(x, ∞) = 0, ∀n≥3. 5.4.1 USF state: simulations Here we consider the sheared system in order to check numerically our theoretical predictions. Throughout this section, we employ the values of the parameters ν= 20 and a= 5. Note that ξ= 0, since there is no stochastic forcing. First, in figure 5.3, we show the evolution of the distribution and the H-functional from a Gaussian initial condition with the steady velocity profile u(x, 0) = us(x) but a higher temperature, T(t= 0) = 7 Ts. In the top panel, we depict the velocity distribution at x= 1/4 for several times. All of them are Gaussian, which agrees with the theoretical prediction of the kinetic equation: when the initial velocity profile coincides with the steady one and only the temperature is perturbed, an initially Gaussian PDF remains Gaussian for all times. Indeed, we can see in the inset how the excess kurtosis, κ=h[v−u(x)]4i/h[v−u(x)]2i2−3,(5.39) only fluctuates around zero at the considered position, consistently with the Gaussian shape. In the bottom panel, it is neatly observed that the H-functional monotonically decreases in time. i i “thesis” — 2018/9/7 — 8:45 — page 108 — #116 i i i i i i 108 Chapter 5. Global stability, derivation of the H-theorem Figure 5.3: Relaxation towards the USF state. The initial distribution is Gaussian, with average velocity in us(x) and temperature T= 7 Ts. (Top) Velocity distribution function at x= 1/4. Four different times are shown, as labeled. In the inset, we present the evolution of the excess kurtosis. Solid lines correspond to the (theoretical) Gaussian distributions for the plotted times, except for the longest in which it represents the theoretical steady distribution. (Bottom) Relaxation of the Hfunctional, which is clearly monotonically decreasing to zero. System size is N= 660, parameters are ν= 20 and a= 5, and curves are averaged over 6000 runs. i i “thesis” — 2018/9/7 — 8:45 — page 109 — #117 i i i i i i Section 5.4. H-theorem for some specific NESS 109 Second, we look into the relaxation to the USF state from another initial preparation, for which the velocity profile u(x, 0) is different from the stationary but T(x, 0) = Ts. The numerical results are shown in figure 5.4, and for the sake of simplicity we use again an initial Gaussian distribution. Specifically, we have that u(x, 0) = us(x)+4.4 sin(2πx). Here, the departure from the Gaussian shape is evident, and thus we have not plotted the kurtosis. Consistently with our theoretical prediction, we get again a monotonous relaxation of Htowards its null stationary value. Finally, we consider situations for which the above presented proof is not rigorously applicable. As stated before, the Gram-Charlier series does not converge when the tails of the distribution decay to zero slower than the square root of the Gaussian. Nevertheless, when all the coefficients γndefined in (5.32) exist and are finite, we still expect the H- theorem to hold. We illustrate this situation with an initial exponential distribution; specifically, we consider P1(v;x, 0) = 1 p2T(t= 0) exp "√2|v−u(x, t = 0)| pT(t= 0) #,(5.40) with u(x, t = 0) = us(x) + 4.4 sin(2πx) and T(t=0)=0.1Ts. In agreement with our expectation, the H-functional also monotonically decreases in figure 5.5. 5.4.2 Numerical results in the uniformly heated system To conclude, we put forward the results of simulations for the uniformly heated system. Specifically, our simulations have been done for ν= 20, a= 0 (no shear) and ξ= 50. In order not to be repetitive, we only present the more complex case in figure 5.6: the relaxation towards the steady state from an initial exponential distribution, as given by (5.40). In particular, we consider that u(x, t =0)=4.4 sin(2πx) and T(t=0)=0.1Ts. Note that the perturbation from the steady values is the same as in figure 5.5 for the sheared case. Once more, we observe the monotonic relaxation of Htowards its stationary value, in neat agreement with our theoretical result, even for a initial distribution for which the Gram-Charlier series does not converge. i i “thesis” — 2018/9/7 — 8:45 — page 110 — #118 i i i i i i 110 Chapter 5. Global stability, derivation of the H-theorem Figure 5.4: The same plots as in figure 5.3, but starting from a different initial condition. Now, the initial PDF is a Gaussian with average velocity u(x, 0) = us(x)+4.4 sin(2πx) and temperature T(t= 0) = Ts. In the top panel, solid lines correspond to the theoretical PDFs for the initial time and the steady state. In the bottom panel, Hdecreases again monotonically towards its steady value, consistently with our theory. i i “thesis” — 2018/9/7 — 8:45 — page 111 — #119 i i i i i i Section 5.4. H-theorem for some specific NESS 111 Figure 5.5: The same plots as in figure 5.3, but starting from an initial PDF with a divergent Gram- Charlier series. Specifically, they correspond to an exponential initial distribution with average velocity u(x, t = 0) = us(x)+4.4 sin(2πx) and temperature T(t=0)=0.1Ts. i i “thesis” — 2018/9/7 — 8:45 — page 112 — #120 i i i i i i 112 Chapter 5. Global stability, derivation of the H-theorem Figure 5.6: Numerical results for the uniformly heated system. Similar to what we did in figure 5.3, we show the time evolution of the PDF (top) and the time evolution of the H-functional (bottom). The system is initially prepared with an exponential PDF with average velocity u(x, t =0)=4.4 sin(2πx) and temperature T(t=0)=0.1Ts. System size is N= 330 and the curves have been averaged over 3000 runs. i i “thesis” — 2018/9/7 — 8:45 — page 113 — #121 i i i i i i Chapter 6 Memory effects in athermal systems The objectives of this chapter are two-fold. First, we present a general, rigorous, derivation of the linear response expression for the Kovacs hump, which was described in the introduction, specifically in section 1.2.3. This is done for systems with a realistic continuous time dynamics, at both the mesoscopic and macroscopic levels of description. In the former, our starting point is the master equation for the probability distribution, which is always linear. In the latter, we begin by considering the hierarchy of equations for the moments, which in general is nonlinear. The presented proof is valid for both molecular and athermal systems, since no hypothesis is needed with regard to the form of either the stationary probability distribution or the fluctuation-dissipation relation. The organization of this chapter is as follows. To start with, we focus on the linear response regime in both the mesoscopic and macroscopic levels of description in sections 6.1.1 and 6.1.2, respectively. Afterwards, section 6.2 is devoted to the particularization of this linear theory to an even more generalized version of our simple model of a granular gas. Here, we consider a general collision rate that makes it necessary to introduce the typical Sonine expansion of kinetic theory in order to close the evolution equation of the temperature. Finally, we test our theoretical predictions by comparison with Monte Carlo simulations in section 6.3. 6.1 Linear theory for the Kovacs memory effect 6.1.1 Master equation approach First, we consider a general system, the state of which is completely characterized by a vector xwith Mcomponents, x={x1, x2, . . . , xM}. For example, in a one-dimensional Ising chain of Nspins, xi=σi=±1, and M=N; for a gas comprising Nparticles with positions riand velocities vi,M= 6N, and x={r1,v1,...,rN,vN}. Henceforth, for the sake of simplicity, a notation suitable for systems in which the states can be labeled with a discrete index α, 1 ≤α≤Ω, is employed. For example, this is the case of the Ising system above, where Ω = 2N. The generalization for a continuous index is straightforward, by 113 i i “thesis” — 2018/9/7 — 8:45 — page 120 — #128 i i i i i i 120 Chapter 6. Memory effects in athermal systems In the nonlinear first Sonine approximation, the evolution equation for the temperature is readily obtained d dtT=−ζ0T1+ β 21 + β(2 + β) 16 a2+β(2 + β)(2 −β)(4 −β) 1024 a2 2+ξ, (6.33a) where ζ0=π−1 /221+βΓ3+β 2. Unless β= 0 (the Maxwell molecules we have considered in previous chapters), the equation for the temperature is not closed. Then, we write down the equation for a2: again, after a lengthy, but straightforward calculation, we derive Td dta2=−2ξa2−ζ0 3β T1+ β 21 + 56 + β(6 + β) 16 a2 −(2 + β)[384 + (2 −β)β(4 + β)] 1024 a2 2−3(4 −β)(2 −β)(2 + β) 512 a3 2.(6.33b) The evolution equations in the standard first Sonine approximation are obtained just neglecting nonlinear terms in a2in (6.33), that is, d dtT=−ζ0T1+ β 21 + β(2 + β) 16 a2+ξ, (6.34a) Td dta2=−ζ0 3β T1+ β 21 + 56 + β(6 + β) 16 a2−2ξa2.(6.34b) For ξ6= 0, the steady solution of these equations is Ts=  ξ ζ0h1 + β(2+β) 16 as 2i  2 2+β , as 2=−16β 96 + 56β+ 6β2+β3.(6.35) Note that (i) 0 ≤ |as 2| ≤ 0.133 for 0 ≤β≤2, which makes it reasonable to use the first Sonine approximation, and (ii) as 2is independent of the driving intensity ξ. This will be useful in the linear response analysis, to be developed below. A sudden change in the driving only changes the stationary value of the temperature, but not that of the excess kurtosis. If ξ= 0, the system evolves towards the homogeneous cooling state, in which the excess kurtosis tends to the value aHCS 2=−16 56 + β(6 + β),(6.36) as predicted by (6.34b), and the temperature decays following Haff’s law, dT/dt ∝ −T1+ β 2. From now on, we use reduced temperature and time, θ=T Ts , s =ζ0Tβ/2 st. (6.37) i i “thesis” — 2018/9/7 — 8:45 — page 121 — #129 i i i i i i Section 6.2. Kovacs hump in the model of granular media 121 The steady temperature Tsplays the role of a natural energy (or granular temperature) unit. In reduced variables, the evolution equations are d dsθ= 1 −θ1+ β 2+β(2 + β) 16 as 2−a2θ1+ β 2,(6.38a) θd dsa2=κ1a2−aHCS 21−θ1+ β 2−κ2(a2−as 2),(6.38b) where we have introduced two parameters of the order of unity, κ1=−β 3aHCS 2 , κ2=−β 3as 2 ,(6.39) 0≤κ1≤3 and 2 ≤κ2≤5 for 0 ≤β≤2. The evolution equations in the first Sonine approximation, (6.34) or (6.38), are the particularization of the equations for the moments (6.20) to our model: J= 2, and z1=T (or θ), z2=a2. Consistently, they are nonlinear, although here, due to the simplifications introduced in the first Sonine approximation, only nonlinear in θ. When the system is close to the NESS, (6.38) can be linearized around it by writing θ= 1+∆θ,a2=as 2+∆a2, d ds ∆θ ∆a2=M·∆θ ∆a2,M= −2(2+β)(12+β) 48+4β+β2−β(2+β) 16 −κ11 + β 2as 2−aHCS 2−κ2!.(6.40) Of course, the general solution of this linear system for arbitrary initial conditions ∆θ(0) and ∆a2(0) can be immediately written, but we omit it here. 6.2.2 Kovacs Hump in Linear Response Now, we look into the Kovacs hump in the linear response approximation. Following the discussion leading to (6.18), first we have to calculate the relaxation function φTfor the granular temperature. The system is at the steady state corresponding to a driving ξ0 for t < 0; at t= 0, the driving is instantaneously changed to ξ, and only the linear terms in ∆ξ=ξ−ξ0are retained. We choose the normalization of φT(s) in such a way that φT(0) = 1, that is, φT(s)≡lim ∆T(0)→0 ∆T(s) ∆T(0) = lim ∆θ(0)→0 ∆θ(s) ∆θ(0).(6.41) Since Tschanges with ξ, but a2does not, we have to solve (6.40) for ∆a2(0) = 0 and arbitrary (small enough) ∆θ(0). The solution is φT(s) = c+eλ+s+c−eλ−s,(6.42a) c+=M11 −λ− λ+−λ− , c−=λ+−M11 λ+−λ− ,(6.42b) where Mij is the (i, j) element of the matrix Mand λ±its eigenvalues, λ±=Tr M±p(Tr M)2−4 det M 2=Tr M±p(M11 −M22)2+ 4M12M21 2.(6.43) i i “thesis” — 2018/9/7 — 8:45 — page 122 — #130 i i i i i i 122 Chapter 6. Memory effects in athermal systems Both eigenvalues λ±are negative, since Tr M<0 and det M>0 for all β > 0. Therefore, |λ+|<|λ−|, and it is λ+that dominates the relaxation of the granular temperature for long times. Moreover, c±>0, and thus, the linear relaxation function φT(s) is always positive and decays monotonically to zero. Next, we consider a Kovacs-like experiment: the system was at the NESS corresponding to a driving ξ0, with granular temperature Ts,0for t < 0; the driving is suddenly changed to ξ1at t= 0 so that the system starts to relax towards a new steady temperature Ts,1for 0 ≤t≤tw, and this relaxation is interrupted at t=tw, because the driving is again suddenly changed to the value ξsuch that the stationary granular temperature Tsequals its instantaneous value at tw. The time evolution of the granular temperature for t≥twis given by the particularization of (6.18) and (6.19) to our situation, that is, KT(s) = ξ0−ξ1 ξ0−ξφT(s)−ξ−ξ1 ξ0−ξφT(s−sw),ξ−ξ1 ξ0−ξ1 =φT(sw),(6.44) where we have made use of the normalization φT(0) = 1. In the linear response approximation, the jumps in the driving values can be substituted by the corresponding jumps in the stationary values of the granular temperature. The structure of the linear relaxation function φT(s), as a linear combination of decreasing exponentials exp(λ±t), λ±<0, with positive weights c±, assures that the Kovacs behavior is normal: (i) KT(s) is always positive and bounded from above by φT(s) and (ii) there is only one maximum at a certain time sk> sw[107]. The anomalous behavior found in the uniformly heated hard-sphere granular for large enough inelasticity is thus not present here. This is consistent with the quasi-elastic limit we have introduced to simplify the collision operator. 6.2.3 Nonlinear Kovacs Hump Here, we consider the Kovacs hump for arbitrary large driving jumps. In our model, we can make use of the smallness of a2, which is assumed in the first Sonine approximation, in order to introduce a perturbative expansion of (6.38) in powers of as 2. The procedure is completely analogous to that performed in [110,111] for a dilute gas of inelastic hard spheres, and thus, we omit the details here. We start by writing a2=as 2A2, with A2of the order of unity, and θ(s) = θ0(s) + as 2θ1(s) + . . . , A2(s) = A20(s) + as 2A21(s) + . . . . (6.45) These expansions are inserted into (6.38), which have to be solved with initial conditions θ(sw) = 1, A2(sw) = Aini 2. To the lowest order, θ0(s) = 1, whereas A20(s) decays exponentially to one, A20(s)−1∼Aini 2−1e−κ2(s−sw).(6.46) In order to describe the Kovacs hump, we compute θ1(s) that verifies the evolution equation dθ1 ds =−1 + β 2θ1+β(2 + β) 16 (A20 −1) ,(6.47) which gives θ(s)−1∼aini 2−as 2β(2 + β) 8(2 + β−2κ2)he−κ2(s−sw)−e−(1+ β 2)(s−sw)i, s ≥sw.(6.48) i i “thesis” — 2018/9/7 — 8:45 — page 123 — #131 i i i i i i Section 6.3. Numerical results 123 The structure of this result is completely analogous to those in [110,111], and thus, the conclusions can also be drawn in a similar way. In particular, we want to highlight that (i) the factor that controls the size of the hump is proportional to aini 2−as 2and (ii) the shape of the hump is codified in the factor between brackets that only depends on β. Note that (aini 2−as 2)>0 for the considered cooling protocols (ξ1< ξ < ξ0). Thus, no anomalous Kovacs hump is expected in the nonlinear regime either. 6.3 Numerical results Here, we compare the theoretical approach above to numerical results. Specifically, we focus on the case β= 1 that gives a collision rate linear in the relative velocity and thus similar to that of hard-spheres. All simulations have been carried out with a restitution coefficient α= 0.999, which corresponds to the quasi-elastic limit in which our simplified kinetic description holds. Furthermore, we have set ω= 1 without loss of generality. 6.3.1 Validation of the first Sonine approximation First of all, we check the validity of the first Sonine approximation, as given by (6.34), to describe the time evolution of our system. In order to do so, we compare several relaxation curves between the NESS corresponding to two different noise strengths. In particular, we always depart from the stationary state corresponding to χ0= 1 and afterwards let the system evolve with χ={0.2,0.6,0.8,1}for t > 0. In Figure 6.1, we compare the Monte Carlo simulations—see once more appendix C for details—with the numerical solution of the evolution equations in the first Sonine approximation (6.34). In addition, we also have plotted the analytical solution of the linear response system, (6.40). The agreement is complete between simulation and theory, and it can be observed how the linear response result becomes more accurate as the temperature jump decreases. In order not to clutter the plot in Figure 6.1, we have not shown the results for the nonlinear first Sonine approximation (6.33). The relative error between their numerical solution and that of the standard first Sonine approximation (6.34) is at most of order 0.1%, for all the cases we have considered. Henceforth, we always use the latter, which is the usual approach in kinetic theory. 6.3.2 Kovacs Hump Since the numerical integration of the first Sonine approximation perfectly agrees with Monte Carlo simulations, we compare the analytical results for the Kovacs hump with the former. Specifically, we work in reduced variables, and therefore, we integrate numerically (6.38). Linear Response It is convenient to rewrite the expression for the Kovacs hump in an alternative form to compare our theory with numerical results. We take advantage of the simple structure of the relaxation function in the first Sonine approximation, which is the sum of two i i “thesis” — 2018/9/7 — 8:45 — page 124 — #132 i i i i i i 124 Chapter 6. Memory effects in athermal systems Figure 6.1: Direct relaxation of the granular temperature Tfor different final noise amplitudes. All curves start from the stationary state corresponding to χ0= 1. Monte Carlo simulations for a system of N= 100 sites (symbols) are compared with the numerical solution of the first Sonine approximation (6.34) (solid lines), and the analytic solution of the linear response system (6.40) (dashed lines). exponentials, to introduce the factorization [107] KT(s) = K0(sw)K1(s−sw),(6.49a) where K0(sw) = c+c− eλ+sw−eλ−sw 1−φT(sw), K1(s−sw) = eλ+(s−sw)−eλ−(s−sw).(6.49b) Firstly, this factorization property shows that the position skof the maximum relative to the waiting time sw, that is, sk−sw, is controlled by the function K1. Thus, sk−sw does not depend on the waiting time, but only on the two eigenvalues λ±. Namely, sk−sw=1 λ+−λ− ln λ− λ+≃ β=1 0.442.(6.50) Secondly, the height of the maximum Kmax does depend on the waiting time swdue to the factor K0(sw). Specifically, it can be shown that Kmax is a monotonically decreasing function of the waiting time swthat vanishes in the limit as sw→ ∞. In order to check the above results, we have fixed the initial and final drivings in the Kovacs protocol χ0and χand changed the intermediate driving value χ1. We do so to simplify the comparison, because the time scale sinvolves the steady value of the temperature; see (6.37). Note that the smaller χ1is, the shorter the waiting time becomes. Therefore, one expects to get a Kovacs hump whose maximum remains at s−sw≃0.44, i i “thesis” — 2018/9/7 — 8:45 — page 125 — #133 i i i i i i Section 6.3. Numerical results 125 Figure 6.2: Kovacs hump in linear response. Initial and final drivings, χ0= 1.05 and χ= 1 are fixed. We have considered four values for the intermediate driving, namely χ1={0.5,0.8,0.95,0.99}. The linear response result (6.49) (solid line) perfectly agrees with the numerical solution of the first Sonine approximation (6.38) (symbols). In addition, we plot the theoretical prediction for the maximum (6.50), which again agrees with the numerics (dotted line). but raises as χ1decreases. This is shown in figure 6.2, where the numerical solution of the first Sonine approximation (6.38) and the analytical result (6.49) are compared. Their agreement is almost perfect for the two lowest curves, corresponding to χ1= 0.99 and χ1= 0.95, as expected, but is still remarkably good for the two topmost ones, corresponding to the not-so-small jumps for χ1= 0.8 and χ1= 0.5. Nonlinear Regime Furthermore, we explore the Kovacs effect out of the linear regime. Figure 6.3 is similar to figure 6.2, but for larger temperature (or driving) jumps. We have also fixed the initial and final values of the driving, χ0= 10 and χ= 1. The intermediate values of the driving are the same as in the linear case except for the largest one, χ1= 0.99, which we have omitted for the sake of clarity—its hump is too small in the scale of the figure. Now, the linear response theory results just provide the qualitative behavior of the hump, correctly predicting the position of the maximum, but not its height. On the one hand, and consistently with the numerical results in an active matter model [114], the Kovacs hump out of the linear response regime is larger than the prediction of linear response theory. On the other hand, the position of the maximum remains basically unchanged, and its height still increases as χ1decreases. We also compare our analytical expansion in as 2with the numerical solutions of (6.38) for large jumps. Specifically, in order to make things as simple as possible, we choose χ1= 0. If the waiting time is long enough, the system reaches the homogeneous cooling i i “thesis” — 2018/9/7 — 8:45 — page 126 — #134 i i i i i i 126 Chapter 6. Memory effects in athermal systems Entropy 2017,19, x 15 of 19 Figure 4. (Color online) Kovacs hump out of the linear regime. The initial driving is much higher than that in Figure 3, χ0= 10, whereas the final and intermediate values of the driving are again χ= 1 and χ1={ 0.5, 0.8, 0.95 } . The linear response theoretical expression (48) (solid line) remains quite below the numerical solutions of the first Sonine approximation (42) (symbols). The theoretical expression for the maximum in linear response (54) (dotted line) still gives a good description thereof, see also Figure 5. We also compare our analytical expansion in as 2 with the numerical solutions of Equations (42) for large jumps. Specifically, in order to make things as simple as possible, we choose χ1= 0. If the waiting time is long enough, the system reaches the homogeneous cooling state and a2(sw) = aHCS 2 , which is then the initial condition for the final stage of the Kovacs protocol. Moreover, we can compute the location of the maximum of the hump from Equation (52), obtaining sk−sw=2 κ2−2−βln 2κ2 2+β≃ β=10.437. (55) This result is numerically indistinguishable from that of linear response, as given by Equation (54) , since the relative error is around 1%. In Figure 5, we put forward a comparison between the Kovacs hump obtained from the numerical solution of the first Sonine approximation equations and our theoretical expression for the nonlinear regime, Equation (52) . Fixing ξ1= 0 and ξ= 1, as ξ0 increases (as the waiting time is increased), the hump approaches Equation (52) with aini 2=aHCS 2 . Moreover, our theory perfectly reproduces all the numerical curves when we substitute the actual values of aini 2into Equation (52) . Figure 6.3: Kovacs hump in the nonlinear regime. The initial driving is much higher than that in figure 6.2, χ0= 10, whereas the final and intermediate values of the driving are the same, χ= 1 and χ1= {0.5,0.8,0.95}. The linear response theoretical expression (6.44) (solid line) seriously underestimates the actual humps, given by the numerical solutions of the first Sonine approximation (6.38) (symbols). However, the theoretical prediction for the maximum position in linear response (6.50) (dotted line) is still quite accurate, see also figure 6.4. state and a2(sw) = aHCS 2, which is then the initial condition for the final stage of the Kovacs protocol. Moreover, we can compute the location of the maximum of the hump from (6.48), obtaining sk−sw=2 κ2−2−βln 2κ2 2 + β≃ β=1 0.437.(6.51) This result is numerically indistinguishable from that of linear response, as given by (6.50), since the relative error is around 1%. In figure 6.4, we put forward a comparison between the Kovacs hump obtained from the numerical solution of the first Sonine approximation equations and our theoretical expression for the nonlinear regime, as given by (6.48). Fixing ξ1= 0 and ξ= 1, as ξ0 increases (as the waiting time is increased), the hump approaches (6.48) with aini 2=aHCS 2. Moreover, our theory perfectly reproduces all the simulations when we substitute the actual values of aini 2into (6.48). 6.3.3 Time evolution of the H-Functional The nonmonotonicity in the relaxation of the granular temperature that is brought about by the Kovacs protocol is not automatically transferred to other relevant physical mag- i i “thesis” — 2018/9/7 — 8:45 — page 127 — #135 i i i i i i Section 6.3. Numerical results 127 Figure 6.4: Kovacs hump in the nonlinear regime. We have considered the following values of the drivings: χ0={2,10,50},χ1= 0 and χ= 1. Symbols stand for the numerical solutions of the first Sonine approximation (6.38), whereas lines correspond to the theoretical expression stemming from a perturbative expansion in as 2, as given by (6.48). For the solid line, aini 2=aHCS 2, whereas we have used the value of aini 2in the numerical solution for the dashed lines. An almost perfect agreement is observed. Finally, we also mark the theoretical prediction for the maximum position in nonlinear response (6.51) (dotted line), which also shows an excellent agreement with the numerics. nitudes. Specifically, here, we deal with the H-functional H(t) = Z+∞ −∞ dvP1(v, t) log "P1(v, t) P(s) 1(v)#,(6.52) which we have analyzed in chapter 5. Therein, we have analytically proven that H(t) is a Lyapunov functional in our system for the Maxwell collision rule, which corresponds to the case β= 0. However, it is precisely for this case that the Kovacs effect vanishes, since the evolution equation for the temperature is closed. Therefore, we look into the time evolution of H(t) for the hard particle case β= 1 in this section. Note that, in addition to our rigorous proof in the model for the Maxwell case, there is strong numerical evidence of H(t) being a Lyapunov functional for granular systems [96,97]. We have computed H(t) numerically from (6.52) within the first Sonine approximation, that is, we have substituted both P1(v, t) and P(s) 1(v) by their expressions in the first Sonine approximation and calculated the integral numerically. This has been done for the Kovacs protocols considered in figures 6.2 and 6.3. The results are shown in figure 6.5 and make it clear that H(t) still monotonically decreases for the Kovacs-like protocols, in both the linear (top panel) and nonlinear (bottom panel) regimes. At the time of the maximum in the hump, s−sw≃0.44, no special signature is observed in the “entropy”, which would be given by −kBH(t). i i “thesis” — 2018/9/7 — 8:45 — page 128 — #136 i i i i i i 128 Chapter 6. Memory effects in athermal systems Figure 6.5: Time evolution of the H-functional in the Kovacs experiments. The top and bottom panels correspond, respectively, to the protocols in figures 6.2 (filled symbols) and 6.3 (open symbols), that is, to the linear and nonlinear regimes. The vertical dotted line marks the maximum of the hump in the corresponding regime. i i “thesis” — 2018/9/7 — 8:45 — page 129 — #137 i i i i i i Chapter 7 Conclusions This thesis intends to be a showcase for the relevance of simple models as useful tools for improving our understanding of complex systems. We have analyzed, mainly, two models of very different nature with the fundamental tools of the statistical mechanics. On the one hand, we have studied a model capable of predicting the deterministic unfolding pathway in modular biomolecules. On the other hand, we have thoroughly explored a model with nonconservative interactions, motivated within the field of soft matter and granular media. In this final chapter we enumerate the major conclusions of this thesis as well as some future perspectives. Regarding the first part devoted to the biophysical model: 1. We have put forward a 1d model of modular biomolecules that comprises several consecutive nonharmonic units, which is submitted to a mechanical pull that tries to control the total length of the chain. Langevin dynamics is assumed for the endpoints of the different units. 2. Dynamical equations are solved, in the limit of negligible thermal noise, by means of a perturbative approach in both the pulling velocity and the asymmetry between units. Within this quite drastic approach, the dynamics becomes deterministic. With no thermal activation being allowed, the only mechanism for unfolding is reaching the stability threshold, where the folded state ceases to exist. 3. We derive the critical values of the pulling velocity at which the unfolding pathway changes. For slow enough pulling, it is the weakest unit the first to unfold. By contrast, for fast enough pulling, the first unit to unfold is the closest to the pulled end. Intermediate situations are also possible, depending on the details of the system. Comparison between numerical solutions of particular systems and theoretical predictions shows a remarkably good agreement. 4. The applicability of the model has been carefully studied. The pulling velocity has to lie within a very specific range: low enough to allow our perturbative approach, and, at the same time, high enough to prevent thermal activation and make the 129 i i “thesis” — 2018/9/7 — 8:45 — page 136 — #144 i i i i i i 136 Bibliography [66] J. J. Brey, P. Maynar, and M. I. Garc´ıa de Soria, Phys. Rev. E 79, 051305 (2009). [67] L. Bertini, A. De Sole, D. Gabrielli, G. Jona-Lasinio, and C. Landim, Phys. Rev. Lett. 87, 040601 (2001). [68] C. Kipnis and C. Landim, Scaling Limits of Interacting Particle Systems, (Springer- Verlag, Berlin, 1999). [69] C. Kipnis, C. Marchioro, and E. Presutti, J. Stat. Phys. 27, 65 (1982). [70] P. I. Hurtado and P. L. Garrido, Phys. Rev. Lett. 102, 250601 (2009). [71] P. I. Hurtado and P. L. Garrido, Phys. Rev. E 81, 041102 (2010). [72] P. I. Hurtado and P. L. Garrido, J. Stat. Mech. Theory Exp. 2009, P02032 (2009). [73] P. I. Hurtado and P. L. Krapivsky, Phys. Rev. E 85, 060103 (2012). [74] Y. Srebro and D. Levine, Phys. Rev. Lett. 93, 240601 (2004). [75] A. Prados, A. Lasanta, and P. I. Hurtado, Phys. Rev. E 86, 031134 (2012). [76] A. Prados, A. Lasanta, and P. I. Hurtado, Phys. Rev. Lett. 107, 140601 (2011). [77] P. I. Hurtado, A. Lasanta, and A. Prados, Phys. Rev. E 88, 022110 (2013). [78] A. Lasanta, P. I. Hurtado, and A. Prados, Eur. Phys. J. E 39 (2016). [79] A. Lasanta, A. Manacorda, A. Prados, and A. Puglisi, New J. Phys. 17, 083039 (2015). [80] V. Garz´o and A. Santos, Kinetic Theory of Gases in Shear Flows: Nonlinear Transport, (Springer, Netherlands, 2003). [81] P. K. Haff, J. Fluid Mech. 134, 401430 (1983). [82] M. Ernst, Phys. Rep. 78, 1 (1981). [83] J. J. Brey, M. J. Ruiz-Montero, and D. Cubero, Phys. Rev. E 54, 3664 (1996). [84] J. J. Brey, A. Prados, M. I. Garc´ıa de Soria, and P. Maynar, J. Phys. A 40, 14331 (2007). [85] S. McNamara, Phys. Fluid. Fluid Dynam. 5, 3056 (1993). [86] V. Garz´o, Phys. Rev. E 73, 021304 (2006). [87] A. M. Lyapunov, Int. J. Control 55, 531 (1992). [88] L. Boltzmann, Lectures on Gas Theory, (Dover Publications, New York, 1995). [89] J. L. Lebowitz, Physica A 194, 1 (1993). [90] J. L. Lebowitz, Phys. Today 46, 32 (1993). i i “thesis” — 2018/9/7 — 8:45 — page 137 — #145 i i i i i i Bibliography 137 [91] J. L. Lebowitz, Physica A 263, 516 (1999). [92] I. Prigogine, Physica A 263, 528 (1999). [93] D. Ruelle, Physica A 263, 540 (1999). [94] S. Chapman, T. Cowling, D. Burnett, and C. Cercignani, The Mathematical Theory of Non-uniform Gases: An Account of the Kinetic Theory of Viscosity, Thermal Conduction and Diffusion in Gases, (Cambridge University Press, Cambridge, New York, 1990). [95] N. G. Van Kampen, Stochastic Processes in Physics and Chemistry, (North-Holland, Amsterdam, 1992). [96] U. M. B. Marconi, A. Puglisi, and A. Vulpiani, J. Stat. Mech. Theory Exp. 2013, P08003 (2013). [97] M. I. Garc´ıa de Soria, P. Maynar, S. Mischler, C. Mouhot, T. Rey, and E. Trizac, J. Stat. Mech. Theory Exp. 2015, P11009 (2015). [98] C. Villani, J. Stat. Phys. 124, 781 (2006). [99] I. Bena, F. Coppex, M. Droz, P. Visco, E. Trizac, and F. van Wijland, Physica A 370, 179 (2006). [100] G. A. Bird, Phys. Fluids 6, 1518 (1963). [101] G. Bird, The DSMC Method, (CreateSpace Independent Publishing Platform, USA, 2013). [102] A. J. Kovacs, J. J. Aklonis, J. M. Hutchinson, and A. R. Ramos, J. Polym. Sci., Polym. Phys. Ed. 17, 1097 (1979). [103] E. M. Bertin, J. P. Bouchaud, J. M. Drouffe, and C. Godrche, J. Phys. A 36, 10701 (2003). [104] A. Buhot, J. Phys. A 36, 12367 (2003). [105] S. Mossa and F. Sciortino, Phys. Rev. Lett. 92, 045504 (2004). [106] G. Aquino, L. Leuzzi, and T. M. Nieuwenhuizen, Phys. Rev. B 73, 094205 (2006). [107] A. Prados and J. J. Brey, J. Stat. Mech. Theory Exp. 2010, P02009 (2010). [108] G. Diezemann and A. Heuer, Phys. Rev. E 83, 031505 (2011). [109] M. Ruiz-Garc´ıa and A. Prados, Phys. Rev. E 89, 012140 (2014). [110] A. Prados and E. Trizac, Phys. Rev. Lett. 112, 198001 (2014). [111] E. Trizac and A. Prados, Phys. Rev. E 90, 012204 (2014). [112] J. J. Brey, M. I. Garc´ıa de Soria, P. Maynar, and V. Buz´on, Phys. Rev. E 90, 032207 (2014). i i “thesis” — 2018/9/7 — 8:45 — page 138 — #146 i i i i i i 138 Bibliography [113] Y. Lahini, O. Gottesman, A. Amir, and S. M. Rubinstein, Phys. Rev. Lett. 118, 085501 (2017). [114] R. K¨ursten, V. Sushkov, and T. Ihle, Phys. Rev. Lett. 119, 188001 (2017). [115] G. Hummer and A. Szabo, Biophys. J. 85, 515 (2003). [116] G. Arad-Haase, S. G. Chuartzman, S. Dagan, R. Nevo, M. Kouza, B. K. Mai, H. T. Nguyen, M. S. Li, and Z. Reich, Biophys. J. 99, 238 (2010). [117] C. M. Bender and S. A. Orszag, Advanced mathematical methods for scientists and engineers I: Asymptotic methods and perturbation theory, (Springer, New York, 1999). [118] F. Rico, L. Gonzalez, I. Casuso, M. Puig-Vidal, and S. Scheuring, Science 342, 741 (2013). [119] O. K. Dudko, G. Hummer, and A. Szabo, Phys. Rev. Lett. 96, 108101 (2006). [120] C. A. Plata, F. Cecconi, M. Chinappi, and A. Prados, J. Stat. Mech. Theory Exp. 2015, P08003 (2015). [121] R. Berkovich, S. Garcia-Manyes, M. Urbakh, J. Klafter, and J. M. Fernandez, Biophys. J. 98, 2692 (2010). [122] R. Berkovich, R. I. Hermans, I. Popa, G. Stirnemann, S. Garcia-Manyes, B. J. Berne, and J. M. Fernandez, Proc. Natl. Acad. Sci. U.S.A. 109, 14416 (2012). [123] R. Berkovich, S. Garcia-Manyes, M. Urbakh, J. Klafter, and J. M. Fernandez, Biophys. J. 98, 2692 (2010). [124] Q. Li, Z. N. Scholl, and P. E. Marszalek, Angew. Chem. Int. Ed. 53, 13429 (2014). [125] I. Schwaiger, A. Kardinal, M. Schleicher, A. A. Noegel, and M. Rief, Nat. Struct. Mol. Biol. 11, 81 (2004). [126] H. Grubm¨uller, B. Heymann, and P. Tavan, Science 271, 997 (1996). [127] H. Lu, B. Isralewitz, A. Krammer, V. Vogel, and K. Schulten, Biophys. J. 75, 662 (1998). [128] J. N. Israelachvili, Intermolecular and Surface Forces (Third Edition), (Academic Press, San Diego, 2011). [129] D. Zuckerman, Statistical Physics of Biomolecules: An Introduction, (CRC Press, Boca Raton, Florida, 2010). [130] D. Frenkel and B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, (Academic Press, San Diego, 2001). [131] D. A. Parry, R. B. Fraser, and J. M. Squire, J. Struct. Biol. 163, 258 (2008). i i “thesis” — 2018/9/7 — 8:45 — page 139 — #147 i i i i i i Bibliography 139 [132] O. J. Rackham, M. Madera, C. T. Armstrong, T. L. Vincent, D. N. Woolfson, and J. Gough, J. Mol. Biol. 403, 480 (2010). [133] L. Truebestein and T. A. Leonard, BioEssays 38, 903 (2016). [134] M. Aittaleb, R. Rashid, Q. Chen, J. R. Palmer, C. J. Daniels, and H. Li, Nat. Struct. Biol. 10, 256 (2003). [135] W. Humphrey, A. Dalke, and K. Schulten, J. Mol. Graph. 14, 33 (1996). [136] J. C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chipot, R. D. Skeel, L. Kal, and K. Schulten, J. Comput. Chem. 26, 1781 (2005). [137] R. B. Best, G. Hummer, and W. A. Eaton, Proc. Natl. Acad. Sci. U.S.A. 110, 17874 (2013). [138] C. A. Plata, Z. N. Scholl, P. E. Marszalek, and A. Prados, J. Chem. Theory Comput. 14, 2910 (2018). [139] H. Dietz, F. Berkemeier, M. Bertz, and M. Rief, Proc. Natl. Acad. Sci. U.S.A. 103, 12724 (2006). [140] Y. Gao, G. Sirinakis, and Y. Zhang, J. Am. Chem. Soc. 133, 12749 (2011). [141] A. Baldassarri, U. M. B. Marconi, and A. Puglisi, EPL 58, 14 (2002). [142] G. Grinstein, D.-H. Lee, and S. Sachdev, Phys. Rev. Lett. 64, 1927 (1990). [143] P. L. Garrido, J. L. Lebowitz, C. Maes, and H. Spohn, Phys. Rev. A 42, 1954 (1990). [144] H. Spohn, Rev. Mod. Phys. 52, 569 (1980). [145] A. Manacorda, C. A. Plata, A. Lasanta, A. Puglisi, and A. Prados, J. Stat. Phys. 164, 810 (2016). [146] A. Manacorda, Lattice Models for Fluctuating Hydrodynamics in Granular and Active Matter, Springer Theses, (Springer, 2018). [147] E. Ben-Naim and P. L. Krapivsky, in T. P¨oschel and N. Brilliantov (editors), Granular Gas Dynamics, volume 624 of Lecture Notes in Physics, 65–94, (Springer, Berlin, 2003). [148] M. H. Ernst, E. Trizac, and A. Barrat, EPL 76, 56 (2006). [149] A. Baldassarri, A. Puglisi, and A. Prados, Phys. Rev. E 97, 062905 (2018). [150] C. A. Plata, A. Manacorda, A. Lasanta, A. Puglisi, and A. Prados, J. Stat. Mech. Theory Exp. 2016, 093203 (2016). [151] T. P. C. van Noije, M. H. Ernst, R. Brito, and J. A. G. Orza, Phys. Rev. Lett. 79, 411 (1997). i i “thesis” — 2018/9/7 — 8:45 — page 140 — #148 i i i i i i 140 Bibliography [152] M. I. Garc´ıa de Soria, P. Maynar, G. Schehr, A. Barrat, and E. Trizac, Phys. Rev. E77, 051127 (2008). [153] A. W. Lees and S. F. Edwards, J. Phys. C 5, 1921 (1972). [154] A. Santos and V. Garz´o, J. Stat. Mech. Theory Exp. 2007, P08021 (2007). [155] A. Santos, V. Garz´o, and J. W. Dufty, Phys. Rev. E 69, 061303 (2004). [156] P. Resibois and M. de Leener, Classical Kinetic Theory of Fluids, (Wiley, New York, 1977). [157] M. I. Garc´ıa de Soria, P. Maynar, and E. Trizac, Phys. Rev. E 85, 051301 (2012). [158] M. I. Garc´ıa de Soria, P. Maynar, and E. Trizac, Phys. Rev. E 87, 022201 (2013). [159] J. J. Brey, A. Dom´ınguez, M. I. Garc´ıa de Soria, and P. Maynar, Phys. Rev. Lett. 96, 158002 (2006). [160] V. V. Prasad, S. Sabhapandit, A. Dhar, and O. Narayan, Phys. Rev. E 95, 022115 (2017). [161] J. J. Brey and A. Prados, Phys. Rev. E 47, 1541 (1993). [162] J. J. Brey and A. Prados, Phys. Rev. B 49, 984 (1994). [163] J. Brey, A. Prados, and M. Ruiz-Montero, J. Non-Cryst. Solids 172-174, 371 (1994). [164] M. O. Vlad and J. Ross, J. Phys. Chem. B 101, 8756 (1997). [165] M. O. Vlad, F. Moran, and J. Ross, J. Phys. Chem. B 102, 4598 (1998). [166] A. Prados, J. Brey, and B. S´anchez-Rey, Physica A 284, 277 (2000). [167] B. Earnshaw and J. Keener, SIAM J. Appl. Dyn. Syst. 9, 220 (2010). [168] J. M. Montanero and A. Santos, Granul. Matter 2, 53 (2000). [169] M. I. Garc´ıa de Soria, P. Maynar, and E. Trizac, Mol. Phys. 107, 383 (2009). [170] I. Pagonabarraga, E. Trizac, T. P. C. van Noije, and M. H. Ernst, Phys. Rev. E 65, 011303 (2001). [171] P. Maynar, M. d. Garc´ıa de Soria, and E. Trizac, Eur. Phys. J. Spec. Top. 179, 123 (2009). [172] V. V. Prasad, S. Sabhapandit, and A. Dhar, EPL 104, 54003 (2013). [173] P. L. Chebyshev, Acta Math. 14, 305 (1890). [174] C. V. L. Charlier, Ark. Math. Astr. och Phys. 2, 1 (1905-06). [175] F. Y. Edgeworth, Math. Proc. Cambridge Philos. Soc. 20, 36 (1905). [176] D. L. Wallace, Ann. Math. Stat. 29, 635 (1958). i i “thesis” — 2018/9/7 — 8:45 — page 141 — #149 i i i i i i Bibliography 141 [177] M. Abramowitz and I. A. Stegun (editors), Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, (Dover, Washington DC, USA, 1972). [178] H. Cram´er, in Proceedings of the Sixth Scandinavian Congress of Mathematicians, Copenhagen, 399–425 (1925). [179] G. Szeg¨o, Orthogonal Polynomials, (American Mathematical Society, New York, 1939). [180] J. Brey and A. Prados, Physica A 197, 569 (1993). [181] I. M¨uller and P. Villaggio, Arch. Ration. Mech. Anal. 65, 25 (1977). [182] G. Puglisi and L. Truskinovsky, Cont. Mech. Thermodyn. 14, 437 (2002). [183] W. Dreyer, J. Jamnik, C. Guhlke, R. Huth, J. Mokon, and M. Gaberek, Nat. Mater. 9, 448 (2010). [184] W. Dreyer, C. Guhlke, and M. Herrmann, Cont. Mech. Thermodyn. 23, 211 (2011). [185] I. Benichou, E. Faran, D. Shilo, and S. Givli, Appl. Phys. Lett. 102, 011912 (2013). [186] L. Bertini, A. De Sole, D. Gabrielli, G. Jona-Lasinio, and C. Landim, Rev. Mod. Phys. 87, 593 (2015). [187] A. Bortz, M. Kalos, and J. Lebowitz, J. Comput. Phys. 17, 10 (1975). [188] A. Prados, J. J. Brey, and B. S´anchez-Rey, J. Stat. Phys. 89, 709 (1997). i i “thesis” — 2018/9/7 — 8:45 — page 142 — #150 i i i i i i 142 Bibliography i i “thesis” — 2018/9/7 — 8:45 — page 143 — #151 i i i i i i Appendix A Stability threshold To first order in ξ, the extension xi,b such that a00 i(xi,b) = 0 verifies a000(`b)(xi,b −`b) + ξδf0 i(`b)=0,(A.1) that is, xi,b =`b−ξδf0 i(`b) a000(`b).(A.2) The corresponding force at the stability threshold is obtained from (2.11). To the lowest order in the deviations, Fi,b ≡a0 i(xi,b)∼a0 i(`b) = Fb+ξ δfi(`b),(A.3) because the next term, a000(`b)(xi,b −`b)2/2, is of the order of ξ2. Therefore, the i-th module reaches its limit of stability at the time for which xi=x(0) i+ξδxi=xi,b, that is, when the length per unit `has the value `iverifying `i+ξδf(`i)−δfi(`i) a00(`i)=`b−ξδf0 i(`b) a000(`b),(A.4) or, equivalently, `i−`b=ξδfi(`i)−δf(`i) a00(`i)−ξδf0 i(`b) a000(`b).(A.5) We know that `→`bwhen ξ→0. But a00(`b) = 0 and thus we cannot substitute `bon the rhs of (A.5). On the other hand, this means that the dominant balance for ξ→0 involves the lhs and the first term on the rhs of (A.5). Therefore, making use of a00(`)∼a000(`b)(`−`b), we get (`i−`b)2∼ξδfi(`b)−δf(`b) a000(`b).(A.6) Since a000(`b)<0 (see figure 2.2), this means that only the units with δfi(`b) smaller than the average (that is, weaker than average) reach the limit of stability in the limit as 143 i i “thesis” — 2018/9/7 — 8:45 — page 144 — #152 i i i i i i 144 Appendix A. Stability threshold vp→0. In fact, it is the weakest unit, that is, the unit with smallest δfi(`b), that unfolds first. It is interesting to note that, in order to obtain (A.6), we have completely neglected the last term on the rhs of (A.5). Since, in turn, this term stems from the last term on the rhs of (A.2), to the lowest order we are solving the equation xi=`b. In other words, to the lowest order the stability threshold can be considered to be given by the non-disordered, zero-asymmetry case, free energy a(x). For the sake of concreteness and simplicity, we have stuck to the asymmetry contribution δxiin this appendix, but the same condition xi=`bwould still be valid, had we taken into account the kinetic contribution ∆xi derived in section 2.2.2. The reason is that there is also a factor a00(`) in the denominator of ∆xi, see (2.24a), and thus both the terms coming from δxiand ∆xiare dominant against the last term on the rhs of (A.2). i i “thesis” — 2018/9/7 — 8:45 — page 145 — #153 i i i i i i Appendix B Fluctuating expression for the dissipation Let us consider the dissipation dl,p at site land at time p. Its main part is obtained by averaging (4.6) with respect to the fast variables yp, i.e. ¯ dl,p =α2−1 4L∆2 l,p + ∆2 l−1,p<0.(B.1) This is the expression that we have to analyze in the fluctuating hydrodynamic description, since the amplitude of the dissipation noise scales as L−3[145,146]. If we consider the average of the dissipation field, it is readily obtained that dav l,p = (α2−1)Tl,p/L, which gives (4.13) in the continuum limit by using d(x, t) = L3dl,p and the definition of ν. Therefore, it is consistent to write at the fluctuating level that ¯ dl,p =α2−1 Lθl,p,(B.2) by defining the fluctuating temperature as θl,p =1 4∆2 l,p + ∆2 l−1,p=v2 l−1,p + 2v2 l,p +v2 l+1,p 4−vl,p vl+1,p +vl−1,p 2.(B.3) The first term on the rhs, (v2 l−1,p +2v2 lp +v2 l+1,p)/4, reduces to el,p plus terms of the order of L−2, which are neglected. Our main goal is to obtain a correct expression for vl,pvl±1,p at the fluctuating level. In general, we have for the average correlations hvl,pvl0,pi=El,pδll0+Cl,l0−l;p(1 −δll0),(B.4) with the definition Cl,l0−l;p=hvl,pvl0,pifor l06=l. (B.5) 145