E Bond correlations and internal distances F Partial structure factor fits

Chapter 1 Introduction

This thesis is about the temperature dependence of dynamics in polymer melts, and how this mobility may change when polymer melts are confined to thin films. Let us start by discussing what polymer melts are about.

Polymer melts are systems that contain many macromolecules in form of a liquid [START_REF] Doi | The theory of polymer dynamics[END_REF]2]. In the simplest case, macromolecules are linear chains comprised of 𝑁 monomeric repeat units (number 𝑁 is usually called a degree of polymerization or, simply, the chain length). The chains can be up to 𝑁 ≲ 10 5 long making them quite big objects with the gyration radius 𝑅 g ≲ 1 𝜇m, even though the size of a monomer is only ∼ 5 Å. Such a large span of length scales is pronounced in the structural and dynamical properties of polymer melts. Being liquids, they show an amorphous dense packing of the monomers on a local scale and low compressibility on a macroscopic scale [START_REF] Hansen | Theory of simple liquids[END_REF]. On the intermediate scales (up to ∼ 𝑅 g ) the structure of the melt is determined by the intra-and intermolecular contributions. The former manifests in the fact that long polymers are self-similar fractal objects [2,[START_REF] De Gennes | Scaling concepts in polymer physics[END_REF] that have the ability to adopt different conformations. Such a chain structure allows the other chains to penetrate inside the volume it occupies leading to on average √ 𝑁 other chains inside the pervaded volume by one chain. These neighboring chains screen the intrachain excluded volume interactions which in a dilute solution would lead to swelling of the polymer [START_REF] Doi | The theory of polymer dynamics[END_REF]2,[START_REF] De Gennes | Scaling concepts in polymer physics[END_REF]. Thus, to a good approximation the chains in a melt have a random-coil-like conformation with the size scaling as 𝑅 g ∼ √ 𝑁 with the chain length. The dynamics of the melt also varies with the length scale: the (longest) relaxation time of a polymer 𝜏 𝑁 is much larger than the monomer relaxation time 𝜏 1 . For moderate 𝑁 it is commonly described by the Rouse model [START_REF] Doi | The theory of polymer dynamics[END_REF]2,[START_REF] De Gennes | Scaling concepts in polymer physics[END_REF][START_REF] Rouse | A theory of the linear viscoelastic properties of dilute solutions of coiling polymers[END_REF], that predicts 𝜏 𝑁 = 𝜏 1 𝑁 2 .

CHAPTER 1. INTRODUCTION

Due to the slow relaxations, polymer melts are more viscous than low-molecular-weight liquids. Moreover, since the chains cannot cross each other, a network of topological constraints appears and persists during significantly long times when the chain length exceeds a certain value 𝑁 e . These entanglements slow down the dynamics of the melt even more. The so-called "reptation" model describes this effect and predicts 𝜏 𝑁 = 𝜏 1 𝑁 3 /𝑁 e [START_REF] Doi | The theory of polymer dynamics[END_REF]2,[START_REF] De Gennes | Scaling concepts in polymer physics[END_REF]. In addition to slowing down the relaxation in the time window 𝜏 1 < 𝑡 < 𝜏 𝑁 , the increase of the chain length leads to a pronounced viscoelastic behavior of the melt even at high temperatures.

Upon cooling to low temperatures, polymer liquids undergo a transition to a solid state, which is either semicrystalline or glassy (amorphous) depending on the structure of the chain. Only polymers with side groups having regular tacticity [2] (e.g. isotactic/syndiotactic polystyrene) or without side groups (ex. polyethylene) can form crystalline lamellae. Still, the full crystallization is almost never achieved and the lamellar sheets exist along with amorphous regions [START_REF] Strobl | The physics of polymers: concepts for understanding their structures and behavior[END_REF][START_REF] Muthukumar | Nucleation in polymer crystallization[END_REF]. For atactic (homo)polymers (e.g. atactic polystyrene) or random copolymers (𝑐𝑖𝑠-𝑡𝑟𝑎𝑛𝑠 polybutadiene, etc.) the crystalline phase does not exist. Upon cooling these polymer melts are in general good glass formers [START_REF] Mckenna | Comprehensive polymer science[END_REF][START_REF] Donth | The glass transition[END_REF].

Similarly to nonpolymeric glass-forming liquids, with decreasing the temperature towards a characteristic value 𝑇 g , all structural relaxation times drastically increase over about 14 orders of magnitude (from ∼ 1 ps at high 𝑇 to ∼ 100 s at 𝑇 g ) in a super-Arrhenius fashion. Below 𝑇 g the relaxation is so slow, that the system is indistinguishable from an amorphous solid. On the contrary, static two-point correlation functions, such as the static structure factor 𝑆(𝑞), change mildly and below 𝑇 g the system has the structure of a "frozen" liquid, that has stopped to flow [START_REF]Polymer glasses[END_REF]. In addition to slowing down, another characteristic feature of the glassy dynamics is a two-step relaxation of the time correlation functions 𝐶(𝑡) = ⟨𝐴 * (𝑡)𝐴(0)⟩ 1 of the structure-sensitive variables 𝐴, with the first step attributed to the weekly temperature dependent microscopic relaxation and the second step corresponding to structural relaxation. The latter one is stretched over large time windows and has a strong dependence on temperature. Typical experimentally observed quantities that show this behavior are the coherent/incoherent density fluctuations measured in dynamic neutron scattering experiments [START_REF] Arbe | Insight into the structure and dynamics of polymers by neutron scattering combined with atomistic molecular dynamics simulations[END_REF]. These effects are also observed as broadening of the dominant loss peak in the frequency domain when measuring the relaxation of the 1 𝐴 * denotes the complex conjugation and the angular brackets ⟨⟩ denote the thermal average.

local electric dipoles in dielectric spectroscopy measurements [START_REF] Lunkenheimer | Dielectric spectroscopy of glass-forming materials: -relaxation and excess wing[END_REF].

Even though the behavior of glass forming polymer melts outlined above has been observed over many decades already [START_REF]Polymer glasses[END_REF], the description of the glassy dynamics at the present time still has rather an empirical character and the mechanisms of the glass transition phenomenon are not yet fully understood. An overview of theoretical perspectives on the glass formation can be found in Refs. [START_REF] Berthier | Theoretical perspective on the glass transition and amorphous materials[END_REF][START_REF] Langer | Theories of glass formation and the glass transition[END_REF]. The situation is even more complicated for the polymer glasses in confinement to the nanoscale dimension, that have extensive industrial application (coatings, semiconductors, polymer nanocomposites) [START_REF] Napolitano | Glass transition of polymers in bulk, confined geometries, and near interfaces[END_REF]. With the current trend to the miniaturization, it is important to understand the effects of confinement on mechanical stability. In particular, many studies on ∼ 10-100 nm thick polymer films reported significant shifts in the average 𝑇 g with respect to the bulk value (see [START_REF]Polymer glasses[END_REF] and references therein and Fig. 1.1). These shifts were found to depend on the geometry of the film (capped/supported/freestanding), the chemical structure of the polymer and even the experimental technique employed in the investigation. Interpretation of a huge amount of data seems to reach a consensus in the fact, that these shifts are the average effect of local 𝑇 g gradients created by the confinement. In particular, the presence of a free surface for the supported and freestanding films is commonly assumed to speed up the dynamics of the close-by polymer layer leading to a decrease of the local 𝑇 g . The presence of a substrate for the capped/supported geometries is considered to have either a slowing down (increase of 𝑇 g ) or a speeding up (decrease of 𝑇 g ) effect depending on the nature of its interaction with the polymer and the roughness of the substrate. Yet, there is still an ongoing debate whether the results of the non-equilibrium and equilibrium experimental techniques can in principle be compared, since they might probe different properties [START_REF] Kremer | Glassy dynamics and glass transition in nanometric layers and films: a silver lining on the horizon[END_REF].

While it is relatively straightforward to measure the film averaged properties using a wide amount of experimental techniques (ellipsometry, calorimetry, dielectric spectroscopy etc.), it is extremely challenging to access experimentally the dynamical gradients created by the surface effects, especially on a scale of the monomer size relevant to the glass transition process. Insights on the magnitude of the surface effects and how deep into the film the dynamics is perturbed (along with the general mechanisms of the structural relaxation) can be obtained from the molecular dynamics computer simulations [START_REF] Allen | Computer simulation of liquids[END_REF].

Moreover, they allow to measure many different observables on the same sample providing Figure 1.1: Number of publications vs. shift in 𝑇 g reported therein for the thinnest studied layer of polystyrene (the figure is taken from Ref. [START_REF] Kremer | Glassy dynamics and glass transition in nanometric layers and films: a silver lining on the horizon[END_REF]). The collection is restricted to molecular weights in the range of 102-104 kg/mol and thicknesses between 5 and 50 nm; the employed experimental method and type of support are indicated. In total 92 publications are included (see [START_REF] Kremer | Glassy dynamics and glass transition in nanometric layers and films: a silver lining on the horizon[END_REF] and references therein). The number of publications on each type of support is indicated in brackets after the type of support; 5 papers employing more than one type of support are recorded in each corresponding category. The increment in the y-axis representing one publication is always one which means that for a publication covering several methods the representative area is split into as many parts.

a link between different experimental techniques and allowing to interpret the vast variety of results.

Concerning the modelling of polymer systems, multiple approaches exist that differ in the level of chemical details. The most generic approach grasps only the basic features of polymers: chain connectivity, excluded volume interactions and, possibly, some chain stiffness [START_REF] Baschnagel | Computer simulations of supercooled polymer melts in the bulk and in confined geometry[END_REF][START_REF] Hsu | A coarse-grained polymer model for studying the glass transition[END_REF]. The polymers are modeled as chains of beads (connected by springs), each of them representing a part of the real chain, coarse-grained to the scale of a Kuhn segment [2]. The size of the segment is the only intrinsic length scale in the system. Being thus independent of chemical details, these generic bead-spring models are used to study universal properties of polymers [START_REF] Kremer | Dynamics of entangled linear polymer melts: a molecular-dynamics simulation[END_REF][START_REF] Svaneborg | Characteristic time and length scales in melts of kremer-grest bead-spring polymers with wormlike bending stiffness[END_REF] including the mechanisms of glass formation [START_REF] Bennemann | Molecular-dynamics simulations of the thermal glass transition in polymer melts: 𝛼-relaxation behavior[END_REF][START_REF] Bennemann | Molecular-dynamics simulation of a glassy polymer melt: rouse model and cage effect[END_REF][START_REF] Peter | Thickness-dependent reduction of the glass-transition temperature in thin polymer films with a free surface[END_REF]. The more thorough quantitative comparison to real polymers requires a more detailed approach in the simulations using atomistic models [START_REF] Paul | Structure and dynamics of amorphous polymers: computer simulations compared to experiment and theory[END_REF][START_REF] Smith | United atom force field for molecular dynamics simulations of 1,4polybutadiene based on quantum chemistry calculations on model molecules[END_REF] that treat every atom present in the polymer (all-atom models) or combine few of them (typically carbon-hydrogen groups) together into united atoms (united atom models). In this case the details of chemical interactions matter, and in addition to the generic features mentioned above, the models include intramolecular potentials for the torsional angles and even sometimes the long range electrostatic interactions.

At the present moment most of the state of the art realistic simulations are able to follow the dynamics only up to few microseconds of real time. The bead-spring models win additional few orders of magnitude due to their simplicity. Although this time scale is much smaller than the relaxation time near 𝑇 g (which is about 100 s), simulations are appreciated [START_REF] Ediger | Dynamics near free surfaces and the glass transition in thin polymer films: a view to the future[END_REF] because they can provide a tomographic picture of the glass former and allow direct access to its properties on the relevant length scale. This ability is heavily used in the current research on the glass transition in general [START_REF] Berthier | Theoretical perspective on the glass transition and amorphous materials[END_REF] and has also been exploited

for systems in confinement [START_REF]Polymer glasses[END_REF][START_REF] Baschnagel | Computer simulations of supercooled polymer melts in the bulk and in confined geometry[END_REF][START_REF] Peter | Thickness-dependent reduction of the glass-transition temperature in thin polymer films with a free surface[END_REF][START_REF] Hanakata | Interfacial mobility scale determines the scale of collective motion and relaxation rate in polymer films[END_REF][START_REF] Yelash | Slow process in confined polymer melts: layer exchange dynamics at a polymer solid interface[END_REF][START_REF] Solar | The dielectric axation in polymer films: a comparison between experiments and atomistic simulations[END_REF].

In general the simulation results corroborate the experimental findings [START_REF] Peter | Thickness-dependent reduction of the glass-transition temperature in thin polymer films with a free surface[END_REF][START_REF] Solar | The dielectric axation in polymer films: a comparison between experiments and atomistic simulations[END_REF][START_REF] Phan | Dynamic gradients, mobile layers, 𝑇 g shifts, role of vitrification criterion, and inhomogeneous decoupling in free-standing polymer films[END_REF].

Smooth, repulsive or weakly attractive substrates or free surfaces enhance the relaxation relative to the bulk. Whereas strongly attractive substrates cause the retarding effect.

These interface-induced perturbations smoothly transition from the boundaries to the interior of the confined liquid. The range of these gradients grows on cooling so that the perturbations can propagate across the entire film for sufficiently low 𝑇 and strong confinement. Contributions from all layers in these gradients determine the average behavior of the films and so the shift of 𝑇 g . Qualitatively such a behavior is observed in the simulations using both realistic and coarse-grained models. However, the quantitative description obtained from the realistic models has been more and more often reported to disagree with the results extracted from generic simulations. The reasons of the disagreements have been attributed to the importance of the intramolecular forces, that are present in the realistic models in addition to the bead-spring-like interactions [START_REF] Colmenero | Are polymers standard glass-forming systems? The role of intramolecular barriers on the glass-transition phenomena of glass-forming polymers[END_REF]. In particular, extensive molecular dynamics simulations have recently been carried out on glassy polymer films using an atomistic model for polybutadiene (PBD) embedded between two graphite walls [START_REF] Solar | The dielectric axation in polymer films: a comparison between experiments and atomistic simulations[END_REF][START_REF] Solar | Polymer dynamics in a polymer-solid interphase: molecular dynamics simulations of 1,4-polybutadiene at a graphite surface[END_REF][START_REF] Solar | Dielectric -relaxation of 1,4-polybutadiene confined between graphite walls[END_REF][START_REF] Solar | Relaxation processes and glass transition of confined polymer melts: a molecular dynamics simulation of 1,4-polybutadiene between graphite walls[END_REF]. By calculating the dielectric response the simulations find, in qualitative agreement with experiments of the group of F. Kremer [START_REF] Mapesa | Segmental and chain dynamics in nanometric layers of poly(cis-1, 4-isoprene) as studied by broadband dielectric spectroscopy and temperature-modulated calorimetry[END_REF][START_REF] Tress | Glassy dynamics and glass transition in nanometric thin layers of polystyrene[END_REF], that confinement has no effect on the film dynamics, except in a layer of about 1 nm thickness near the substrate. It is argued that this result is divergent relative to the bead-spring CHAPTER 1. INTRODUCTION models and is caused by the dominance of intramolecular barriers. Since such barriers are present in all real polymers, references [START_REF] Solar | The dielectric axation in polymer films: a comparison between experiments and atomistic simulations[END_REF][START_REF] Colmenero | Are polymers standard glass-forming systems? The role of intramolecular barriers on the glass-transition phenomena of glass-forming polymers[END_REF][START_REF] Solar | Dielectric -relaxation of 1,4-polybutadiene confined between graphite walls[END_REF][START_REF] Solar | Relaxation processes and glass transition of confined polymer melts: a molecular dynamics simulation of 1,4-polybutadiene between graphite walls[END_REF] represent an important warning sign, challenging simulation efforts based on the generic models and the unifying physical picture emerging from these simulations for the glassy dynamics of polymers.

Clearly, this warning sign calls for a systematic study of how the intramolecular barriers and substrate interactions influence the glassy dynamics in atomistic simulations.

Our work is focused on one particular type of the intramolecular interactions, the torsion (dihedral) barriers, that control the conformations of a polymer chain.

We start from simulations using the chemically realistic united atom model of PBD [START_REF] Smith | United atom force field for molecular dynamics simulations of 1,4polybutadiene based on quantum chemistry calculations on model molecules[END_REF][START_REF] Smith | Molecular dynamics of a 1,4-polybutadiene melt. comparison of experiment and simulation[END_REF] (described in Chapter 2) and make a crossover to the model where all the dihedral barriers were switched off, complementing previous studies of how they influence the structure and the glassy dynamics of the PBD melt [START_REF] Krushev | The role of internal rotational barriers in polymer melt chain dynamics[END_REF][START_REF] Krushev | Intramolecular caging in polybutadiene due to rotational barriers[END_REF][START_REF] Bedrov | Molecular dynamics simulation study of the 𝛼-and 𝛽-relaxation processes in a realistic model polymer[END_REF] (Chapter 3). Then we employ these models in the simulations of supported PBD films (Chapter 4). Comparing between them will reveal how the interplay of intermolecular packing, intrachain barriers and surface effects influences the glass-forming properties of the PBD film. Chapter 5 provides a synthesis of the key results and conclusions of this work. The French summary of the thesis follows this conclusion chapter. Finally, a few appendices provide further technical details.

Chapter 2

Simulation model and computational details

Introduction

The computational technique used in our work to simulate polybutadiene (PBD) is classical molecular dynamics (MD) simulations [START_REF] Allen | Computer simulation of liquids[END_REF]. In its essence, the algorithm integrates Newton's equations of motion for all particles in the studied system: the positions and velocities of the particles at the next timestep are calculated from the positions and velocities at the precedent timestep using the interaction potentials defined by the user.

This set of interaction potentials along with the values of their parameters is commonly referred to as the force field. The velocities and/or simulation domain dimensions can be coupled to additional dynamic variables in order to model the effects of a thermostat and/or a barostat [START_REF] Lammps | Large-scale Atomic/Molecular Massively Parallel Simulator[END_REF]. This determines the thermodynamic ensemble in which the coordinates and velocities are generated. The simulations last until a sufficient number of timesteps has been reached. During a simulation the coordinates and velocities of the particles are saved at selected timestep intervals into the so-called simulation trajectory, which is used during the post-analysis to compute the quantities of interest.

Our simulations were conducted using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [START_REF] Lammps | Large-scale Atomic/Molecular Massively Parallel Simulator[END_REF][START_REF] Plimpton | Fast parallel algorithms for short-range molecular dynamics[END_REF], which is a high performance computer code that handles parallel computing, provides extensive force field customization and a variety of tools for efficient computation of physical quantities during ongoing simulations. Additional benefits of using LAMMPS were the broad expertise and plenty of trajectory analysis 14 CHAPTER 2. SIMULATION MODEL AND COMPUTATIONAL DETAILS tools that were developed at the research site.

There are 4 versions of the polymer model used in this work:

• The Chemically Realistic Chain (CRC) model is a modified version of the quantum chemistry-based united-atom model (from here on it is referred to as the original model) of 1,4-polybutadiene [START_REF] Smith | United atom force field for molecular dynamics simulations of 1,4polybutadiene based on quantum chemistry calculations on model molecules[END_REF] developed by G. D. Smith and W. Paul in 1998 for MD simulations of polymer melts. The original model was tuned to reproduce the structure and dynamics of bulk PBD very well by extensive comparison with experimental data [START_REF] Smith | Molecular dynamics of a 1,4-polybutadiene melt. comparison of experiment and simulation[END_REF][START_REF] Smith | NMR spin-lattice relaxation and conformational dynamics in a 1,4-polybutadiene melt[END_REF][START_REF] Smith | A molecular-dynamics simulation study of dielectric relaxation in a 1,4-polybutadiene melt[END_REF]. Later it was also used for simulations of confined films [START_REF] Yelash | Slow process in confined polymer melts: layer exchange dynamics at a polymer solid interface[END_REF][START_REF] Solar | The dielectric axation in polymer films: a comparison between experiments and atomistic simulations[END_REF][START_REF] Solar | Dielectric -relaxation of 1,4-polybutadiene confined between graphite walls[END_REF][START_REF] Solar | Relaxation processes and glass transition of confined polymer melts: a molecular dynamics simulation of 1,4-polybutadiene between graphite walls[END_REF][START_REF] Yelash | Three-step decay of time correlations at polymersolid interfaces[END_REF] using the GROMACS computer code [START_REF] Abraham | Gromacs: high performance molecular simulations through multi-level parallelism from laptops to supercomputers[END_REF]. However, due to the specifities of the simulations explained later in this chapter, this model cannot be used directly for the supported film geometry having free interfaces. Also, the direct transfer of some of the model features into LAMMPS is not possible. The CRC model developed by the author incorporates the modifications to account for these issues while maintaining the ability to reproduce the equation of state, the dynamics and the structure of polybutadiene well.

• The CRC2 and CRC4 models are derived from the CRC model by reducing the dihedral potentials by a factor of 2 and 4 correspondingly.

• The Freely Rotating Chain (FRC) model is also obtained from the CRC model by completely disabling the dihedral potentials.

Two types of polymer systems were studied in this work: a rather small system (∼ 5 × 10 3 united atoms) with periodic boundary conditions that models a bulk PBD and a quite large supported PBD film (∼ 10 nm thick with ∼ 10 5 united atoms) on top of a graphite substrate. All 4 models were used for the simulations of bulk PBD. In the case of films, only the CRC and the FRC models were employed, since such simulations are much more time and storage space consuming.

The following Section 2.2 presents the definition of the models, and Section 2.3

describes the simulation methodology. The parameters and comparison to the original CRC model are discussed in Section 2.4. The last Section 2.5 gives some information about technical details of processing the simulation trajectories. 116 UAs (only one chain length was studied in this work). For simplicity, no distinction between CH 2 groups within the chain and CH 3 groups at both ends of the chain is made.

Our PBD chains are random copolymers with 55% 𝑡𝑟𝑎𝑛𝑠 and 45% 𝑐𝑖𝑠 monomers.

The force field of PBD is composed of four potentials associated with bond stretching, bond angle bending, torsional rotations and pair interactions.

The stretching of the bond length ℓ is modeled by a harmonic potential

𝑈 bond (ℓ) = 𝐾 bond (ℓ -ℓ 0 ) 2 , (2.1)
where 𝐾 bond is (half) the force constant and ℓ 0 is the equilibrium bond length. There are CHAPTER 2. SIMULATION MODEL AND COMPUTATIONAL DETAILS Fig. 2.2 with the parameters given in Table 2.1 1 .

In the simulations the bond angle 𝜃 is calculated by the scalar product of the two consecutive bond vectors. Computationally, it is therefore convenient to express the bending potential in terms of cos 𝜃. We used the harmonic in cos 𝜃 bending potential 𝑈 bend (𝜃) = 𝐾 bend (cos 𝜃 -cos 𝜃 0 ) 2 , (

with the 𝐾 bend being (half) the bending stiffness constant and 𝜃 0 the equilibrium angle.

There are 2 types of bond angles in the system that correspond to different possible combinations of the consecutive pair of bond vectors along the PBD chain: CH 2 -CH-CH and CH-CH 2 -CH 2 (cf. Fig. 2.3 and Table 2.2). The scale of the potentials in Fig. 2.3 is quite big with respect to the value of the thermal energy 𝑘 B 𝑇 ≈ 0.7 kcal/mol at the highest studied temperature. Thus the fluctuations of the bond angle only happen in a small region around the equilibrium values. We will discuss this later in Section 2.4.2.

The torsion (or dihedral) angle 𝜙 is the angle between two intersecting planes defined by three successive bonds along the chain backbone. The associated potential 𝑈 tor (𝜙)

accounts for steric interactions between UAs separated by these three bonds. Since 𝑈 tor (𝜙)

is an even and periodic function of 𝜙 (0 °≤ 𝜙 ≤ 360°), it can be expressed as an expansion in powers of cos 𝜙. We employ the Ryckaert-Bellemans function [START_REF] Ryckaert | Molecular dynamics of liquid 𝑛-butane near its boiling point[END_REF] in our simulations 𝐴 𝑛 (cos 𝜙) 𝑛-1 .

(

There are 5 types of dihedral potentials in the system, that are depicted in Fig. 2.4

with the parameters {𝐴 𝑛 } 𝑛=1,…,6 given in Table 2.3. We have adopted the same naming as in Ref. [START_REF] Smith | United atom force field for molecular dynamics simulations of 1,4polybutadiene based on quantum chemistry calculations on model molecules[END_REF].

The double bond 𝑐𝑖𝑠 and 𝑡𝑟𝑎𝑛𝑠 potentials (Fig. 2.4b) determine the rotation along the double (CH-CH) bond axis, and therefore, the configuration of the PBD monomer.

Per definition, there is only one minimum, that occurs at 𝜙 min = 0°(or 360°) for the 𝑐𝑖𝑠 and at 𝜙 min = 180°for the 𝑡𝑟𝑎𝑛𝑠 potential.

The 𝛽 torsional potential determines the rotation along the axis of the CH 2 -CH 2 bond between two adjacent monomers. It has a shape familiar from hydrocarbon chains, such as polyethylene [2]. There are three minima, the primary being the 𝑡𝑟𝑎𝑛𝑠 state at 𝜙 = 180°and the secondary minima being the 𝑔𝑎𝑢𝑐ℎ𝑒-𝑝𝑙𝑢𝑠 and 𝑔𝑎𝑢𝑐ℎ𝑒-𝑚𝑖𝑛𝑢𝑠 states at 𝜙 ≃ 180 ∘ ± 120 ∘

The 𝛼 𝑐𝑖𝑠 and 𝛼 𝑡𝑟𝑎𝑛𝑠 potentials correspond to the CH-CH 2 bond that is next to the double bond. Qualitatively, the shape of 𝑈 tor (𝜙) for these dihedrals is the mirror image of the 𝛽 potential, having minima at 𝜙 = 0 ∘ , 360 ∘ and 𝜙 ≃ 180 ∘ ± 60 ∘ .

The 𝑐𝑖𝑠 and 𝑡𝑟𝑎𝑛𝑠 potentials (both double bond and 𝛼) are randomly assigned to the monomers in a PBD chain with the proportion of 45% 𝑐𝑖𝑠 and 55% 𝑡𝑟𝑎𝑛𝑠.

The values of the parameters {𝐴 𝑛 } 𝑛=1,…,6 in Table 2.3 correspond to our chemically with 𝜆𝑈 tor (𝜙) where 𝜆 = 0, 1/4, 1/2, 1. We refer to these models as:

𝜆 = 0 ∶ freely rotating chain (FRC) model , 𝜆 = 1 4 ∶ CRC4 model , 𝜆 = 1 2 ∶ CRC2 model , (2.4
)

𝜆 = 1 ∶ CRC model .
The rescaled (𝜆 = 1 2 , 1 4 ) or disabled (𝜆 = 0) torsional potentials is the only difference between the models. The FRC model was already introduced and studied in Refs. [START_REF] Krushev | The role of internal rotational barriers in polymer melt chain dynamics[END_REF][START_REF] Krushev | Intramolecular caging in polybutadiene due to rotational barriers[END_REF].

Due to the absence of torsional potential the model enables free rotation around the double bonds, thereby eliminating the difference between 𝑐𝑖𝑠 and 𝑡𝑟𝑎𝑛𝑠 conformers. Within the FRC model PBD is therefore no longer a random copolymer of 𝑐𝑖𝑠 and 𝑡𝑟𝑎𝑛𝑠 monomers, but a regular copolymer of CH and CH 2 units. While it would have been possible to preserve the random 𝑐𝑖𝑠/𝑡𝑟𝑎𝑛𝑠 copolymer character by maintaining the dihedral potentials for the double bonds only [START_REF] Smith | Relationship between the 𝛼-and 𝛽-relaxation processes in amorphous polymers: insight from atomistic molecular dynamics simulations of 1,4 polybutadiene melts and blends[END_REF], we have adopted the FRC model here because it resembles extensively studied glass-forming bead-spring models with intramolecular forces resulting only from bond-length and bending potentials 2 [START_REF] Colmenero | Are polymers standard glass-forming systems? The role of intramolecular barriers on the glass-transition phenomena of glass-forming polymers[END_REF][START_REF] Bernabei | Dynamic arrest in polymer melts: competition between packing and intramolecular barriers[END_REF][START_REF] Bernabei | The role of intramolecular barriers on the glass transition of polymers: computer simulations versus mode-coupling theory[END_REF][START_REF] Bernabei | From caging to rouse dynamics in polymer melts with intramolecular barriers: a critical test of the mode coupling theory[END_REF][START_REF] Schnell | Simulated glass-forming polymer melts: glass transition temperature and elastic constants of the glassy state[END_REF]. The CRC4 and CRC2 models are introduced because they interpolate between the FRC and CRC models.

Pair interactions between united atoms of a chain separated by four or more bonds and between united atoms of different chains are modeled by a 12-6 Lennard-Jones (LJ) potential with an additional switching function 𝑆(𝑟) commonly used in the GROMACS code [START_REF] Gromacs | GROningen MAchine for Chemical Simulation[END_REF], but also implemented in LAMMPS (see pair_style lj/gromacs command [START_REF] Lammps | Large-scale Atomic/Molecular Massively Parallel Simulator[END_REF]).

The potential reads:

𝑈 pair (𝑟) = ⎧ { ⎨ { ⎩ 4𝜖 [( 𝜎 𝑟 ) 12 -( 𝜎 𝑟 ) 6 ] + 𝑆(𝑟) if 𝑟 < 𝑟 c , 0 if 𝑟 ≥ 𝑟 c , (2.5)
with

𝑆(𝑟) = ⎧ { ⎨ { ⎩ 𝐶 if 𝑟 < 𝑟 1 , 𝐴 3 (𝑟 -𝑟 1 ) 3 + 𝐵 4 (𝑟 -𝑟 1 ) 4 + 𝐶 if 𝑟 1 ≤ 𝑟 < 𝑟 c , (2.6)
2 Though, the bond length of the FRC model is much shorter than in the bead-spring models.

20

CHAPTER 2. SIMULATION MODEL AND COMPUTATIONAL DETAILS and

𝐴 = -3𝑈 ′ term (𝑟 c ) + (𝑟 c -𝑟 1 )𝑈 ′′ term (𝑟 c ) (𝑟 c -𝑟 1 ) 2 , 𝐵 = 2𝑈 ′ term (𝑟 c ) + (𝑟 c -𝑟 1 )𝑈 ′′ term (𝑟 c ) (𝑟 c -𝑟 1 ) 3 , 𝐶 = -𝑈 term (𝑟 c ) + 1 2 (𝑟 c -𝑟 1 )𝑈 ′ term (𝑟 c ) - 1 12 (𝑟 c -𝑟 1 ) 2 𝑈 ′′ term (𝑟 c ) .
(2.7)

where 𝑈 term (𝑟 c ) stands either for the repulsive term, 4𝜖(𝜎/𝑟) 12 , or for the attractive term, -4𝜖(𝜎/𝑟) 6 , of the LJ potential, 𝜎 is the zero-crossing distance (𝑈 pair (𝑟 = 𝜎) = 0), 𝜖 is the value of the minimum of the potential at 𝑈 pair (𝑟 = 2

1
6 𝜎) = -𝜖. Prime and double prime symbols denote respectively the first and second derivatives with respect to 𝑟. The switching function obeys the following boundary conditions:

𝑆 ′ (𝑟 1 ) = 𝑆 ′′ (𝑟 1 ) = 0 , 𝑆(𝑟 𝑐 ) = -𝑈 term (𝑟 c ) , 𝑆 ′ (𝑟 𝑐 ) = -𝑈 ′ term (𝑟 c ) , 𝑆 ′′ (𝑟 𝑐 ) = -𝑈 ′′ term (𝑟 c ) . (2.8)
For 𝑟 < 𝑟 1 Eqs. (2.5) and (2.6) imply that 𝑈 pair (𝑟) is just a shifted LJ potential, while for 𝑟 1 ≤ 𝑟 < 𝑟 c the LJ potential and its 1 st and 2 nd derivatives are progressively smoothed to zero due to the form of 𝑆(𝑟) and the boundary conditions.

We choose identical values, 𝑟 1 = 9 Å and 𝑟 c = 12 Å, for all pairs of UAs. The LJ parameters used in the present study are summarized in Table 2.4 and the corresponding potentials are plotted in Fig. 2.5

In the case of supported film, the united atoms of PBD interact with a substrate, that is modeled using the crystalline lattice of graphite atoms with mass 𝑚 CH 2 = 12.011 u.

The LJ parameters for the graphite were taken from the work [START_REF] Yelash | Slow process in confined polymer melts: layer exchange dynamics at a polymer solid interface[END_REF], which in turn has taken the values reported in [START_REF] Steele | The physical interaction of gases with crystalline solids: i. gas-solid energies and properties of isolated adsorbed atoms[END_REF]. They are also given in Table 2.4 and Fig. 2.5.

Only from the pair interactions point of view PBD corresponds to a binary mixture (see Section 3.3). The values in for better glass-forming ability [START_REF] Ninarello | Models and algorithms for the next generation of glass transition studies[END_REF][START_REF] Ingebrigtsen | Crystallization instability in glass-forming mixtures[END_REF]. In the case of the interaction of PBD with the graphite, the Lorentz-Berthelot combining rules were used.

Our force field is different from the one used in the original works [START_REF] Smith | United atom force field for molecular dynamics simulations of 1,4polybutadiene based on quantum chemistry calculations on model molecules[END_REF][START_REF] Yelash | Slow process in confined polymer melts: layer exchange dynamics at a polymer solid interface[END_REF][START_REF] Solar | The dielectric axation in polymer films: a comparison between experiments and atomistic simulations[END_REF][START_REF] Solar | Dielectric -relaxation of 1,4-polybutadiene confined between graphite walls[END_REF][START_REF] Solar | Relaxation processes and glass transition of confined polymer melts: a molecular dynamics simulation of 1,4-polybutadiene between graphite walls[END_REF][START_REF] Smith | Molecular dynamics of a 1,4-polybutadiene melt. comparison of experiment and simulation[END_REF][START_REF] Smith | NMR spin-lattice relaxation and conformational dynamics in a 1,4-polybutadiene melt[END_REF][START_REF] Smith | A molecular-dynamics simulation study of dielectric relaxation in a 1,4-polybutadiene melt[END_REF][START_REF] Yelash | Three-step decay of time correlations at polymersolid interfaces[END_REF], that used the GROMACS code, in two aspects. In these works the bonds were constrained and the plain LJ potential truncated at 𝑟 c = 9 Å with long-range tail corrections to energy and pressure was used for pair interactions. The reasons for these differences are discussed in Section 2.4.

The simulation methodology 2.3.1 General information

We perform molecular dynamics (MD) simulations of PBD with the Large Scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [START_REF] Allen | Computer simulation of liquids[END_REF][START_REF] Lammps | Large-scale Atomic/Molecular Massively Parallel Simulator[END_REF][START_REF] Plimpton | Fast parallel algorithms for short-range molecular dynamics[END_REF] pressure equal to 𝑝 = 0 atm was imposed under the NPT conditions and the corresponding volume 𝑉 was imposed under the NVT conditions. In this way, the melt can be compared to the simulations of the supported films in contact with vacuum at the free surface.

The initial configurations at 𝑇 = 353 K for the bulk and the film systems were taken from the work [START_REF] Solar | Relaxation processes and glass transition of confined polymer melts: a molecular dynamics simulation of 1,4-polybutadiene between graphite walls[END_REF], which employed the implementation of the realistic model (with all dihedrals) in GROMACS code [START_REF] Gromacs | GROningen MAchine for Chemical Simulation[END_REF]. These configurations were converted to LAMMPS input files by means of a custom-made script using a conversion 1 kcal/mol = 4184 J/mol for the energy units. The bulk system was used to establish the parameters of our models in LAMMPS and compare to the original GROMACS simulations (see Sections 2.3.2 and 2.4). The film systems were simulated using our models from the beginning. We

Bulk details

The bulk system contains 𝑁 c = 40 PBD chains in a cubic box of about 50 Åin size (Fig. 2.6). Periodic boundary conditions were applied in all directions. The initial configuration at 𝑇 = 353 K was converted from GROMACS to LAMMPS format by means of a custom-made script, that properly converts the units and replaces constrained bonds by the harmonic potential.

In order to verify if the model was transferred from GROMACS to LAMMPS correctly, our first simulations were conducted using the truncated LJ potential and tail corrections used in both codes. This allowed us to make a comparison to the original GROMACS simulations, which used such a potential in the first place. From now on, we refer to this model as the 'original model'. to compare with the works [START_REF] Krushev | The role of internal rotational barriers in polymer melt chain dynamics[END_REF][START_REF] Krushev | Intramolecular caging in polybutadiene due to rotational barriers[END_REF], two working temperatures were chosen: 𝑇 = 353 K and 𝑇 = 273 K and the equilibration time was set to 100 ns. At each working temperature, the following simulation protocol was adopted (Fig. 2.7).

• Starting from the initial configuration, the equilibration is performed under the NPT • In the subsequent deformation stage we start from the instantaneous volume of the final configuration of the NPT run and isotropically change the size of the simulation box until the equilibrium volume 𝑉 (𝑇 ) is reached (fix deform command in LAMMPS [START_REF] Lammps | Large-scale Atomic/Molecular Massively Parallel Simulator[END_REF]). The deformation is carried out over a short time of 10 ps.

• Due to the fast deformation, it is necessary to relax residual stresses possibly created by the quenching of the box volume, even though the deformation is small (≲ 0.1 Å). This is done over another 100 ns of the NVT run at fixed 𝑉 (𝑇 ), which suffices to stabilize the average pressure at 𝑝 = 0 atm.

• After equilibration we perform a NVT production run that lasts at least 1 𝜇s (10 9 timesteps).

At 𝑇 = 353 K the initial (𝑡 = 0) configuration is the one converted from the work [START_REF] Solar | Relaxation processes and glass transition of confined polymer melts: a molecular dynamics simulation of 1,4-polybutadiene between graphite walls[END_REF].

In order to obtain the initial configuration at 𝑇 = 273 K, the system was continuously cooled down to 𝑇 = 113 K at imposed pressure 𝑝 = 0 atm according to the ramp

𝑇 (𝑡) = 𝑇 0 -Γ𝑡 , (2.9)
CHAPTER 2. SIMULATION MODEL AND COMPUTATIONAL DETAILS with the cooling rate Γ = 20 K / 6 ns = 3.33 K/ns (Fig. 2.7). Then the configuration that corresponded to the time point with 𝑇 = 273 K was used as the starting one (𝑡 = 0) for constant temperature NPT equilibration.

In order to check the density of PBD, the NPT equilibration at 𝑇 = 353 K was continued up to 1 𝜇s and the configurations at 𝑡 = 200 ns, 400 ns, 600 ns, 800 ns, 1 𝜇𝑠 were used to start 9 additional cooling runs: 4 runs with Γ = 3.33 K/ns and 5 runs with Γ = 0.83 K/ns. The 𝑡 = 200 ns point was also used to launch a cooling run of the system where all dihedral potentials were disabled 3 . This is motivated by the results of Honnell 𝑒𝑡 𝑎𝑙. [START_REF] Honnell | On the pressure equation for chain molecules[END_REF], who derived that the pressure equation for chain molecules does not explicitly depend on the angular potentials, and the works [START_REF] Krushev | The role of internal rotational barriers in polymer melt chain dynamics[END_REF][START_REF] Krushev | Intramolecular caging in polybutadiene due to rotational barriers[END_REF], where the dynamics without the dihedrals was reported to be orders of magnitude faster than with the dihedrals while the structure is identical. Thus, we expect the density at 𝑝 = 0 atm to be indistinguishable between the realistic model and the model without the dihedrals within the fluctuations at studied temperatures. We also use this argument later for the initial configurations of our CRC2, CRC4 and FRC models.

From one of the slower cooling runs, we then followed the protocol described above of supported films and developed the modified pair potential to yield equivalent result in inhomogeneous systems (free surface!) without tail correction (see Section 2.4.3 for detailed explanation). We used the system at 𝑇 = 273 K to launch the runs with modified pair interactions to develop our CRC model. The aim of the modifications was to match the densities of the systems at the studied temperatures without having to resort to the tail corrections. The new pair potential was verified to recover desired density by simulating a system at 𝑇 = 353 K and performing the same set of cooling runs using our CRC model.

Due to the indistinguishable densities and for time saving purposes, we used the equilibrium volumes 𝑉 (𝑇 ) obtained using the original model to continue our simulation protocol and performed the NVT equilibrations and the production runs using already our CRC/CRC2/CRC4/FRC models. As a last check, we also performed long (1 𝜇s) NPT runs at each temperature using all four models to make sure that the densities are indeed the same.

Additionally, the system at 𝑇 = 213 K was simulated using all our models in a similar manner as the higher temperatures, with an exception that the NPT equilibration has been performed using our CRC model. In addition, the systems at 𝑇 = 153, 113 K were simulated using the FRC model following the protocol described above..

Film details

The initial film configuration was taken from the work [START_REF] Solar | Relaxation processes and glass transition of confined polymer melts: a molecular dynamics simulation of 1,4-polybutadiene between graphite walls[END_REF], where the film comprising 720 PBD chains (83520 UAs) was confined between two layers of graphite atoms at the bottom of the simulation box and one layer at the top with periodic boundary conditions in all directions. In order to make the film supported on one side with free surface on the other side, we have added about 10 nm of vacuum above the top graphite layer, then moved the whole system of PBD and graphite up by 3.4 Å (the distance between the graphite layers). Finally the top layer of graphite was moved to the bottom of the simulation box, so all three layers are together and separated by 3.4 Å. A snapshot of the box is presented in Fig. 2.8.

All film simulations are under NVT conditions. The presence of the free surface of the film in contact with vacuum effectively corresponds to the simulations at zero pressure.

The periodic boundary conditions were used in all directions. In the lateral dimensions the simulations box is commensurable with the (0,0,1) graphite crystalline surface. During the simulations all graphite atoms are fixed at their equilibrium positions in the crystal lattice. In the 𝑧-dimension, the 10 nm of vacuum is sufficient for the united atoms of PBD not to interact with their periodic images.

The simulations of the films were performed using our CRC and FRC models. During the conversion of the initial configuration, the pair interactions were replaced by the ones from our models. The film simulation protocol is similar to the bulk one and is presented schematically in Fig. 2.9

• The initial configuration at 𝑇 = 353 K is equilibrated over 100 ns.

• Then the cooling run is performed according to the ramp given by Eq. (2.9) with the cooling rate Γ = 0.83 K/ns.

• The cooling configurations that corresponded to the time points of the working temperatures 𝑇 = 273, 253, 240, 225, 213 K were used as starting ones for the constant temperature equilibration and production runs.

In fact, there is one long constant temperature run at each working temperature.

The equilibration and production parts are distinguished simply by the time point, from which we start measuring the quantities of interest. The same is valid for the systems at 𝑇 = 353 K: the configuration after 100 ns of equilibration using the CRC model was copied to launch the cooling (and also the FRC equilibration) and the run was continued. The constant temperature runs lasted 1 𝜇s in total in the case of the CRC model. The FRC runs at the same working temperatures lasted 700 ns. As in the bulk, the two additional temperatures 𝑇 = 153 K and 𝑇 = 113 K were performed using the FRC model, that lasted also 1 𝜇s. The detailed map of our simulations is given in Fig.

A.2 in Appendix A.

Model parameters and previous works

General information

As it was mentioned previously, our CRC model has several differences in the force field with respect to the model used in recent works [START_REF] Yelash | Slow process in confined polymer melts: layer exchange dynamics at a polymer solid interface[END_REF][START_REF] Solar | The dielectric axation in polymer films: a comparison between experiments and atomistic simulations[END_REF][START_REF] Solar | Polymer dynamics in a polymer-solid interphase: molecular dynamics simulations of 1,4-polybutadiene at a graphite surface[END_REF][START_REF] Solar | Dielectric -relaxation of 1,4-polybutadiene confined between graphite walls[END_REF][START_REF] Solar | Relaxation processes and glass transition of confined polymer melts: a molecular dynamics simulation of 1,4-polybutadiene between graphite walls[END_REF][START_REF] Yelash | Three-step decay of time correlations at polymersolid interfaces[END_REF]. These works were focused on polymer films and, in turn, also employed an adapted version of the original potential functions [START_REF] Smith | United atom force field for molecular dynamics simulations of 1,4polybutadiene based on quantum chemistry calculations on model molecules[END_REF][START_REF] Smith | Molecular dynamics of a 1,4-polybutadiene melt. comparison of experiment and simulation[END_REF] for implementation in the GROMACS code. Since our study is also focused on films, we prioritized the recent version of the model [START_REF] Yelash | Slow process in confined polymer melts: layer exchange dynamics at a polymer solid interface[END_REF][START_REF] Solar | The dielectric axation in polymer films: a comparison between experiments and atomistic simulations[END_REF][START_REF] Solar | Polymer dynamics in a polymer-solid interphase: molecular dynamics simulations of 1,4-polybutadiene at a graphite surface[END_REF][START_REF] Solar | Dielectric -relaxation of 1,4-polybutadiene confined between graphite walls[END_REF][START_REF] Solar | Relaxation processes and glass transition of confined polymer melts: a molecular dynamics simulation of 1,4-polybutadiene between graphite walls[END_REF][START_REF] Yelash | Three-step decay of time correlations at polymersolid interfaces[END_REF] when implementing our simulations in LAMMPS for both bulk and film systems. We Figure 2.9: Sketch of the film simulation protocol consisting of a cooling run according to Eq. (2.9) followed by the equilibration and production runs.

CHAPTER 2. SIMULATION MODEL AND COMPUTATIONAL DETAILS

Intra-chain interactions

In order to verify if the intra-chain interactions were transferred from GROMACS to LAMMPS correctly, we set up rather short simulations of the bulk PBD using the original GROMACS scripts obtained from the authors of the work [START_REF] Solar | Relaxation processes and glass transition of confined polymer melts: a molecular dynamics simulation of 1,4-polybutadiene between graphite walls[END_REF] using the same initial configuration and compare them to our LAMMPS simulations. Following the results of [START_REF] Shirts | Lessons learned from comparing molecular dynamics engines on the sampl5 dataset[END_REF], our objective was not to obtain an exact match of the energies and particle movements during both simulations (it is impossible due to the difference in the software features discussed later), but rather to get a qualitative agreement of the mean values of the energy terms and a quantitative agreement of the representative structural and dynamical properties of PBD. At the moment of performing the comparison the pair interactions had not been modified yet. The truncated LJ potential with tail corrections was used with the same parameters in both MD engines.

The main difference in the intra-chain interactions between our LAMMPS simulations and the original GROMACS simulations is in the bond potential. All bonds between united atoms of PBD were constrained to their equilibrium values ℓ 0 in GROMACS.

LAMMPS does not allow to constrain the backbone bonds of a chain in order to achieve good parallel performance. However, we wanted to stay close to the original simulations with ℓ = ℓ 0 , therefore we used a harmonic bond potential instead with a large value for the force constant 𝐾 bond = 2500 kcal/(mol Å 2 (cf. Eq. (2.1) and Table 2.1). This value is a compromise between spectroscopic data [START_REF] Nallasamy | Vibrational spectra and assignments of cis-and trans-1,4-polybutadiene[END_REF] and not too long computation times.

With this value one finds √ ⟨(ℓ -ℓ 0 ) 2 ⟩ = √𝑘 B 𝑇 /(2𝐾 bond ) ≈ 0.01 Å for 𝑇 = 353 K, the highest studied temperature. Thus, bond length fluctuations are less then 1% of ℓ 0 . The estimation of the smallest bond oscillation time in the system as the period of oscillation of a bond between two CH united atoms (the atoms with the smallest mass) results in and 𝜈 C=C ≈ 1650 cm -1 [START_REF] Nallasamy | Vibrational spectra and assignments of cis-and trans-1,4-polybutadiene[END_REF]). For the latter, the peaks at ℓ 0 = 1.5 Å and ℓ 0 = 1.53 Å are merged into one peak in between. Even though this is probably closer to the experiment, our intention at this point was to introduce as little deviation to the model as possible and to reduce the influence on the other potentials for better comparison with the constrained system. We have compared the distribution measured in the simulations to the Boltzmann distribution

𝜏 bond = 2𝜋√𝑚 CH /(2𝐾 bond ) ≈
𝑝 bond (ℓ) = exp(-𝑈 bond (ℓ)/𝑘 B 𝑇 ) ∫ ∞ 0 dℓ exp(-𝑈 bond (ℓ)/𝑘 B 𝑇 ) (2.10)
and obtained an excellent agreement, which was expected because the bond potential is the strongest in the system.

As for the bending angle potentials, our LAMMPS simulations and recent GROMACS works used a harmonic in cos 𝜃 bending potential given by Eq. (2.2) and In this way, 𝐾 bend was determined from 𝑘 𝜃 of Ref. [START_REF] Smith | United atom force field for molecular dynamics simulations of 1,4polybutadiene based on quantum chemistry calculations on model molecules[END_REF]. 𝐾 bend is smaller than 𝐾 bond by a factor of about 37, implying larger fluctuations of 𝜃 than for ℓ. Still, the fluctuations remain very small because √⟨(𝜃 -𝜃 0 ) 2 ⟩ ≈ [𝑘 B 𝑇 /(2𝐾 bend sin 2 𝜃 0 )] 1/2 ≈ 5 ∘ even for 𝑇 = 353 K. (2.13)

The data agree very well with the latter verifying the correct transfer of the bending angle potential from GROMACS to LAMMPS. The bending potential is also a hard variable in the system.

The torsion potential has been transferred from GROMACS to LAMMPS exactly. In both cases it is given by Eq. (2.3) and Table 2.3, while the initially developed model [START_REF] Smith | United atom force field for molecular dynamics simulations of 1,4polybutadiene based on quantum chemistry calculations on model molecules[END_REF][START_REF] Smith | Molecular dynamics of a 1,4-polybutadiene melt. comparison of experiment and simulation[END_REF] utilized

𝑈 tor (𝜙) = 1 2 6 ∑ 𝑛=1 𝑘 𝑛 [1 -cos(𝑛𝜙)] . (2.14)
The coefficients {A_n} in (2.15)

The LAMMPS and GROMACS distributions agree with each other very well indicating a correct implementation of the dihedrals in our model. Both of them, however, differ from the Boltzmann distributions: the dihedral angles are redistributed from 0 ∘ and 360 ∘ to 120 ∘ and 240 ∘ . We examine this in detail in Section 3.1.1. According to the work [START_REF] Smith | United atom force field for molecular dynamics simulations of 1,4polybutadiene based on quantum chemistry calculations on model molecules[END_REF], these deviations come from the influence of the intra-chain pair interactions between united atoms separated by four and more bonds on the 𝛼 𝑐𝑖𝑠 dihedral. In fact, in Ref. [START_REF] Smith | United atom force field for molecular dynamics simulations of 1,4polybutadiene based on quantum chemistry calculations on model molecules[END_REF] the potentials were optimized to account for an effective torsional potential created by the presence of the pair interactions and to reproduce the correlations between the consecutive torsions.

As the last check of the model implementation we compared the structure and the dynamics of the PBD obtained from our LAMMPS simulations to the constrained GROMACS simulations. As a representative structural quantity we used the melt and single chain structure factors 𝑆(𝑞) and 𝐹 (𝑞) (Fig. 2.13):

𝑆(𝑞) = 1 𝑁 c 𝑁 𝑁 c ∑ 𝑖,𝑗=1 𝑁 ∑ 𝑎,𝑏=1 ⟨exp[-𝑖 ⃗ 𝑞( ⃗ 𝑟 𝑎,𝑖 -⃗ 𝑟 𝑏,𝑗 )]⟩ , (2.16
)

𝐹 (𝑞) = 1 𝑁 c 𝑁 𝑁 c ∑ 𝑖=1 𝑁 ∑ 𝑎,𝑏=1 ⟨exp[-𝑖 ⃗ 𝑞( ⃗ 𝑟 𝑎,𝑖 -⃗ 𝑟 𝑏,𝑖 )]⟩ , (2.17)
where 𝑁 c is the number of chains, 𝑁 the number of atoms in a chain, ⃗ 𝑞 the wavevector, ⃗ 𝑟 𝑎,𝑖 the position of atom 𝑎 of the chain 𝑖. For the 𝑆(𝑞) the discrete ⃗ 𝑞 are used, which are commensurable with the simulation box (the components of ⃗ 𝑞 are multiples of 2𝜋/𝐿 with 𝐿 being the box size). The dynamics was compared using the mean-square displacement of united atoms 𝑔 0 (𝑡) (Fig. 2.14)

𝑔 0 (𝑡) = 1 𝑁 c 𝑁 𝑁 c 𝑁 ∑ 𝑛=1 ⟨[ ⃗ 𝑟 𝑛 (𝑡) -⃗ 𝑟 𝑛 (0)] 2 ⟩, (2.18)
where ⃗ 𝑟 𝑛 (𝑡) the position of the 𝑛-th united atom at time 𝑡.

We have also compared our results to the data published in Ref. [START_REF] Krushev | Intramolecular caging in polybutadiene due to rotational barriers[END_REF], which is an older work than Ref. [START_REF] Solar | Relaxation processes and glass transition of confined polymer melts: a molecular dynamics simulation of 1,4-polybutadiene between graphite walls[END_REF], from which we received the initial configurations and GROMACS scripts. The data are in a very good agreement at the studied temperatures.

Pair interactions

As it was already mentioned in the previous chapters, the originally developed simulations of bulk PBD [START_REF] Smith | United atom force field for molecular dynamics simulations of 1,4polybutadiene based on quantum chemistry calculations on model molecules[END_REF][START_REF] Smith | Molecular dynamics of a 1,4-polybutadiene melt. comparison of experiment and simulation[END_REF] and the recent works on the confined films [START_REF] Solar | The dielectric axation in polymer films: a comparison between experiments and atomistic simulations[END_REF][START_REF] Solar | Dielectric -relaxation of 1,4-polybutadiene confined between graphite walls[END_REF][START_REF] Solar | Relaxation processes and glass transition of confined polymer melts: a molecular dynamics simulation of 1,4-polybutadiene between graphite walls[END_REF] utilized the plain Lennard-Jones potential truncated at 𝑟 c = 9 Å for pair interactions

𝑈 pair (𝑟) = ⎧ { ⎨ { ⎩ 4𝜖 [( 𝜎 𝑟 ) 12 -( 𝜎 𝑟 ) 6 ] if 𝑟 < 𝑟 c , 0 if 𝑟 ≥ 𝑟 c , (2.19)
with the parameters for different types of united atoms given in Table 2.5. To account for missing interactions for distances 𝑟 > 9 Å, the long range tail corrections [START_REF] Allen | Computer simulation of liquids[END_REF][START_REF] Sun | Compass: an ab initio force-field optimized for condensed-phase applicationsoverview with details on alkane and benzene compounds[END_REF][START_REF] Frenkel | Understanding molecular simulation[END_REF] to pair energy and pressure were applied. We have also used the truncated LJ potential with tail corrections when setting up our first simulation of bulk PBD in LAMMPS. In practice, Ref. [START_REF] Krushev | Intramolecular caging in polybutadiene due to rotational barriers[END_REF].

these corrections are constant values calculated in the beginning of the simulation 7 , which are added to the values of pair energy and pressure at each timestep. The corrections to the energy, in fact, do not affect the dynamics of the system, because the forces between the particles are proportional to the gradient of the potential, which eliminates the constant.

The pressure correction, however, has a direct influence on the dynamics under the NPT conditions, because it enters the barostating routine. This modifies the volume of the simulation box, and therefore, the density according to the desired pressure. For the bulk 7 see Eq. ( 5) of Ref. [START_REF] Sun | Compass: an ab initio force-field optimized for condensed-phase applicationsoverview with details on alkane and benzene compounds[END_REF] for the analytical expressions of the corrections system under the NVT conditions, the influence of the corrections comes only via the value of the fixed volume. Enabling/disabling the tail corrections in the bulk NVT simulations at the same fixed volume does not change the dynamics of the simulated system whatsoever.

The situation is different for the simulations of the films, which are NVT simulations only. In the confined geometry the density can be adjusted by varying the distance between the walls. The works [START_REF] Solar | The dielectric axation in polymer films: a comparison between experiments and atomistic simulations[END_REF][START_REF] Solar | Dielectric -relaxation of 1,4-polybutadiene confined between graphite walls[END_REF][START_REF] Solar | Relaxation processes and glass transition of confined polymer melts: a molecular dynamics simulation of 1,4-polybutadiene between graphite walls[END_REF] used this method to match the density of PBD in the middle of the film to the bulk value. In the supported films, however, the density is determined by the interactions, and since the tail corrections have no effect in the NVT simulations, the value in the middle of the film corresponds to the one of the bulk simulated under the NPT conditions without the corrections. This issue does not allow us to continue using the truncated LJ potential with tail corrections, because we would not be able to directly compare the bulk and the film simulations under the same conditions.

Yet, the original model has been validated against experiments [START_REF] Smith | Molecular dynamics of a 1,4-polybutadiene melt. comparison of experiment and simulation[END_REF][START_REF] Smith | NMR spin-lattice relaxation and conformational dynamics in a 1,4-polybutadiene melt[END_REF][START_REF] Smith | A molecular-dynamics simulation study of dielectric relaxation in a 1,4-polybutadiene melt[END_REF][START_REF] Byutner | Prediction of the linear viscoelastic shear modulus of an entangled polybutadiene melt from simulation and theory[END_REF][START_REF] Byutner | Viscoelastic properties of polybutadiene in the glassy regime from molecular dynamic simulations[END_REF], a salient feature we want to keep. Fig. 2.15a presents the bulk density of PBD obtained from the cooling runs at zero pressure using our LAMMPS simulations with the truncated LJ potential and corrections enabled. We compared our results to various experimental and simulation references, that we were able to find (see Appendix C for an overview).

The data in the literature are for PBD of various compositions and molecular weight under different pressure conditions. Thus, it only gives a range of densities to compare with.

Due to the presence of the tail corrections, our data lay nicely within the references.

Disabling them leads to a very significant decrease of the density in the NPT simulations at 𝑇 = 353 K due to quite small 𝑟 c = 9 Å, as it is shown in Fig. 2.15b. Increasing 𝑟 c up to 14 Å results in densities closer to the experimental ones, but also leads to more than twice longer CPU times.

In order to avoid the tail corrections completely, we decided to apply the switching function 𝑆(𝑟) to the LJ potential, as it is described in Section 2.2 (cf. Eqs. (2.5) to (2.7)).

The inner cutoff in Eq. (2.5) was kept at 𝑟 1 = 9 Å, so that the shape of the potential at 𝑟 < 𝑟 1 is kept as in the original model. The outer cutoff was set at 𝑟 c = 12 Å. The switching function 𝑆(𝑟) smoothes the potential and its 1 st and 2 nd derivatives to zero over a distance of about one particle diameter 9 Å < 𝑟 < 12 Å.

The modified pair potential Eq. (2.5) is of finite range and is still of shorter range than the original model with tail corrections. This difference affects the density of the CHAPTER 2. SIMULATION MODEL AND COMPUTATIONAL DETAILS melt, which is smaller for the modified model than for the original model and moreover does not fit into the interval of reference data from the literature (Fig. 2.15). Therefore, it is necessary to adjust the potential parameters in order to reproduce the desired density.

To this end, we choose to leave the zero-crossing distance 𝜎 unchanged but to enhance 𝜖 by a factor 𝜆 > 1. This choice is motivated by the generally accepted point of view [START_REF] Hansen | Theory of simple liquids[END_REF] that the local structure of (simple) liquids at high density is dominated by the packing of the molecular repulsive cores (i.e. by 𝜎), whereas the attractive interactions (of strength 𝜖)

determine the cohesive energy and along with that, thermodynamic properties, such as the density or compressibility.

To determine 𝜆 we proceed in the following way. For the original model with truncated LJ potential and tail corrections we carry out an NPT simulation to get the reference value of the density at pressure 𝑝 = 0 atm for two temperatures, 𝑇 = 273 K and 353 K.

As described in Section 2.3. 2.5 (they are the same as in Table 2.4).

To verify the estimated 𝜆 we performed NVT simulations at 𝑇 = 273 K and at 𝑇 = 373 K with the new values of 𝜖 𝑖𝑗 . For both temperatures the densities and pressures agree with the target values of the original model.

We have also performed an inverse cross-check by imposing the 𝑝 = 0 atm and comparing the density of PBD obtained from the cooling runs at two rates Γ = 3.33 K/ns and Γ = 0.83 K/ns using the original and the modified potential. The results are presented in Fig. 2.16.

All curves superimpose well at 𝑇 > 225 K, so our new model is able to represent the density of PBD. Depending on the cooling rate Γ, the curves bend around 𝑇 ≈ 200 K, which suggests a glass transition. The modified model appears to bend slightly earlier upon cooling the system, which is expected due to the increased 𝜖.

In addition, we have compared our results to the density obtained from the first simulations where all the dihedrals had been disabled, which also superimposed well in the The density of PBD obtained from the cooling runs using the original and the modified pair interactions. The purple line represents the cooling at the rate Γ = 3.33 K/ns using the original model. This is the same data as in the Fig. 2.15. The green line is the data obtained from a slower cooling rate of Γ = 0.83 K/ns. The black and the red lines are from the simulations using the modified model. The data described above were obtained by averaging over 5 independent cooling runs (the simulation details are given in Appendix A). The blue line corresponds to a single cooling run using the original model without the dihedrals. region of the working temperatures. This model has been already studied in Refs. [START_REF] Krushev | The role of internal rotational barriers in polymer melt chain dynamics[END_REF][START_REF] Krushev | Intramolecular caging in polybutadiene due to rotational barriers[END_REF], in which much faster dynamics of the PBD was reported for the system without dihedrals while the structure was observed to be unchanged. Thus, even though our cooling run without dihedrals starts from the configuration obtained from the system with dihedrals, we expect all perturbations that may occur by doing so (if any) to be relaxed during first few nanoseconds of the run. This allowed us to omit the NPT equilibration step for the models with modified dihedrals, and to use the densities obtained using our CRC model for the NVT simulations of our models with reduced and disabled dihedral barriers. We will continue this discussion in Section 3.2 where we analyse the cooling runs using all our models.

Details of sampling and analysis of the trajectories

In order to compute the quantities of interest there are trajectory files written during the simulation, in which the coordinates of all particles in the simulated system (usually called "a frame") are saved at selected timesteps. Our simulations write two types of such files. The first type contains frames written at timesteps that are equidistantly spaced (Fig. 2.17a). These files are used for calculating ensemble averages of the static properties.

The sequence of frames is firstly divided into equal blocks, and the quantities are computed and averaged within each block separately. The average of the results over all blocks is then performed along with the standard deviation of mean, that is usually used as the estimation of error bars. The second type of files, which we call "logarithmic" trajectory, also has blocks of frames, but the frames are spaced within each block exponentially at short times, such that there is an equal distribution of points when plotting dynamic quantities logarithmically in time (Fig. 2.17b). These files are used for calculating the time averages with the resolution of the time windows varying over several decades: these "logarithmic blocks" produce the short time average, and the first frame of each block is used for the long time average as in the files of the first type.

The thermodynamic quantites in the simulations were calculated by LAMMPS itself and writted into separate log files. The various static and dynamic quantities presented in this work were calculated using trajectory analysis tools developed in the Theory and Simulations of Polymers (TSP) team at Institute Charles Sadron and the library written by the author in C++ programming language (available at [START_REF]C++ MD trajectory analysis library[END_REF]). The post-analysis, figures and this manuscript itself were done using R programming language for statistical computing and the available packages [START_REF] Core | R: a language and environment for statistical computing[END_REF][START_REF] Wickham | Ggplot2: elegant graphics for data analysis[END_REF][START_REF] Xie | Bookdown: authoring books and technical documents with R markdown[END_REF][START_REF] Xie | Knitr: a comprehensive tool for reproducible research in R[END_REF][START_REF] Xie | R markdown cookbook[END_REF].

The bulk simulations were carried out using LAMMPS build of

Chapter 3

Bulk systems analysis

This chapter presents the results obtained from the simulations of the bulk systems.

We begin with the analysis of the conformational properties in Section 3.1. Then the density and the glass transition are discussed in Section 3.2. Afterwards the static structure of the PBD is analyzed in Section 3.3. The obtained results allow us to check the thermodynamic consistency of the simulations and analyze the heat capacity in Section 3.4. In Section 3.5

we apply the Rouse model to describe the static properties of PBD and give the first insights on the dynamics. A more detailed analysis of the dynamics is made in Section 3.6.

Finally, we discuss our results in Section 3.7.

Conformational properties

Distributions of the torsion angle

A first insight into the influence of reducing the torsional potential on the properties of PBD may be provided by the distribution of the dihedral angle, 𝑝 𝜆 (𝜙), where 𝜆 identifies the different models as defined in Section 2. CRC model is in excellent agreement with Eq. (3.1). This is not unexpected because the dihedral associated with the double bond is a "hard" conformational variable which is kept close to the equilibrium value by its stiff torsional potential. Reducing this potential could change the situation. However, the agreement between the simulated distribution and Eq. (3.1) is still good for both the CRC2 and CRC4 models (Fig. 3.1a). While progressive softening of the torsional potential enhances deviations from the trans configuration, the population at 𝜙 = 0 ∘ or 360 ∘ vanishes, even for the CRC4 model despite the decrease of the potential barrier of the cis state by a factor of 4. This implies that cis/trans isomerization is precluded not only for the CRC model, but also for the CRC2 and CRC4 models. Since the torsional potential for the double bond of the cis conformer is only shifted by 180 ∘ relative to that of trans conformer, the results reported here for the trans conformer are found for the cis conformer, too (not shown).

The difference between the cis and trans conformers disappears when the torsional potential is switched off. Then, Eq. (3.1) predicts a uniform distribution, 𝑝 𝜆=0 (𝜙) = 1/(2𝜋).

The simulated distribution for the FRC model deviates from this prediction (Fig. 3.1a).

The probability of finding 𝜙 near 180 ∘ is increased relative to the uniform distribution, while it is decreased near 𝜙 = 0 ∘ or 360 ∘ . These deviations stem from intrachain LJ interactions between united atoms separated by four and more bonds, which contribute to the conformational energies and thereby create an effective torsional potential. Such an influence of the LJ potential is expected from the force field parametrization developed in Ref. [START_REF] Smith | United atom force field for molecular dynamics simulations of 1,4polybutadiene based on quantum chemistry calculations on model molecules[END_REF], where intramolecular pair interactions were optimized, in conjunction with the other potentials, so as to reproduce correlations between consecutive torsional states (second-order effects between consecutive pairs of torsions and third-order effects involving three consecutive torsions).

For the FRC model this also implies that similar deviations between the simulated and the uniform distributions must be observed for the other dihedrals of the 𝛽, 𝛼 cis and 𝛼 trans bonds. This is illustrated in Figs. UAs adjacent to it, where all dihedral angles were fixed to 𝜙 = 0 ∘ . This geometry results in a distance of only 1.457 Å between the outer UAs. Since these UAs are separated by more than 4 bonds along the chain, the pair interaction between them takes place. According to Eq. (2.5) and Table 2.4, the pair energy at this distance is equal to 𝑈 pair ≈ 10 5 kcal/mol (≈ 10 7 K), which is so high, that the 𝛼 𝑐𝑖𝑠 dihedrals at 𝜙 = 0 ∘ are never adopted during the simulations and are redistributed around the other minima of the potential. This effect is commonly referred to in the literature as the 'steric hindrance' [START_REF] Abe | Conformational energies of n-alkanes and the random configuration of higher homologs including polymethylene[END_REF][START_REF] Ingold | Quantitative study of steric hindrance[END_REF] or 'pentane interference' [START_REF] Smith | Configurational entropy of polyethylene and other linear polymers[END_REF][START_REF] Taylor | Average length and radius of normal paraffln hydrocarbon molecules[END_REF]. These intrachain correlations beyond those of single torsions are clearly important for the conformational properties of PBD [START_REF] Smith | United atom force field for molecular dynamics simulations of 1,4polybutadiene based on quantum chemistry calculations on model molecules[END_REF].

Bond correlations and internal distances

Let ⃗ 𝑟 𝑛 denote the position of the 𝑛th united atom and ⃗ 𝑙 𝑛 = ⃗ 𝑟 𝑛+1 -⃗ 𝑟 𝑛 its bond vector (𝑛 = 1, … , 𝑁). For the united atom 𝑚 = 𝑛 + 𝑠 that is separated by 𝑠 bonds from 𝑛, we define the mean-square internal end-to-end distance 𝑅 2 e (𝑠) and the bond correlation function 𝑃 1 (𝑠),

𝑅 2 e (𝑠) = ⟨( ⃗ 𝑟 𝑛+𝑠 -⃗ 𝑟 𝑛 ) 2 ⟩ , 𝑃 1 (𝑠) = ⟨ ⃗ ℓ 𝑛+𝑠 ⋅ ⃗ ℓ 𝑛 ⟩ 𝑙 2 , (3.2)
where 𝑙 2 = ⟨ ⃗ ℓ 2 𝑛 ⟩ is the mean-square bond length of the PBD chain. Here ⟨…⟩ denotes the thermal average (over all configurations and chains of the system) and also the average over all possible pairs 𝑛, 𝑚 = 𝑛 + 𝑠. (Hence, the statistics will deteriorate for 𝑠 → 𝑁 -1.)

Since ⃗ 𝑟 𝑛+𝑠 -⃗ 𝑟 𝑛 = ∑ 𝑛+𝑠-1 𝑖=𝑛 ⃗ 𝑙 𝑖 , 𝑅 2
e (𝑠) and 𝑃 1 (𝑠) are related to each other by

𝑅 2 e (𝑠) = 𝑙 2 𝑠 + 2𝑙 2 𝑠-1 ∑ 𝑘=1 (𝑠 -𝑘)𝑃 1 (𝑘) . (3.3)
Both quantities are of considerable theoretical importance [2,[START_REF] Semenov | Bond-vector correlation functions in dense polymer systems[END_REF][START_REF] Wittmer | Intramolecular long-range correlations in polymer melts: the segmental size distribution and its moments[END_REF] and have often been studied in polymer melt simulations [START_REF] Svaneborg | Characteristic time and length scales in melts of kremer-grest bead-spring polymers with wormlike bending stiffness[END_REF][START_REF] Wittmer | Intramolecular long-range correlations in polymer melts: the segmental size distribution and its moments[END_REF][START_REF] Wittmer | Long range bond-bond correlations in dense polymer solutions[END_REF][START_REF] Hsu | Static and dynamic properties of large polymer melts in equilibrium[END_REF]. Therefore, we also determine them here.

As an example for the typical behavior of 𝑃 1 (𝑠) found for all models in the interval

𝑃 1 (𝑠) = ⎧ { ⎨ { ⎩ 𝛼 for 𝑠 = 1 , 𝐴 exp(-𝑠/𝑠) for 𝑠 > 1 . (3.4)
where we defined 𝛼 = -⟨cos 𝜃⟩ (> 0). We have analyzed all models in the temperature interval 213 K ≤ 𝑇 ≤ 353 K with Eq. (3.4). The directly measured values for 𝛼 and the fit results for 𝐴 and 𝑠 are compiled in Appendix E.

Insertion of Eq. (3.4) into Eq. (3.3) gives

𝑅 2 e (𝑠) = 𝑙 2 𝑠 [𝐶 ∞ - 2 𝑠 (𝛼 -𝐴𝛼 + 𝐴 𝛼 -𝛼 𝑠+1 (1 -𝛼) 2 )] , (3.5)
where 𝛼 = e -1/𝑠 (< 1) and 𝐶 ∞ is the characteristic ratio of the infinitely long chain,

𝐶 ∞ = lim 𝑁→∞ 𝑅 2 e (𝑁 -1) (𝑁 -1)𝑙 2 = 1 + 2(𝛼 -𝐴𝛼) + 2𝐴𝛼 1 -𝛼 . (3.6)
The characteristic ratio determines the effective bond length 𝑏 e = √𝐶 ∞ 𝑙 [START_REF] Doi | The theory of polymer dynamics[END_REF] and the temperature coefficient 𝜅 of the mean-square end-to-end distance 𝑅 2 e (𝑁 -1) of the chain [START_REF] Flory | Statistical mechanics of chain molecules[END_REF],

𝜅 = 1000 d ln 𝑅 2 e (𝑁 -1) d𝑇 = 1000 d ln 𝐶 ∞ d𝑇 . (3.7)
The factor of 1000 is introduced in this definition because 𝜅 is typically of the order of [START_REF] Smith | United atom force field for molecular dynamics simulations of 1,4polybutadiene based on quantum chemistry calculations on model molecules[END_REF][START_REF] Fetters | Chain dimensions and entanglement spacings[END_REF]. Fig. 3.4a depicts the ratio 𝑅 2 e (𝑠)/𝑠 as a function of 𝑠 for all models at 𝑇 = 353 K. This ratio starts from the mean-square bond length 𝑙 2 (= 2.158 Å 2 ) 1 , then increases first steeply for small 𝑠 and eventually levels off when 𝑠 approaches 𝑁 -1. We find that the dependence 𝑅 2 e (𝑠)/𝑠 on 𝑠 can be well fitted by Eq. (3.5) for all models and temperatures (see Appendix E for a compilation of the fit results). The full line in Fig. 3.4a gives an example for the FRC model. By repeating this analysis at all studied temperatures we obtain the 𝑇 dependence of 𝐶 ∞ (Table 3.1). The results are plotted as ln 𝐶 ∞ versus 𝑇 in Fig. 3.4b. This plot format approximately rectifies the data, yielding negative values for 𝜅 in the range -0.64 K -1 ≲ 𝜅 ≲ -0.52 K -1 for all models studied. A negative value for 𝜅 implies that the chains tend to expand on cooling.

10 -3 K -1
The results for 𝜅 may be compared with experimental values [START_REF] Mays | Temperature dependence of unperturbed dimensions for stereoirregular 1,4-polybutadiene and poly(𝛼-methylstyrene)[END_REF] and Rotational Isomeric State (RIS) calculations [START_REF] Smith | United atom force field for molecular dynamics simulations of 1,4polybutadiene based on quantum chemistry calculations on model molecules[END_REF][START_REF] Abe | Configurational statistics of 1,4-polybutadiene chains[END_REF]. These studies show that the thermal coefficient is sensitive to the stereochemical composition of PBD. While pure cis-PBD has positive 𝜅, negative values are obtained for the pure trans-form. For mixed stereoirregular mi- 1 The stiff bond length potential of the PBD model implies that the bond length is Gaussian distributed around the equilibrium bond length ℓ 0 , leading to a

𝑇 independent ⟨ ⃗ ℓ 2 ⟩ = ℓ 2 0 . Out of 𝑁 -1 bonds, a PBD chain has 𝑁/4 CH=CH bonds of length ℓ 0 = 1.34 Å, 𝑁/4 -1 CH 2 -CH 2 bonds of length ℓ 0 = 1.53 Å, and 𝑁/2 CH-CH 2 bonds of length ℓ 0 = 1.50 Å. The weighted average of ℓ 2 0 gives 𝑙 2 = 2.158 Å 2 .
Eq. (3.5)

𝑙 2
𝑏 2 e (FRC) . (b) Temperature dependence of the characteristic ratio 𝐶 ∞ from Eq. (3.6) for the CRC (red circles), CRC2 (green triangles), CRC4 (purple pluses) and FRC models (blue crosses). The plot format, ln 𝐶 ∞ versus 𝑇, is motivated by Eq. (3.7) when assuming 𝜅 to be constant. The full line presents a linear fit to the FRC data, yielding 𝜅 = -0.52 K -1 . The dashed line shows the fit result to the CRC model, leading to 𝜅 = -0.64 K -1 . The dotted line indicates the experimental result, 𝜅 = -0.10 K -1 (obtained for 298 K ≤ 𝑇 ≲ 373 K), from Table VIII of Ref. [START_REF] Mays | Temperature dependence of unperturbed dimensions for stereoirregular 1,4-polybutadiene and poly(𝛼-methylstyrene)[END_REF].

crostructures similar to our simulation model, RIS calculations [START_REF] Abe | Configurational statistics of 1,4-polybutadiene chains[END_REF] mostly lead to small negative values of 𝜅 (∼ -0.1 K -1 ), in resonance with experimental results (where PBD also contains about 9% of vinyl groups in addition to cis and trans units) [START_REF] Mays | Temperature dependence of unperturbed dimensions for stereoirregular 1,4-polybutadiene and poly(𝛼-methylstyrene)[END_REF]. From this comparison we can conclude that our simulation results are in qualitative accord with available literature data (𝜅 < 0), but the extent of the chain extension on cooling is stronger in the simulation.

Finally, a remarkable observation from the previous analysis is the weak influence of the torsions. Even when switching off all dihedral potentials the conformational properties of the resulting FRC model are very close to those of the CRC model. This observation agrees with the original works [START_REF] Krushev | The role of internal rotational barriers in polymer melt chain dynamics[END_REF][START_REF] Krushev | Intramolecular caging in polybutadiene due to rotational barriers[END_REF]. Qualitatively, we can rationalize this feature by invoking the hindered rotation model that expresses 𝐶 ∞ as [2]:

𝐶 ∞ = 1 + 𝛼 1 -𝛼 1 -⟨cos 𝜙⟩ 1 + ⟨cos 𝜙⟩ . (3.8)
If the torsions were uniformly distributed for the FRC model, we would have ⟨cos 𝜙⟩ = 0 and switching off of the torsional potential would give a big effect on 𝐶 ∞ . However, Fig. 3.1 shows that the distribution of the dihedral angles is not uniform for the FRC model due to the intrachain pair interactions, leading to a finite value of ⟨cos 𝜙⟩. Therefore, we determined ⟨cos 𝜙⟩ for all models. The results are summarized in Table 3.2, together with 𝐶 ∞ from Eq. (3.8) (using the values of 𝛼 from Appendix E). We see that the values of 𝐶 ∞ from Eq. (3.8) for the different models are close to one another at a given 𝑇. This is consistent with the trend shown in Table 3.1. Even at a quantitative level the characteristic ratio from the hindered rotation model is close to the results from Table 3.1. This implies that the dominant contribution to chain stiffness stems on average from local correlations involving three bonds along the chain backbone in PBD.

Density

The comparison of the bulk density of simulated PBD to experimental data has been already presented in Section 2. The colored crosses are the data from 1 𝜇s long NPT runs. The solid black lines represent the fits using Eq. (3.10) with the parameters form Table 3.3. Inset: The fraction of the densities from the cooling (circles) and NVT runs (triangles) divided by the corresponding density from the long NPT runs for each model. rate Γ = 0.83 K/ns, that were, at first, averaged together at each timestep, and then smoothed so that each point is an average of the raw data over a neighboring region of Δ𝑇 = 1 K. The curves overlap at high temperatures before the bending occurs at different temperatures depending on the model. Clearly, such a bending suggests the glass transition that depends on the strength of the dihedrals. We will discuss this later in this section. The overlapping data reveals the mild influence of the dihedrals on the density of PBD in the liquid state.

Honnell et al have studied the equation of state for chain molecules in 1987 [START_REF] Honnell | On the pressure equation for chain molecules[END_REF] and derived that the pressure does not depend on the angular potentials explicitly2 . It is determined by the inter-and intrachain contributions of the site-site correlations, that include only the pair and the bond terms:

( 𝑝𝛽 𝜌 c ) = 𝑁 - 2𝜋𝛽𝜌 c 3 ∑ 𝛼,𝛾 ∫ ∞ 0 𝑔 𝛼,𝛾 inter (𝑟) 𝑑𝑈 pair 𝑑𝑟 𝑟 3 𝑑𝑟 - 4𝜋𝛽 3𝑉 𝑁-1 ∑ 𝑖=1 ∫ ∞ 0 𝑔 𝑖,𝑖+1 intra (𝑟) 𝑑𝑈 bond 𝑑𝑟 𝑟 3 𝑑𝑟 (3.9) - 4𝜋𝛽 3𝑉 𝑁-4 ∑ 𝛼=1 𝑁 ∑ 𝛾=𝛼+4 ∫ ∞ 0 𝑔 𝛼,𝛾 intra (𝑟) 𝑑𝑈 pair 𝑑𝑟 𝑟 3 𝑑𝑟 ,
where 𝜌 c is the number density of chains, 𝑁 is the number of sites in a chain, 𝑔 𝛼,𝛾 inter (𝑟) is the probability to find the sites 𝛼 and 𝛾, that belong to different chains, at a distance 𝑟; 𝑔 𝛼,𝛾 intra (𝑟) is the probability to find the sites 𝛼 and 𝛾, that belong to the same chain, at a distance 𝑟.

Indeed, the pressure components reported by LAMMPS (Fig. 3.6a), that correspond to the contributions arising directly from the bending and torsional potentials are equal to zero for all models. However, the angular potentials may still influence the site-site correlations 𝑔 𝛼,𝛾 inter (𝑟) and 𝑔 𝛼,𝛾 intra (𝑟). In order to estimate these correlations we have calculated the total radial distribution function 𝑔(𝑟) and its inter/intrachain contributions for all models assuming all sites to be equivalent. They are presented in Figs. 3.6b to 3.6d. The curves, that correspond to the intermolecular part, overlap (Fig. 3.6c), so the 1st integral term in Eq. (3.9) is also same for all models. The bonds are very stiff and do not feel the presence of the dihedrals. They result in very high peaks around the equilibrium bond lengths 3.2. DENSITY 51 (they are truncated in Fig. 3.6d), which bring the same contribution to the 2nd integral term in Eq. (3.9) for all models. The intramolecular part has visible differences between the models at 3 Å ≲ 𝑟 ≲ 5 Å (Fig. 3.6d). However, the distances 𝑟 ≲ 4 Å correspond to the correlations only up to the 3rd neighbor along the chain (√𝑅 2 e (𝑠 = 3) ≈ 3.5 Å in Fig. 3.4a). These interactions are excluded in the simulations. The intramolecular part of the pair potential is taken into account starting form the 4th neighbor, which corresponds to the distances 𝑟 ≳ 4 Å (√𝑅 2 e (𝑠 = 4) ≈ 4.5 Å). In this region the data for the 𝑔 intra (𝑟) differs noticeably only close to the peak. We expect these differences to be negligible with The colours denote different models: CRC -red, CRC2 -green, CRC4 -purple, FRCblue.

CHAPTER 3. BULK SYSTEMS ANALYSIS

respect to the other terms of Eq. (3.9). Thus, with sufficient amount of confidence, we can assume the pressure, and therefore, the resulting density of PBD to be the same for all studied models in the melt state.

Consequently, the values of the equilibrium volume obtained from the NPT runs using the CRC model were employed for the NVT simulations of all models at all working temperatures. These data are shown in Fig. 3.5 as black triangles. As a post factum check, we have compared our data to the densities extracted from the 1 𝜇s long NPT runs for each model (colored crosses), that were launched after the NVT runs. These values of the density are obtained with much better averaging statistics and still agree nicely with the rest of the data within the size of the symbol. The inset of Fig. 3.5 shows that the differences between the cooling densities, the densities imposed in the NVT runs and the ones after long NPT runs are less than 1%.

As it was mentioned above, the cooling curves change their slope at different temperatures depending on the model, which suggests an influence of the dihedrals on the glass transition temperature 𝑇 g . There are several ways to extract the 𝑇 g from the evolution of the density upon cooling. A straightforward approach often used in the experiments and simulations is to estimate the 𝑇 g from the intersection point of the linear glassy and liquid regions [START_REF] Schnell | Simulated glass-forming polymer melts: glass transition temperature and elastic constants of the glassy state[END_REF][START_REF] Hutchinson | Determination of the glass transition temperature[END_REF][START_REF] Wunderlich | The nature of the glass transition and its determination by thermal analysis[END_REF]. Even though the choice of the linear regions is ambiguous, this method works reasonably well for the experimental data, because the experimental cooling rates of the order of Γ ∼ 0.001-1 K/s [START_REF] Brüning | Glass transition on long time scales[END_REF][START_REF] Agapov | Does the molecular weight dependence of 𝑇 g correlated to 𝑀 e ?[END_REF][START_REF] Caruthers | Thermo-mechanical signatures of polymeric glasses[END_REF] result in sharply pronounced changes of the slope over only few Kelvins. Still, the ambiguity in the method leads to the problem of reproducibility of the results. In order to solve this problem Dalnoki-Veress et al proposed an alternative method to determine the 𝑇 g [START_REF] Dalnoki-Veress | Molecular weight dependence of reductions in the glass transition temperature of thin, freely standing polymer films[END_REF]. Assuming the thermal expansion coefficients to be constant in the liquid and in the glass (which correspond to the slopes 𝑀 and 𝐺 in the density curves), interpolating across the transition region by a hyperbolic tangent function and integrating the resulting equation, they obtained

𝜌(𝑇 ) = 𝑤 ( 𝑀 -𝐺 2 ) ln [cosh ( 𝑇 -𝑇 g 𝑤 )] + (𝑇 -𝑇 g ) ( 𝑀 + 𝐺 2 ) + 𝜌 g , (3.10)
with 𝜌(𝑇 ) being the density at the temperature 𝑇, 𝑤 the width of the glass transition region, 𝑀 and 𝐺 the slopes d𝜌(𝑇 )/d𝑇 in the melt and glass regions respectively, 𝑇 g the glass transition temperature, 𝜌 g the density at 𝑇 = 𝑇 g . We employed this method in our simulations, which, due to significantly faster cooling rate Γ = 0.83 K/ns are expected there is still an ambiguity in choosing the temperature interval (𝑇 * g -Δ𝑇 , 𝑇 * g + Δ𝑇 ) over which the fit is applied 3 . In order to investigate the influence of the choice of Δ𝑇, we have attempted to fit our smoothed data using Eq. (3.10) over a set of Δ𝑇 ranging from tens of Kelvins to the highest possible value. The results are presented in Fig. 3.7, where the values of the fitted parameters are plotted against Δ𝑇. The lower bound of Δ𝑇 is limited by the amount of smoothing, since it is meaningless to apply the fit over a region of the fluctuation of the data. The upper bound of Δ𝑇 is given by the distance to the closest end of the available data. According to Fig. 3.7, the 𝑇 g , 𝜌 g and the slopes 𝑀 and 𝐺 converge to the corresponding plateaus, but the transition width parameter 𝑤 appears to be the most sensitive to Δ𝑇.

In order to obtain the final values of the parameters, we fixed the 𝑇 g , 𝜌 g , 𝑀 and 𝐺 to their plateau values and refitted the data for the highest Δ𝑇 in order to obtain an 'averaged' 𝑤. These data are presented in Table 3.3 and as dashed lines in Fig. 3.7. In case of the FRC model this resulted in slightly higher value of 𝑤 with respect to the Δ𝑇dependent fit. The fitted curves are presented in Fig. 3.5 as solid black lines.

The most important fitting parameter is, of course, the glass transition temperature 𝑇 g . The value for the CRC model is reasonably higher than the experimental 𝑇 g ≈ 175 K (for high molecular weights) [START_REF] Agapov | Does the molecular weight dependence of 𝑇 g correlated to 𝑀 e ?[END_REF][START_REF] Hintermeyer | Molecular weight dependence of glassy dynamics in linear polymers revisited[END_REF] again due to much faster cooling rate. Weakening the dihedrals significantly shifts the 𝑇 g to lower values. The value for the FRC model qualitatively agrees with the results of the coarse-grained simulations, where typical values of 𝑇 g are of the order of a fraction of 𝜖/𝑘 B (≈ 55 K). [START_REF] Hsu | A coarse-grained polymer model for studying the glass transition[END_REF][START_REF] Schnell | Simulated glass-forming polymer melts: glass transition temperature and elastic constants of the glassy state[END_REF]. 3 The choice of the initial value of 𝑇 * g for the fit interval is arbitrary. We observed no significant influence on the results of the fits after choosing the value of 𝑇 * g visually or iterating after the first fitting attempts. For PBD the CH 2 united atom (UA) is the big particle and CH united atom the small particle. Following commonly employed notations we refer to the big particle as the A particle (= CH 2 ) and to the small particle as the B particle (= CH).

For PBD the LJ interaction parameters are

𝜖 AA = 0.107639 kcal mol , 𝜎 AA = 4.009 Å , 𝜖 BB = 0.114999 kcal mol , 𝜎 BB = 3.385 Å , 𝜖 AB = 0.116723 kcal mol , 𝜎 AB = 3.793 Å , (3.11)
while for the KA mixture they are [START_REF] Kob | Testing mode-coupling theory for a supercooled binary Lennard-Jones mixture: the van Hove correlation function[END_REF] 𝜖 AA = 1.0 , 𝜎 AA = 1 , 𝜖 BB = 0.5 , 𝜎 BB = 0.88 , 𝜖 AB = 1.5 , 𝜎 AB = 0.8 .

(3.12)

These parameters imply that the particle size ratios 𝛿 𝜎 = 𝜎 BB /𝜎 AA ,

𝛿 PBD 𝜎 = 0.844 and 𝛿 KA 𝜎 = 0.88 , (3.13)
are nearly the same for PBD and the KA mixture. Moreover, both systems deviate from the Lorentz rule [START_REF] Hansen | Theory of simple liquids[END_REF], 𝜎 𝑖𝑗 ≠ (𝜎 𝑖𝑖 + 𝜎 𝑗𝑗 )/2, and are thus "nonadditive". Following Eq. ( 8) of

Ref. [START_REF] Ninarello | Models and algorithms for the next generation of glass transition studies[END_REF] we can quantify deviations from additivity by

𝜎 𝑖𝑗 = 1 2 (𝜎 𝑖𝑖 + 𝜎 𝑗𝑗 )[1 -𝑓 𝜎 |𝜎 𝑖𝑖 -𝜎 𝑗𝑗 |] . (3.14) CHAPTER 3. BULK SYSTEMS ANALYSIS
This modified mixing rule implies that identical particles (𝑖 = 𝑗) have unaltered diameters, whereas the cross term of unlike particles deviates from additivity for 𝑓 𝜎 ≠ 0. If 𝑓 𝜎 > 0, the small and large particles overlap more than they could for additive systems, while they overlap less for 𝑓 𝜎 < 0. For PBD and the KA mixture we find

𝑓 PBD 𝜎 = -0.0416 Å -1 ⇒ 𝑓 PBD 𝜎 𝜎 AA ≃ -0.1668 , 𝑓 KA 𝜎 = 1.24 . (3.15)
Compared to the KA mixture, the overlap between unlike particles is decreased for PBD, but the deviation from nonadditivity is fairly weak-weaker by a factor of about 7 relative to the KA mixture.

We can compare the LJ energy parameters of both models in a similar way. First, the LJ energy ratio 𝛿 𝜖 = 𝜖 BB /𝜖 AA have the following values

𝛿 PBD 𝜖 ≃ 1.036 and 𝛿 KA 𝜖 = 0.5 . (3.16)
For PBD the LJ energy between like particles is almost the same. This is not unexpected given that A (= CH 2 ) and B particles (= CH) are chemically very close.

For simple mixtures the interaction energy is often modeled by the Berthelot rule [START_REF] Hansen | Theory of simple liquids[END_REF], [START_REF] Ingebrigtsen | Crystallization instability in glass-forming mixtures[END_REF]. Due to this instability in favor of crystal nucleation by composition fluctuations Ref. [START_REF] Ingebrigtsen | Crystallization instability in glass-forming mixtures[END_REF] concludes that the KA mixture "has a fatal flaw as a glass former".

𝜖 𝑖𝑗 = √ 𝜖 𝑖𝑖 𝜖 𝑗𝑗 . It
These findings for simple liquids invite to explore composition fluctuations also for PBD. This chapter is a first step in this direction by focusing on the partial static structure factors and the related Bhatia-Thornton structure factors. We introduce these structure factors below and apply them to PBD. A spin-off of this analysis is a temperature dependence of the isothermal compressilibity, which we also determine by a volumefluctuation method.

Partial static structure factors: Background for Binary

Mixtures of Simple Liquids

We Let us introduce the collective static structure factor,

𝑆(𝑞) = 1 𝑁 ⟨𝜌( ⃗ 𝑞)𝜌(-⃗ 𝑞)⟩ , (3.19)
which is defined in terms of the coherent density fluctuations for wave vector ⃗ 𝑞,

𝜌( ⃗ 𝑞) = 𝑁 ∑ 𝑗=1 exp (i ⃗ 𝑞 ⋅ ⃗ 𝑟 𝑗 ) (for 𝑞 ≠ 0), (3.20)
where ⟨…⟩ denotes the canonical average and ⃗ 𝑟 𝑗 is the position of particle 𝑗 (= 1, … , 𝑁).

For spatially homogeneous and isotropic systems, the structure factor depends only on the modulus of the wave vector, 𝑞 = | ⃗ 𝑞|. In the limit 𝑞 → 0, 𝑆(𝑞) is small because the fluctuations of the particle number relative to the average ⟨𝑁 ⟩ (= 𝑁) are weak in a dense system (at constant chemical potential),

𝑆(𝑞 → 0) = ⟨𝑁 2 ⟩ -⟨𝑁 ⟩ 2 ⟨𝑁 ⟩ ≪ 1. (3.21)
For a system consisting only of one particle type we also have [START_REF] Hansen | Theory of simple liquids[END_REF]

𝑆(𝑞 → 0) = 𝑘 B 𝑇 𝜌𝜅 𝑇
where 𝑇 is the temperature, 𝜌 is the particle number density (𝜌 = 𝑁 /𝑉 with 𝑉 the volume of the system) and 𝜅 𝑇 is the isothermal compressibility.

Further insight into structural properties of a binary mixture can be obtained from the partial static structure factors

𝑆 𝛼𝛽 (𝑞) = 1 𝑁 ⟨𝜌 𝛼 ( ⃗ 𝑞)𝜌 𝛽 (-⃗ 𝑞)⟩ (𝛼, 𝛽 = A, B) (3.22)
defined by the partial density fluctuations

𝜌 𝛼 ( ⃗ 𝑞) = 𝑁 𝛼 ∑ 𝑗 𝛼 =1 exp (i ⃗ 𝑞 ⋅ ⃗ 𝑟 𝑗 𝛼 ) , (3.23)
where ⃗ 𝑟 𝑗 𝛼 is the position of particle 𝑗 𝛼 of species 𝛼. As 𝜌 A ( ⃗ 𝑞) + 𝜌 B ( ⃗ 𝑞) = 𝜌( ⃗ 𝑞) ≡ 𝜌 n ( ⃗ 𝑞), the collective structure factor can be expressed as

𝑆(𝑞) = 𝑆 nn (𝑞) = 𝑆 AA (𝑞) + 𝑆 BB (𝑞) + 2𝑆 AB (𝑞). (3.24)
While 𝑆 𝛼𝛽 (𝑞) characterize spatial correlations between like or unlike particles, 𝑆(𝑞) describes number-number (nn) correlations (hence the notation 𝑆 = 𝑆 nn ). Eq. (3.24) is not the only physically significant linear combination of the partial structure factors. Since composition (or concentration) fluctuations 𝜌 c ( ⃗ 𝑞) are defined by

𝜌 c ( ⃗ 𝑞) = 𝜌 A ( ⃗ 𝑞) -𝑥 A 𝜌( ⃗ 𝑞) = 𝑥 B 𝜌 A ( ⃗ 𝑞) -𝑥 A 𝜌 B ( ⃗ 𝑞), the structure factor 𝑆 cc (𝑞) = 1 𝑁 ⟨𝜌 c ( ⃗ 𝑞)𝜌 c (-⃗ 𝑞)⟩ = 𝑥 2 B 𝑆 AA (𝑞) + 𝑥 2 A 𝑆 BB (𝑞) -2𝑥 A 𝑥 B 𝑆 AB (𝑞) (3.25) 3.3. STRUCTURE FACTORS, COMPRESSIBILITY 59
represents composition-composition (cc) correlations, and the structure factor between 𝜌 n and 𝜌 c ,

𝑆 nc (𝑞) = 1 𝑁 ⟨𝜌 n ( ⃗ 𝑞)𝜌 c (-⃗ 𝑞)⟩ = 𝑥 B 𝑆 AA (𝑞) -𝑥 A 𝑆 BB (𝑞) + (𝑥 B -𝑥 A )𝑆 AB (𝑞), (3.26)
describes number-composition (nc) correlations. The structure factors 𝑆 nn (𝑞), 𝑆 cc (𝑞) and 𝑆 nc (𝑞) are often referred to as Bhatia-Thornton structure factors [START_REF] Bhatia | Structural aspects of the electrical resistivity of binary alloys[END_REF]. They have been studied extensively for metallic alloys [START_REF] Horbach | Self-diffusion and interdiffusion in al 80 ni 20 melts: simulation and experiment[END_REF][START_REF] Das | Structural relaxation in a binary metallic melt: molecular dynamics computer simulation of undercooled al 80 ni 20[END_REF] or colloidal suspensions [START_REF] Thorneywork | Structure factors in a two-dimensional binary colloidal hard sphere system[END_REF].

For binary mixtures a deviation between 𝑆(𝑞 → 0) and 𝑘 B 𝑇 𝜌𝜅 𝑇 is expected from the work of Bhatia and Thornton [START_REF] Bhatia | Structural aspects of the electrical resistivity of binary alloys[END_REF] and also from the Kirkwood-Buff theory for multicomponent solutions [START_REF] Ben-Naim | Molecular theory of solutions[END_REF]. For 𝑞 → 0 the Bhatia-Thornton structure factors are related to the thermodynamic properties of the binary mixture:

𝑆(𝑞 → 0) = 𝑘 B 𝑇 𝜌𝜅 𝑇 + 𝛿 2 𝑆 cc (𝑞 → 0), (3.27)
𝑆 cc (𝑞 → 0) = 𝑁 𝑘 B 𝑇 (𝜕 2 𝐺/𝜕𝑥 2 A ) 𝑝,𝑇 ,𝑁 , (3.28)
𝑆 nc (𝑞 → 0) = -𝛿 𝑆 cc (𝑞 → 0), (3.29)
where 𝐺 is the Gibbs free energy, 𝑝 the pressure and

𝛿 = 𝜌(𝑣 A -𝑣 B ) (3.30)
is a volumetric factor given by the partial molar volumes

𝑣 A = 𝜕𝑉 𝜕𝑁 A | 𝑝,𝑇 ,𝑁 B and 𝑣 B = 𝜕𝑉 𝜕𝑁 B | 𝑝,𝑇 ,𝑁 A . (3.31)
Eq. (3.27) shows that in a mixture fluctuations of the total particle number, 𝑆(𝑞 → 0), do not only stem from compressibility effects-that is, from the volume response of the system to a pressure variation-but also from composition fluctuations and their coupling to the number density. Since thermodynamic stability requires (𝜕 2 𝐺/𝜕𝑥 2 A ) 𝑝,𝑇 ,𝑁 > 0, we have in general 𝑆 cc (𝑞 → 0) > 0. This implies that the second term in the right-hand-side of Eq. (3.27) is positive, if 𝛿 ≠ 0 or 𝑣 A ≠ 𝑣 B . The molar volumes of the two species can CHAPTER 3. BULK SYSTEMS ANALYSIS be calculated via the Kirkwood-Buff theory from the partial structure factors in the limit 𝑞 → 0 [START_REF] Ben-Naim | Molecular theory of solutions[END_REF],

𝜌𝑣 A = 𝑥 A 𝑆 BB (0) -𝑥 B 𝑆 AB (0) 𝑥 2 A 𝑆 BB (0) + 𝑥 2 B 𝑆 AA (0) -2𝑥 A 𝑥 B 𝑆 AB (0) , (3.32
)

𝜌𝑣 B = 𝑥 B 𝑆 AA (0) -𝑥 A 𝑆 AB (0) 𝑥 2 A 𝑆 BB (0) + 𝑥 2 B 𝑆 AA (0) -2𝑥 A 𝑥 B 𝑆 AB (0) , (3.33)
and similarly the theory also allows to express the compressibility as

𝑘 B 𝑇 𝜌𝜅 𝑇 = 𝑆 AA (0)𝑆 BB (0) -[𝑆 AB (0)] 2 𝑥 2 A 𝑆 BB (0) + 𝑥 2 B 𝑆 AA (0) -2𝑥 A 𝑥 B 𝑆 AB (0) , (3.34)
where 𝑆 𝛼𝛽 (0) is an abbreviation for 𝑆 𝛼𝛽 (𝑞 → 0). S AA (q) + S BB (q)

Polybutadiene

2S AB (q)

CRC: T=353 K S AA (q) + S BB (q)

2S AB (q)

CRC: lines FRC: symbols

𝑆 AA (𝑞 → ∞) + 𝑆 BB (𝑞 → ∞) = 𝑥 A + 𝑥 B = 1. For large 𝑞, say 𝑞 ≳ 8 Å -1
or 𝑞𝜎 AA ≳ 32, the structure factor 𝑆(𝑞) is thus dominated by correlations between identical particles. The sum 𝑆 AA (𝑞) + 𝑆 BB (𝑞) is positive for all 𝑞, whereas 𝑆 AB (𝑞) oscillates around 0 and remains positive for 𝑞 < 𝑞 * with 𝑞 * being the position of the first peak of 𝑆(𝑞). This behavior for small 𝑞 is different from the one typically found for simple binary mixtures [START_REF] Das | Structural relaxation in a binary metallic melt: molecular dynamics computer simulation of undercooled al 80 ni 20[END_REF][START_REF] Nauroth | Quantitative test of the mode-coupling theory of the ideal glass transition for a binary lennard-jones system[END_REF][START_REF] Moreno | Anomalous dynamic arrest in a mixture of large and small particles[END_REF][START_REF] Ruscher | Glassy dynamics of a binary voronoi fluid: a mode-coupling analysis[END_REF]:

For simple liquids 𝑆 AB (𝑞) is negative for 𝑞 ≈ 𝑞 * and goes through a minimum when 𝑞 decreases. This minimum outweighs the positive contribution of 𝑆 AA (𝑞) + 𝑆 BB (𝑞), leading to a dip in 𝑆(𝑞) before 𝑆(𝑞) increases toward a plateau as 𝑞 → 0. Such a dip is not observed for the PBD melt. Here 𝑆(𝑞) continuously decreases-toward the compressibility plateau 𝑆(𝑞 → 0) = 𝑘 B 𝑇 𝜌𝜅 𝑇 , as Fig. 3.9 indicates and will be discussed further.

Using the partial static structure factors and

𝑥 A = 𝑥 B = 1/2 (PBD
𝜌𝑣(𝑞) = 𝜌𝑣 0 + 𝐴 𝑣 𝑞 2 + 𝐵 𝑣 𝑞 4 , (3.35)
where 𝜌𝑣 0 , 𝐴 𝑣 and 𝐵 𝑣 are fit parameters. The values of 𝛿 are obtained using Eq. (3.30).

They are given in Table 3.4.

Similarly, we parameterize the 𝑞 dependence of 𝑆 cc (𝑞) and the 𝑞 dependent compressibility for small 𝑞:

𝑆 fit cc (𝑞) = 𝐴 cc 𝑞 2 + 𝐵 cc 𝑞 4 , (3.36)
k B T ρκ T (q) = 0.188 q (Å -1 ) An alternative-and often employed-method to determine the compressibility is to monitor the volume fluctuations in NPT simulations (as a special case of the general theory of the ensemble dependence of fluctuations [START_REF] Lebowitz | Ensemble dependence of fluctuations with application to machine computations[END_REF][START_REF] Wittmer | Communication: pressure fluctuations in isotropic solids and fluids[END_REF])

1 x A x B
T = 353 K δ CRC = ρ(v A -v B ) = 1.25 δ FRC = ρ(v A -v B ) =
ρv A (q), ρv B (q)
𝜅 𝑇 = - 1 𝑉 𝜕𝑉 𝜕𝑝 | 𝑇 = ⟨𝑉 2 ⟩ -⟨𝑉 ⟩ 2 𝑘 B 𝑇 ⟨𝑉 ⟩ . (3.38)
We employed this method to determine 𝑘 B 𝑇 𝜌𝜅 𝑇 from 1 𝜇s long NPT runs using all our models of PBD. q (Å -1 ) S(q), k B T ρκ are included in Fig. 3.12. The crosses depict MD data for cis-1,4 PBD from Ref. [START_REF] Tsolou | Temperature and pressure effects on local structure and chain packing in cis-1,4-polybutadiene from detailed molecular dynamics simulations[END_REF] and the squares experimental data taken from Ref. [START_REF] Schweizer | Theory of dynamic barriers, activated hopping, and the glass transition in polymer melts[END_REF]. We find reasonable agreement between these literature and the present results with the volume fluctuations data being the most accurate. The data obtained from the partial structure factors appears to deviate significantly, even though we expect them to be more precise than a straightforward extrapolation of 𝑆(𝑞 → 0). Of course, our data for the structure factors at low 𝑞 are still quite noisy even with the averaging statistics corresponding to several 𝜇s of the simulations.

T (q) colour CRC FRC shape S(q) k B T ρκ T (q) (a)
However, we believe there can still be a room for improvement in the choice of the fitting

1 √𝑘 B 𝑇 𝜌𝜅 𝑇 = -𝐴 + 𝐵 𝑇 . (3.39)
Eq. (3.39) was not claimed to be "quantitatively reliable nor rigorous, but merely suggestive of how experimental data for the dimensionless compressibility might be organized" (cf. Appendix A in [START_REF] Schweizer | Theory of dynamic barriers, activated hopping, and the glass transition in polymer melts[END_REF]). However, Eq. (3.39) was found to be practically relevant, since it linearizes the experimental data for numerous polymers, including PBD (cf. Figs. 4 and 5 in [START_REF] Schweizer | Theory of dynamic barriers, activated hopping, and the glass transition in polymer melts[END_REF]). We reproduce the (scanned) PBD data and the fit to Eq. (3.39) from [START_REF] Schweizer | Theory of dynamic barriers, activated hopping, and the glass transition in polymer melts[END_REF] in Fig. 3.12.

More recently, Mirigian and Schweizer derived Eq. (3.39) from a low-temperature expansion of a van-der-Waals (vdW) model for the equation of state [START_REF] Mirigian | Elastically cooperative activated barrier hopping theory of relaxation in viscous fluids: II. Thermal liquids[END_REF]. Eq. (3.39) enters the mapping of a hard-sphere fluid to molecular liquids [START_REF] Mirigian | Elastically cooperative activated barrier hopping theory of relaxation in viscous fluids: II. Thermal liquids[END_REF] and polymer melts [START_REF] Mirigian | Dynamical theory of segmental relaxation and emergent elasticity in supercooled polymer melts[END_REF]. This mapping underlies the elastically cooperative nonlinear Langevin equation (ECNLE) theory for the 𝑇 dependence of the 𝛼 relaxation time [START_REF] Mirigian | Elastically cooperative activated barrier hopping theory of relaxation in viscous fluids: II. Thermal liquids[END_REF][START_REF] Mirigian | Dynamical theory of segmental relaxation and emergent elasticity in supercooled polymer melts[END_REF] and extensions of the ECNLE theory to polymer films [START_REF] Phan | Dynamic gradients, mobile layers, 𝑇 g shifts, role of vitrification criterion, and inhomogeneous decoupling in free-standing polymer films[END_REF][START_REF] Mirigian | Influence of chemistry, interfacial width, and non-isothermal conditions on spatially heterogeneous activated relaxation and elasticity in glass-forming free standing films[END_REF][START_REF] Schweizer | Progress towards a phenomenological picture and theoretical understanding of glassy dynamics and vitrification near interfaces under nanoconfinement[END_REF]. The results for 𝑘 B 𝑇 𝜌𝜅 𝑇 , found for all models are indicated by triangles. The orange crosses present the MD results from Ref. [START_REF] Tsolou | Temperature and pressure effects on local structure and chain packing in cis-1,4-polybutadiene from detailed molecular dynamics simulations[END_REF] for cis-1,4-PBD at 1 atm (systems with 128 C-atoms per chain and 32 chains), which were calculated by scanning the data from Fig. 4 and Fig. 7 of [START_REF] Tsolou | Temperature and pressure effects on local structure and chain packing in cis-1,4-polybutadiene from detailed molecular dynamics simulations[END_REF].

Interestingly, essentially the same vdW approach as employed by Schweizer and coworkers was proposed before by Long and Lequeux [START_REF] Long | Heterogeneous dynamics at the glass transition in van der waals liquids, in the bulk and in thin films[END_REF] to model the pressure-volumetemperature (PVT) behavior of polymer melts. Starting from a vdW-like ansatz for the free energy and imposing equilibrium conditions Long and Lequeux get the following expressions for the density

𝜌 = 𝜌 0 2 1 + √1 -𝑇 /𝑇 c 1 + 𝛽(1 -√1 -𝑇 /𝑇 c )/2 (3.40)
and the (inverse) compressibility As shown in Fig. 3.12, the fits give good agreement, yielding values for the parameters (𝜌 0 , 𝑇 c , 𝛽, 𝛾) that are in reasonable accord with literature results for other hydrocarbon polymers, such as polystyrene or poly(isobutylene) [START_REF] Long | Heterogeneous dynamics at the glass transition in van der waals liquids, in the bulk and in thin films[END_REF].

1 𝜅 𝑇 = 𝑘 B 𝑇 𝜌 0 2𝛾 2 [1 + √1 - 𝑇 𝑇 c ] 3 × [ 1 (1 -√1 -𝑇 /𝑇 c ) 2 - 𝑇 c 𝑇 ] . ( 3

Energies and heat capacity

Let 𝑥 denote the microstate of the PBD melt containing 𝑁 c polymers with 𝑁 united atoms each. The hamiltonian is written as

ℋ(𝑥) = 𝐾(𝑥) + 𝑈 tot (𝑥) , (3.42)
where 𝐾 is the total kinetic energy and 𝑈 tot is the total potential energy consisting of valence terms (bond, bond angle, dihedral) and pair (Lennard-Jones) interactions.

As explained in Section 2. For 𝑝 = 0 the enthalpy is given by 𝑈 (𝑇 , 𝑝).

In the following sections we discuss the contributions to 𝑈 (𝑇 , 𝑝) from the potential energy of PBD before we turn to heat capacity.

Hard variables: Bond length

The bond length ℓ is subjected to the harmonic bond potential defined by

𝑈 bond (ℓ) = 𝐾 bond (ℓ -ℓ 0 ) 2 , (3.45)
where 𝐾 bond is (half) the force constant and ℓ 0 the equilibrium bond length. The parameters for these constants are given in Table 2.1.

In Section 2.4.2 we observed an agreement of the bond length distribution obtained from the simulations with the one for independent bonds (i.e. with the distribution determined only by Eq. (3.45)). This allows us to assume that the bonds are independent of each other and also independent of the potential energies for bond angles, torsion angles and pair interactions. Then, the partition function for a bond reads

𝑍 bond (𝛽) = 4𝜋 ∫ ∞ 0 dℓ ℓ 2 exp [-𝛽𝐾 bond (ℓ -ℓ 0 ) 2 ] ≃ 4𝜋ℓ 2 0 √ 𝜋 𝛽𝐾 bond , (3.46)
where 𝛽 = 1/(𝑘 B 𝑇 ).

The last approximate equality was obtained by utilizing Eq. (3.462/7) from Ref. [START_REF] Gradshteyn | Table of integrals, series, and products[END_REF] and the fact that the dimensionless quantity (𝑘 B = 1.9872067 × 10 -3 kcal/(mol,K))

𝛽𝐾 bond ℓ 2 0 = 2500 1.9872067 × 10 -3
⟨𝑈 bond ⟩ = - 𝜕 ln 𝑍 bond (𝛽) 𝜕𝛽 ≃ 1 2 𝑘 B 𝑇 , (3.48)
as expected for a one-dimensional harmonic oscillator according to the equipartition theorem.

PBD has three types of bonds:

• 1 st bond (b1) = CH-CH double bond: There is one such bond per monomer. So the total number of b1 bonds is given by

𝑁 b1 = # of monomers = 𝑛 = 𝑁 4 ,
where 𝑛 (= 29) is the number of monomers per chain and 𝑁 (= 116) is the number of united atoms per chain.

• 2 nd bond (b2) = CH-CH 2 bond: There are two such bonds per monomer. So the total number of b2 bonds is given by

𝑁 b2 = 2𝑛 = 𝑁 2 .
• 3 rd bond (b3) = CH 2 -CH 2 bond: This bond links consecutive monomers. So the total number of b3 bonds is given by

𝑁 b3 = 𝑛 -1 = 𝑁 4 -1 .
When adding up the energies for all bonds in the melt and using Eq. (3.48), we find for the average extensive bond energy (𝐸 bond ) of the system

𝐸 bond (𝑇 , 𝑝, 𝑁 c 𝑁 ) ≃ 𝑁 c 𝑁 (1 - 1 𝑁 ) 1 2 𝑘 B 𝑇 , (3.49)
where 𝑁 c (= 40) is the number of chains in the system. However, this deviation from equilibrium is much weaker than for the other potential energies, as we will see below.

Hard variables: Bond angle

An approach analogous to that of the bond length can be used to estimate the average energy of the bond angle. The bond angle 𝜃 is subjected to the potential

𝑈 bend (𝜃) = 𝐾 bend (cos 𝜃 -cos 𝜃 0 ) 2 , (3.50)
where the values for 𝐾 bend and 𝜃 0 are given in justifies the assumption of the angles being independent of each other and of the other variables. The partition function of an angle is then given by

𝑍 bend (𝛽) = ∫ 𝜋 0 d𝜃 sin 𝜃 exp [-𝛽𝐾 bend (cos 𝜃 -cos 𝜃 0 ) 2 ] = 1 2 √ 𝜋 𝑞 [ erf ([1 -𝑥 0 ]𝑞) + erf ([1 + 𝑥 0 ]𝑞)] , (3.51)
where erf(𝑥) is the error function, 𝑞 = √𝛽𝐾 bend and 𝑥 0 = cos 𝜃 0 . With the values from For the average energy of a bond angle this implies

⟨𝑈 bend ⟩ = - 𝜕 ln 𝑍 bend (𝛽) 𝜕𝛽 ≃ 1 2 𝑘 B 𝑇 . (3.53)
The PBD model has two types of bond angles:

• 1 st angle (a1) = CH 2 -CH-CH angle: There are two angles per monomer, leading to

𝑁 a1 = 2𝑛 = 𝑁 2
for the total number of a1 angles.

• 2 nd angle (a2) = CH-CH 2 -CH 2 angle: A single monomer has no a2 angle, two monomers have 2 a2 angles, three monomers have 4 a2 angles, etc. So the total number of a2 angles is given by

𝑁 a2 = 2(𝑛 -1) = 𝑁 2 (1 - 4 𝑁 ) .
Adding up the energies for all angles in the system and using Eq. (3.53), we get for the average extensive bond-angle energy

𝐸 bend (𝑇 , 𝑝, 𝑁 c 𝑁 ) ≃ 𝑁 c 𝑁 (1 - 2 𝑁 ) 1 2 𝑘 B 𝑇 . (3.54)
This is the same result as for the bond length, bearing the (small 1/𝑁) finite-chain length correction.

Fig. 3.14 shows the results for the bond angle; they are similar to those for bond length.

The intensive energies 𝑒 bend (= 𝐸 bend /(𝑁 c 𝑁 )) from the cooling runs and equilibrium simulations agree with each other for 𝑇 ≥ 213 K and also with Eq. (3.54). However, for 𝑇 → 0 systematic deviations from Eq. (3.54) are clearly visible and much larger than for the bond length. The onset of these deviations occur in the glass transition zone around 𝑇 g , as illustrated for the CRC model.

Soft variables: Dihedral angles

As a start to the analysis of the dihedral angles, we can make the same assumption as for the bond lengths and bond angles: we take the dihedral angles to be independently distributed according to the Boltzmann weight of their potential energies. Section 3.1.1 compares the predicted distributions to those from the MD simulations and concludes the following :

• Double bond cis and trans dihedrals

The simulated distributions agree with the predictions from the Boltzmann weights for the CRC model and the reduced barrier models (Fig. 3.1a). This is expected, since the torsions of the double bonds are in fact hard variables. Their potentials are very stiff, allowing only for weak vibrations around the planar configuration at all temperatures studied. The vibrations contribute 𝑘 B 𝑇 /2 per double bond to the average dihedral energy.

• 𝛼 trans and 𝛽 dihedrals

The potential energies of these dihedrals have pronounced minima but the barriers between them are more than an order magnitude smaller than well depth of the double bond dihedrals. Due to this softening of the potentials some deviations from the predicted distributions must be expected. Indeed, Fig. 3.1b show that the predictions from the Boltzmann weights are not as good as for the double bonds.

Still, the agreement is semiquantitative.

• 𝛼 cis dihedrals Fig. 3.1d show strong deviations between predicted and simulated distributions. In particular, the simulated distributions tend to 0 for 𝜙 → 0 ∘ and 𝜙 → 360 ∘ , whereas the predicted distributions have a maximum for these limiting angles. Clearly, the assumption of independent dihedrals fails here. Rotation around the 𝛼 cis bond leads to strong steric hindrance due to nonbonded pair interactions for 𝜙 → 0 ∘ and 𝜙 → 360 ∘ , which suppresses the occurrence of these limiting angles. This influence of the nonbonded interactions is similar to the "pentane effect" for polyethylene [START_REF] Flory | Statistical mechanics of chain molecules[END_REF].

A harmonic expansion around the minima of the dihedral potentials cannot be sufficient to describe the 𝑇 dependence of the average torsion energy 𝐸 tor (𝑇 , 𝑝, 𝑁 c 𝑁 ) in the polymer liquid. Anharmonic effects need to be taken into account. A possible scheme • 𝑛 -1 beta dihedrals with 3 minima,

• 2𝑛𝑓 trans alpha-trans dihedrals with 3 minima.

Then, the harmonic contribution to 𝐸 tor (𝑇 , 𝑝, 𝑁 c 𝑁 ) reads

𝑁 c 𝑁 4 [1 + 3(1 -4/𝑁 ) + 6𝑓 trans ] 1 2 𝑘 B 𝑇 = 𝑁 c 𝑁 𝐵 approx 𝜙 𝑇 ⇒ 𝐵 approx 𝜙 = 1.788 × 10 -3 kcal/mol K ,
which is close to the fit value for 𝐵 𝜙 from Fig. 3.15.

The impact of the glass transition is clearly visible from the main panel of Fig. 3.15.

Below 213 K the cooling process is too fast for the melt to stay at equilibrium. Accordingly, the dihedral energy at low 𝑇 is larger than the extrapolated liquid curve. For 𝑇 → 0 a harmonic approximation provides an excellent description of 𝐸 tor (𝑇 ) (cf. solid line in Fig. 3.15). However, the fitted ground state energy, 𝐴 g 𝜙 = -1.104315 kcal/mol, is much larger than the one expected from torsional potential, which is given by the energy minimum at 180 ∘ of the 𝛽 dihedral (= -5.28281 kcal/mol). This exemplifies the nonequilibrium character of the glass created by the cooling process.

Soft variables: Pair interactions

From the viewpoint of the pair interactions the PBD model corresponds to a binary Lennard-Jones (LJ) mixture (Section 3.3.1). It is therefore interesting to survey in the literature how the excess energy of LJ and other simple liquids has been described.

In this respect, an influential theory was developed by Rosenfeld and Tarazona (RT) [START_REF] Rosenfeld | Density functional theory and the asymptotic high density expansion of the free energy of classical solids and fluids[END_REF]. The RT theory starts from an accurate approximation for the excess free energy functional of densely packed hard spheres (so-called fundamental measure functionals) and uses thermodynamic perturbation theory to extend the theory to continuous ("soft") potentials. First-order perturbation theory around 𝜂 = 1 packing fraction gives for the excess internal energy (𝐸 ex ) of a three-dimensional classical bulk fluid with 𝑁 particles 3.4. ENERGIES AND HEAT CAPACITY 77 the following constant volume (particle density 𝜌) and constant temperature expression (cf. Eq. ( 34) in Ref. [START_REF] Rosenfeld | Density functional theory and the asymptotic high density expansion of the free energy of classical solids and fluids[END_REF]):

𝐸 ex (𝜌, 𝑇 ) 𝑁 = 𝐴 RT (𝜌) + 𝐵 RT (𝜌)𝑇 3/5 + ⋯ .

The leading 𝑇 3/5 term is expected to be accurate for large densities (near freezing) and predominanly repulsive interactions (see the detailed discussion in [START_REF] Rosenfeld | Density functional theory and the asymptotic high density expansion of the free energy of classical solids and fluids[END_REF]). It implies that the excess isochoric heat capacity (𝐶 ex 𝑉 = 𝜕𝐸 ex /𝜕𝑇 | 𝑉 ,𝑁 ) increases with decreasing temperature as 𝐶 ex 𝑉 ∝ 𝑇 -2/5 . A detailed simulation study of the 𝑇 dependence of 𝐶 ex 𝑉 was carried out in Ref. [START_REF] Ingebrigtsen | Communication: the rosenfeldtarazona expression for liquids' specific heat: a numerical investigation of eighteen systems[END_REF] for 18 model liquids with different stoichiometric compositions, molecular topologies and interactions. It was found that the RT expression 𝐶 ex 𝑉 ∝ 𝑇 -2/5 provides a better approximation for liquids with strong correlations between equilibrium fluctuations of virial and potential energy, i.e., for so-called "Roskilde-simple" liquids [START_REF] Ingebrigtsen | What is a simple liquid?[END_REF]. The Kob-Andersen binary LJ mixture [START_REF] Ingebrigtsen | Crystallization instability in glass-forming mixtures[END_REF][START_REF] Kob | Testing mode-coupling theory for a supercooled binary Lennard-Jones mixture: the van Hove correlation function[END_REF][START_REF] Kob | Testing mode-coupling theory for a supercooled binary Lennard-Jones mixture: intermediate scattering function and dynamic susceptibility[END_REF][START_REF] Andersen | Molecular dynamics studies of heterogeneous dynamics and dynamic crossover in supercooled atomic liquids[END_REF] is a Roskilde-simple liquid (for 𝜌 = 1.2 in LJ units) [START_REF] Ingebrigtsen | Communication: the rosenfeldtarazona expression for liquids' specific heat: a numerical investigation of eighteen systems[END_REF].

It is therefore tempting to test whether the following ansatz,

𝐸 pair (𝑇 , 𝑝, 𝑁 c 𝑁 ) ≃ 𝑁 c 𝑁[𝐴 RT (𝑝) + 𝐵 RT (𝑝)𝑇 3/5 ] , (3.56)
can provide a good description for the pair interactions of the PBD model. The impact of the glass transition is clearly visible in the main panel of Fig. 3.16. For the CRC model we find, as for the dihedral angles, that the cooling process is too fast for melt to stay at equilibrium below 213 K. Accordingly, the pair energy at low 𝑇 is larger than the extrapolated liquid curve. For 𝑇 → 0 a harmonic approximation provides again an excellent description of the energy (cf. dash-dotted line in Fig. 3.16), as expected for a solid.

Heat capacity of the polymer liquid: results for the CRC model

The thermodynamic relation between the (extensive) heat capacities at constant pressure 𝐶 𝑝 and at constant volume 𝐶 𝑉 is given by

𝐶 𝑝 = 𝐶 𝑉 + 𝑇 𝑉 𝛼 2 𝜅 𝑇 = 𝐶 𝑉 + 𝑁 c 𝑁 𝑇 𝛼 2 𝜌𝜅 𝑇 , (3.57)
where 𝛼 is the thermal expansion coefficient,

𝛼 = 1 𝑉 𝜕𝑉 𝜕𝑇 | 𝑝,𝑁 c 𝑁 ,

ENERGIES AND HEAT CAPACITY

79

𝜅 𝑇 the isothermal compressibility,

𝜅 𝑇 = - 1 𝑉 𝜕𝑉 𝜕𝑝 | 𝑇 ,𝑁 c 𝑁
, and 𝜌 = 𝑁 c 𝑁 /𝑉 is the number density of united atoms in the system (we recall that 𝑁 = 4𝑛 where 𝑛 is the number of monomers per chain).

If we consider the heat capacity as the ability of a material to take up energy and to statistically distribute this energy over all degrees of freedom, it is natural to report the specific heat (𝑐 = 𝐶/(𝑁 c 𝑁 )) divided by 𝑘 B because this quantity can be interpreted as the number of degrees of freedom per particle (united atom) contributing to the heat capacity at temperature 𝑇.

Fig. 3.17 shows such a plot for the CRC model. The specific heat at constant pressure 𝑐 𝑝 was obtained from the fluctuations of the enthalpy 𝐻 = 𝑈 + 𝑝𝑉 in the NPT runs (circles) 4 . Alternatively, 𝑐 𝑝 can also be calculated by

𝑐 𝑝 = 1 𝑁 c 𝑁 𝜕𝐻 𝜕𝑇 | 𝑝,𝑁 c 𝑁 = 1 𝑁 c 𝑁 𝜕𝑈 𝜕𝑇 | 𝑝,𝑁 c 𝑁 , (3.58)
where 𝑈 (𝑇 , 𝑝, 𝑁 c 𝑁 ) is the internal energy. The 𝑝𝑉 term of 𝐻 does not contribute because the NPT simulations are performed at 𝑝 = 0. The internal energy is given by 𝑈 = 𝐸 kin + 𝐸 bond + 𝐸 bend + 𝐸 tor + 𝐸 pair (3.59)

with 𝐸 kin being the (average extensive) kinetic energy,

𝐸 kin = ⟨𝐾⟩ = 3 2 𝑁 c 𝑁 𝑘 B 𝑇 .
Adding 𝐸 kin and the energies for the bond length [Eq. [START_REF] Grebowicz | Heat capacities of 1,4-polybutadienes[END_REF]. The green squares depict the results from Eq. (3.67) using the experimental data (divided by 4𝑅) for the solid state from Fig. 9 of [START_REF] Grebowicz | Heat capacities of 1,4-polybutadienes[END_REF] as a proxy to account for contributions from quantum mechanical vibrations in the experimental system.

For the "soft" dihedral and pair interactions we find from Eq. (3.55) and Eq. (3.56) the pair interactions, as seen from dash-dotted line in Fig. 3.17.

1 𝑁 c 𝑁 𝜕𝐸 dih 𝜕𝑇 | 𝑝,𝑁 c 𝑁 = 𝐵 𝜙 -2|𝐶 𝜙 |𝑇 , (3.61
The specific heat at constant volume can be calculated from the fluctuations of the hamiltonian ℋ in the NVT simulations

𝑐 𝑉 = 1 𝑁 c 𝑁 1 𝑘 B 𝑇 2 [⟨ℋ 2 ⟩ -⟨ℋ⟩ 2 ] .
(3.63)

The squares in Fig. 3.17 show the result of this calculation. As expected, 𝑐 𝑉 < 𝑐 𝑝 .

Alternatively, 𝑐 𝑉 can be obtained from 𝑐 𝑝 via Eq. (3.57)

𝑐 𝑉 = 1 𝑁 c 𝑁 𝜕𝐻 𝜕𝑇 | 𝑝,𝑁 c 𝑁 - 𝑇 𝛼 2 𝜌𝜅 𝑇 . (3.64)
We know the first term. The second term can be calculated independently from the simulations or from the model proposed by Long and Lequeux [START_REF] Long | Heterogeneous dynamics at the glass transition in van der waals liquids, in the bulk and in thin films[END_REF]. In a previous chapter we showed that the Long-Lequeux model provides a very good description of the 𝑇 dependence of 𝜌 and 𝜅 𝑇 . Therefore, we employ this model here. With this choice the right-hand side of Eq. (3.64) can be calculated analytically. The dashed line in Fig. 3.17

shows that the result of this calculation agrees with Eq. (3.63), as it should be for thermally equilibrated systems. Here this agreement illustrates again the internal consistency of the MD simulations.

Heat capacity of the polymer liquid: comparison with experiment

The inset of Fig. 3.17 compares the simulation results for 𝑐 𝑝 (circles with solid line) with the experimental heat capacity of liquid 1,4 PBD (dashed line) [START_REF] Grebowicz | Heat capacities of 1,4-polybutadienes[END_REF]. The CHAPTER 3. BULK SYSTEMS ANALYSIS experimental data refer to cis-PBD and trans-PBD. Both polymers crystallize, but have identical 𝑐 𝑝 in the liquid phase. To parametrize the 𝑇 dependence of liquid PBD Ref. [START_REF] Grebowicz | Heat capacities of 1,4-polybutadienes[END_REF] recommends the following equation for the specific heat per monomer: How can this difference be explained? A comprehensive discussion of the thermal analysis of polymers is provided in the book by Wunderlich [START_REF] Wunderlich | Thermal analysis of polymeric materials[END_REF]. Section 2.3.9 of [START_REF] Wunderlich | Thermal analysis of polymeric materials[END_REF] reports that the heat capacity of many liquid polymers increases linearly with temperature.

𝑐 𝑝 =
Eq. (3.65) thus represents a typical behavior found in experiment. The linear 𝑇 dependence is interpreted as resulting from the superposition of three contributions [START_REF] Pyda | Computation of heat capacities of liquid polymers[END_REF]:

𝐶 𝑝 = 𝐶 vib + 𝐶 conf + 𝑇 𝑉 𝛼 2 𝜅 𝑇 , (3.66)
where 𝐶 vib denotes the vibrational heat capacity at constant volume and 𝐶 conf the contribution to 𝐶 𝑝 due to conformational rotations of the chain backbone. The analysis for several hydrocarbon polymers in [START_REF] Pyda | Computation of heat capacities of liquid polymers[END_REF] shows that, even in the liquid, the major contribution to 𝐶 𝑝 comes from 𝐶 vib . As explained in [START_REF] Grebowicz | Heat capacities of 1,4-polybutadienes[END_REF][START_REF] Wunderlich | Thermal analysis of polymeric materials[END_REF][START_REF] Pyda | Computation of heat capacities of liquid polymers[END_REF], 𝐶 vib can be calculated from the vibrational spectrum of a polymer by considering two modes of vibration: "skeletal modes", involving torsion vibrations, bond-angle vibrations as well as collective vibrations along the chain backbone, and "group modes", involving vibrations of small groups of atoms, e.g. stretching vibrations of the CC and CH bonds, bending vibrations of the angle in CH 2 , etc. (see Fig. 2.48 on p. 123 in [START_REF] Wunderlich | Thermal analysis of polymeric materials[END_REF] and Table 4 in [START_REF] Grebowicz | Heat capacities of 1,4-polybutadienes[END_REF]). To calculate 𝐶 vib both modes are treated quantum mechanically, based on the Debye and Einstein theories for the heat capacity.

The Debye and Einstein theories predict the heat capacity to increase toward the classical limit with increasing 𝑇. Applications reveal that this classical limit is not yet reached in the experimentally studied 𝑇 regime of the polymer liquid (cf. [START_REF] Pyda | Computation of heat capacities of liquid polymers[END_REF] and p. 138 in [START_REF] Wunderlich | Thermal analysis of polymeric materials[END_REF]). Therefore, 𝐶 vib still increases with 𝑇. Although 𝐶 conf and 𝑇 𝑉 𝛼 2 /(𝜌𝜅 𝑇 ) are found to decrease upon heating, this decrease is weaker than the increase of 𝐶 vib . Therefore, the 3.5. ROUSE THEORY ANALYSIS 83 sum of the three contributions still leads to an increase of 𝐶 𝑝 with 𝑇, which can be fit by linear 𝑇 dependence.

Clearly, our classical united-atom model must fail to reproduce this behavior: Hydrogen atoms are not explicitly taken into account and quantum mechanical effects are absent. In our case, 𝑐 vib is a constant given by the classical result as a multiple of 𝑘 B .

If quantum mechanical skeletal and group vibrations are really the cause of the discrepancy between the simulated and experimental 𝑐 𝑝 , can we approximately correct for these missing effects? As an attempt to do so, we scanned the heat capacity data for the solid state from Fig. 9 of [START_REF] Grebowicz | Heat capacities of 1,4-polybutadienes[END_REF] and used them as a proxy for the vibrational contribution because it accounts for contributions to 𝑐 𝑝 from dihedral rotations and pair interactions in the polymer liquid.

(
The squares in the inset of Fig. 3.17 show that Eq. (3.67) is indeed reasonable. Using 𝑐 vib (𝑇 ) from experiments inverts the 𝑇 dependence originally found from the simulation, now yielding semiquantitative agreement with experiment (squares in the inset). Our findings are similar to those published recently by Honguy et al [START_REF] Gao | Comparing simulated specific heat of liquid polymers and oligomers to experiments[END_REF], where an agreement of the heat capacity extracted from MD simulations with the experimental data was achieved after taking into account the corrections for quantum effects.

Rouse theory analysis 3.5.1 Introduction

Being one of the simplest theoretical polymer models, the Rouse model [START_REF] Rouse | A theory of the linear viscoelastic properties of dilute solutions of coiling polymers[END_REF] is still extensively used for studying polymer systems in experiments (typically neutron scattering [START_REF] Richter | From rouse dynamics to local relaxation: a neutron spin echo study on polyisobutylene melts[END_REF][START_REF] Richter | Direct observation of the crossover from α-relaxation to rouse dynamics in a polymer melt[END_REF][START_REF] Riedel | Rouse-model-based description of the dielectric relaxation of nonentangled linear 1,4-cis-polyisoprene[END_REF]) and simulations [START_REF] Bennemann | Molecular-dynamics simulation of a glassy polymer melt: rouse model and cage effect[END_REF][START_REF] Krushev | The role of internal rotational barriers in polymer melt chain dynamics[END_REF][START_REF] Krushev | Intramolecular caging in polybutadiene due to rotational barriers[END_REF][START_REF] Mondello | Dynamics of n-alkanes: comparison to rouse model[END_REF][START_REF] Pérez-Aparicio | Chain dynamics of poly(ethylene-alt-propylene) melts by means of coarse-grained simulations based on atomistic molecular dynamics[END_REF][START_REF] Kalathi | Rouse mode analysis of chain relaxation in homopolymer melts[END_REF][START_REF] Hsu | Detailed analysis of rouse mode and dynamic scattering function of highly entangled polymer melts in equilibrium[END_REF]. The model considers a polymer chain as a sequence of beads connected by harmonic springs, that is surrounded by a (structureless)

Langevin thermostat imposing in addition to temperature 𝑇 also a friction coefficient (both parameters are coupled to a random force via the fluctuation-dissipation theorem) [START_REF] Doi | The theory of polymer dynamics[END_REF]. It is expected to provide a suitable description of the properties of unentangled polymer melts on length scales larger than the monomer size and on the time scales larger than the associated monomer relaxation time. Yet, since the excluded volume and chain stiffness are not taken into account, deviations from real polymers are also expected at small scales [START_REF] Mondello | Dynamics of n-alkanes: comparison to rouse model[END_REF].

Analytically the Rouse model provides a framework of describing the motion of a chain using the discrete Rouse modes ⃗ 𝑋 𝑝 (𝑡), defined as [5] The modes are theoretically considered to be orthogonal with each other [START_REF] Kreer | Monte carlo simulation of long chain polymer melts: crossover from rouse to reptation dynamics[END_REF], so that the properties of the polymer are given only by their autocorrelations 5

⃗ 𝑋 𝑝 (𝑡) = 1 𝑁 𝑁 ∑ 𝑛=1 ⃗ 𝑟 𝑛 (𝑡) cos [ (𝑛 -1 2 )𝑝𝜋 𝑁 ] (𝑝 = 0, … , 𝑁 -1) (3
𝐶 𝑝𝑝 (𝑡) = ⟨ ⃗ 𝑋 𝑝 (𝑡) ⃗ 𝑋 𝑝 (0)⟩ ⟨ ⃗ 𝑋 2 𝑝 (0)⟩ . ( 3.70)
The functions 𝐶 𝑝𝑝 (𝑡) allow direct observation of the chain relaxations on the corresponding length scales. Also, many static and dynamic observables, which are relevant in the experiments and simulations, can be expressed in terms of 𝐶 𝑝𝑝 , so by knowing them one can fully characterize the polymer within the framework of the Rouse model.

As for the PBD, the original CRC and FRC models have been already studied within the framework of the Rouse model in Refs. [START_REF] Krushev | The role of internal rotational barriers in polymer melt chain dynamics[END_REF][START_REF] Krushev | Intramolecular caging in polybutadiene due to rotational barriers[END_REF] with ⃗ 𝑟 𝑛 (𝑡) in Eq. (3.68) being the position of united atom 𝑛 at time 𝑡 of the PBD chain. In the following sections we compare 5 In practice the cross-correlations ⟨ ⃗ 𝑋 𝑝 ⃗ 𝑋 𝑞≠𝑝 ⟩ are of the order of 1% of the autocorrelations ⟨ ⃗ 𝑋 2 𝑝 ⟩, so they are usually neglected. Our calculations and the results reported in Ref. [START_REF] Krushev | The role of internal rotational barriers in polymer melt chain dynamics[END_REF] show that this is indeed the case for PBD. and extend this analysis using our models across wider range of temperatures. We also relate our results to the literature data for common polymers.

Static Rouse modes

The static Rouse modes autocorrelations are defined as

𝐶 𝑝𝑝 = ⟨ ⃗ 𝑋 𝑝 (0) ⃗ 𝑋 𝑝 (0)⟩ = ⟨ ⃗ 𝑋 2 𝑝 (0)⟩ . (3.71)
Ref. [START_REF] Krushev | The role of internal rotational barriers in polymer melt chain dynamics[END_REF] reported no noticeable differences in 𝐶 𝑝𝑝 between the models and temperatures up to 𝑝 = 40. Indeed, the results of our calculations appear to overlap with the data scanned from Fig. 3a in Ref. [START_REF] Krushev | The role of internal rotational barriers in polymer melt chain dynamics[END_REF], as presented in Fig. 3.18. The visible deviations appear for larger 𝑝.

Interestingly, for large 𝑝 one observes three spikes in the data, which are present for all models and temperatures. A closer look on the data (inset in Fig. 3.18) revealed, that

𝐶 𝑝𝑝 = 𝑏 2 e 8𝑁 1 sin 2 ( 𝜋𝑝 2𝑁 ) = 𝐶 ∞ 𝑙 2 8𝑁 1 sin 2 ( 𝜋𝑝 2𝑁 ) . (3.72)
However, due to the chain stiffness and excluded volume effects, only the first few modes follow such a scaling even in the bead-spring simulations [START_REF] Bennemann | Molecular-dynamics simulation of a glassy polymer melt: rouse model and cage effect[END_REF][START_REF] Kreer | Monte carlo simulation of long chain polymer melts: crossover from rouse to reptation dynamics[END_REF]. Although the corrections that account for these effects were proposed [START_REF] Meyer | Static rouse modes and related quantities: corrections to chain ideality in polymer melts[END_REF] and applied to the original CRC and FRC models of PBD in Ref. [START_REF] Krushev | The role of internal rotational barriers in polymer melt chain dynamics[END_REF], we are able to use our parametrization of the conformational

𝐶 𝑝𝑝 = 𝑙 2 8𝑁 [ 𝐶 ∞ sin 2 ( 𝜋𝑝 2𝑁 ) - 4𝐴𝛼 1+𝛼 1-𝛼 1 -2𝛼 cos( 𝜋𝑝 𝑁 ) + 𝛼 2 -4(𝛼 -𝐴𝛼) + 𝑂(𝑁 -1 )] . (3.73)
Fig. 3.19 tests our parametrization. The scaling of the axes is inspired by Eq. (3.72), so that for 𝑝 𝑁 → 0, the value of 2(𝜋𝑝) 2 𝑁 𝐶 𝑝𝑝 → 𝑙 2 𝐶 ∞ = 𝑏 2 e as it is shown in the inset. We present the data at 𝑇 = 353 K, but the result is valid for all studied temperatures. Using the values of the parameters reported in Appendix E, our predictions for the 𝐶 𝑝𝑝 are in a good agreement with the simulation data. The differences between the models are very small, since the parameters also do not differ significantly. This shows once more that the static properties of PBD on the scales from the size of the statistical segment to the size of the chain are mildly sensitive to the presence of the dihedrals. e for each model.

Dynamic Rouse modes

The dynamic Rouse modes autocorrelations are defined in Eq. (3.70). The Rouse model predicts them to decay exponentially [START_REF] Rouse | A theory of the linear viscoelastic properties of dilute solutions of coiling polymers[END_REF] 𝐶

𝑝𝑝 (𝑡) = exp [ - 𝑡 𝜏 𝑝 ] (3.74)
with 𝜏 𝑝 = 𝜏 R /𝑝 2 and 𝜏 R = 𝜏 𝑝=1 being the Rouse time, which is the longest relaxation time of the polymer chain corresponding to the mode of the largest length scale (𝑝 = 1).

However, it is known from extensive simulation studies, that due to the excluded volume, chain stiffness and "viscoelastic hydrodynamic interactions" [START_REF] Farago | Hydrodynamic and viscoelastic effects in polymer diffusion[END_REF][START_REF] Morhenn | Collective intermolecular motions dominate the picosecond dynamics of short polymer chains[END_REF], the 𝐶 𝑝𝑝 (𝑡) deviates from simple exponential decay in Eq. (3.74) (especially for large 𝑝) and, instead, can be well described by the stretched exponential (Kohlrausch-Williams-Watts function) [START_REF] Bennemann | Molecular-dynamics simulation of a glassy polymer melt: rouse model and cage effect[END_REF][START_REF] Pérez-Aparicio | Chain dynamics of poly(ethylene-alt-propylene) melts by means of coarse-grained simulations based on atomistic molecular dynamics[END_REF][START_REF] Kalathi | Rouse mode analysis of chain relaxation in homopolymer melts[END_REF][START_REF] Hsu | Detailed analysis of rouse mode and dynamic scattering function of highly entangled polymer melts in equilibrium[END_REF]]

𝐶 𝑝𝑝 (𝑡) = 𝐴 𝑝 exp [ -( 𝑡 𝜏 𝑝 ) 𝛽 𝑝 ] . (3.75)
In addition, for a polymer melt close to the glass transition a two-step relaxation was reported [START_REF] Bennemann | Molecular-dynamics simulation of a glassy polymer melt: rouse model and cage effect[END_REF]. Depending on the mode number, the emergence of a plateau between the initial short-time decay and the final long-time decay can be observed. Ref. [START_REF] Krushev | Intramolecular caging in polybutadiene due to rotational barriers[END_REF] documented that this is the case for the original CRC model of PBD. By contrast, the original FRC model showed faster, single stretched exponential relaxation. an increase of the relaxation times is typically observed for glass forming systems upon decreasing the temperature. It is usually explained in the bead-spring model simulations [START_REF] Bennemann | Molecular-dynamics simulation of a glassy polymer melt: rouse model and cage effect[END_REF] as a collective effect, where a particle becomes trapped within a cage created by its nearest neighbors. While such picture pertains for the FRC model, in case of the CRC(2,4) models, the trapping has been demonstrated to originate from the dihedral barriers ("intramolecular caging") [START_REF] Krushev | Intramolecular caging in polybutadiene due to rotational barriers[END_REF]: when a particle cannot overcome the barrier it participates in caging another particle.

In order to characterize the emergence of the plateaus and the stretching, we attempted to fit the 𝐶 𝑝𝑝 (𝑡) using Eq. (3.75). However, having 4 models, 7 working temperatures and 115 modes means 3220 fits to perform! Another problem is that a unique subset of the data should be used for each fit due to the differently emerging plateaus for each model.

In order to overcome these issues we use an automated fitting with the following approach:

• A typical fitting situation is presented in Fig. 3.22a. At first, the values of the relaxation times have been extracted previously in Fig. 3.21 as points where 𝐶 𝑝𝑝 (𝑡 = 𝜏 𝑝 ) = 1/𝑒 (horizontal dashed line in Fig. 3.22a). They were fixed in all the fits. The small 𝑝 modes that had not relaxed over a time of the simulation at low temperatures were omitted from the fit.

• Since there is a prefactor 𝐴 𝑝 , Eq. (3.75) has to be corrected, so that the definition of the relaxation times 𝐶 𝑝𝑝 (𝑡 = 𝜏 𝑝 ) = 1/𝑒 holds

𝐶 𝑝𝑝 (𝑡) = 𝐴 𝑝 exp [ -(ln 𝐴 𝑝 + 1)( 𝑡 𝜏 𝑝 ) 𝛽 𝑝 ] . (3.76)
• The values of 𝐶 𝑝𝑝 (𝑡) were limited in time from above at 𝑡 < 10𝜏 𝑝 in order to omit the noisy tails of the autocorrelation.

• For long times 𝑡 > 2 × 106 fs, where the data is spaced linearly, only the time points selected in logarithmic fashion were used in order to make the fit less biased by the tails, that are densely populated with data (black symbols in Fig. 3.22a).

• We are interested only in the second step of the relaxation. Since the initial relaxation appear to overlap for all our models, we can filter it out by considering only the points where the differences between the models start to occur. In At the same time, the stretching exponent 𝛽 𝑝 is weakly dependent on the temperature and decreases with the mode number 𝑝 starting from 𝛽 ≈ 0.8 down to 𝛽 ≈ 0.5, which means that the relaxation is slower than the exponential (the curves are stretched). While the decrease from 𝛽 ≈ 0.8 at small 𝑝 agrees with the results reported in the literature [START_REF] Pérez-Aparicio | Chain dynamics of poly(ethylene-alt-propylene) melts by means of coarse-grained simulations based on atomistic molecular dynamics[END_REF][START_REF] Kalathi | Rouse mode analysis of chain relaxation in homopolymer melts[END_REF], our data is qualitatively different at large 𝑝, where the bead-spring simulations report a minimum of 𝛽 ≈ 0.5 at 𝑁 /𝑝 ≈ 3.

In principle, for 𝑁 /𝑝 ≳ 10 it would possible to fix 𝛽 𝑝 ≈ 0.8 for all 𝑇 and obtain a fit of similar quality. This fixings implies time-temperature superposition, an important principle in glass physics and polymer physics. However, for 𝑁 /𝑝 ≲ 10, the principle For the CRC4 and the FRC models, the stretched exponential, does not actually represent well the shape of the relaxation at 𝑝 ≳ 10 (as it can be seen in Figs. 3.20c to 3.20f). In fact, the values of 𝑝, from which the fits fail, can be observed in Fig. 3.22b where the amplitudes 𝐴 𝑝 come back to unity after the initial decrease. This means that the fitting routine reached the upper bound in 𝐴 𝑝 , so the stretching parameter 𝛽 𝑝 was the only adjustable variable for the fit 7 . While for the CRC4 model 𝛽 𝑝 has a similar decreasing trend as of the CRC2 and CRC models, it is qualitatively different for the FRC model. In the case of the FRC model, the stretching exponent shows a minimum of 𝛽 𝑝 ≈ 0.5 at 𝑁 /𝑝 ≈ 3. Interestingly, similar minima are reported in the bead-spring model simulations [START_REF] Pérez-Aparicio | Chain dynamics of poly(ethylene-alt-propylene) melts by means of coarse-grained simulations based on atomistic molecular dynamics[END_REF][START_REF] Kalathi | Rouse mode analysis of chain relaxation in homopolymer melts[END_REF]. In general, the failure of the stretched exponential for the CRC4 3.6. MEAN-SQUARE DISPLACEMENTS 93 and FRC models can be explained by the fact that the studied region of temperatures 213 ≤ 𝑇 ≤ 353 K is quite far from the corresponding glass transition temperatures 𝑇 g . We can observe such a temperature dependence in Fig. 3.22b by the fact that the amplitudes 𝐴 𝑝 for these models reach the cap of unity earlier for higher temperatures.

Mean-square displacements

Having seen in the previous section, that the dihedral barriers in PBD lead to an enormous change in the relaxation times of the Rouse modes, it is important to investigate how they influence the motion of the united atoms. In the simulations this motion can be observed directly by looking at their mean-square displacements (MSD). In particular, it is interesting to look at

𝑔 0 (𝑡) = 1 𝑁 𝑁 ∑ 𝑛=1 ⟨ [ ⃗ 𝑟 𝑛 (𝑡) -⃗ 𝑟 𝑛 (0)] 2 ⟩ , (3.77)
which describes the MSD averaged over all united atoms of a chain;

𝑔 1 (𝑡) = ⟨ [ ⃗ 𝑟 𝑁/2 (𝑡) -⃗ 𝑟 𝑁/2 (0)] 2 ⟩ , (3.78)
which is the MSD of the innermost united atom of a chain;

𝑔 2 (𝑡) = 1 𝑁 𝑁 ∑ 𝑛=1 ⟨ [ ⃗ 𝑟 𝑛 (𝑡) -⃗ 𝑅 cm (𝑡) -⃗ 𝑟 𝑛 (0) + ⃗ 𝑅 cm (0)] 2 ⟩ , (3.79)
which describes an averaged MSD of the united atoms with respect to the chain's center of mass ⃗ 𝑅 cm (𝑡);

𝑔 3 (𝑡) = ⟨ [ ⃗ 𝑅 cm (𝑡) -⃗ 𝑅 cm (0)] 2 ⟩ , (3.80)
which is the MSD of the chain's center of mass itself;

𝑔 4 (𝑡) = ⟨ [ ⃗ 𝑟 end (𝑡) -⃗ 𝑟 end (0)] 2 ⟩ , (3.81)
which describes the MSD of the united atoms at the ends of the chain. temperature 𝑇 g , the MSD shows three distinctive regimes [2]. At short times, the initial motion of particles just after the beginning of the measurement is referred to as ballistic, since it is determined by Newtonian dynamics with 𝑔 0 (𝑡) ∼ 𝑡 2 . In the long times regime the particles diffuse with 𝑔 0 (𝑡) ∼ 𝑡. The intermediate regime is characteristic of polymer systems. It shows a behavior of 𝑔 0 (𝑡) ∼ 𝑡 𝑥 with the exponent 𝑥 < 1 and, therefore, is called subdiffusive. The value of 𝑥 is typically between 0.5 < 𝑥 < 0.75 8 . With lowering the temperature, the situation becomes slightly different. In particular, the ballistic regime is followed by an emerging plateau, that shifts the subdiffusive regime to longer times. Such a plateau is common for the glass forming systems: it indicates the trapping of the particles in a cage, that leads to the dynamical arrest.

Fig. 3.24 presents the MSD averaged over all united atoms of a chain 𝑔 0 (𝑡) for PBD simulated using all our models at different temperatures as empty colored symbols. The corresponding Rouse predictions are given by the colored lines with black dots

𝑔 R 0 (𝑡) = 𝑔 R 2 (𝑡) + 𝑔 3 (𝑡). (3.83)
The general behavior of the curves follows the description above with the regimes denoted in the figure. However, the data for different models shows significantly different temperature dependences. At our working temperatures 213 < 𝑇 < 353 K the FRC model demonstrates the fastest motion. The curves are very close to each other and have clearly distinct regimes. With introducing the dihedrals and making the barriers higher the dynamics becomes slower, so that the diffusive regime cannot be reached over few 𝜇s of the simulations. In the case of the CRC model we observe the caging plateaus that make a dramatic effect on the dynamics with lowering the temperatures. As it was mentioned in the previous section, the presence of the dihedrals imposes constraints on the movement of the particles when the thermal fluctuations are comparable to the barrier height. For PBD those are the 𝛼 and 𝛽 dihedrals. When a particle becomes trapped within the barriers, it also cannot anymore make room for the rearrangements of its neighbors.

Considering many such particles, this leads to the structural caging. Since the height of the dihedral barriers is larger than pair energy barriers, one observes the emergence of the cage at higher temperatures when comparing to the FRC model. We have performed t (fs) . Yet, this inconsistency rather reflects the poor quality of the extrapolation (cf. Fig. 3.27), than having a physical meaning. The VFT equation works well at the temperatures close to 𝑇 g , but reaching the relaxation times at these temperatures requires much longer simulations than few 𝜇s.

g 0 (Å 2 ) 113 K 153 K 213 K 225 K 240 K 253 K 273 K 293 K 353 K (d)
As an aside, we compared our results for the CRC model to various experimental references and original model simulations (Fig. 3.26). Our data fits within the experiments, confirming that our simulations represent the dynamics of PBD well in the studied temperature interval.

Discussion

By now we have adapted the united atom chemically realistic (CRC) model of 1,4polybutadiene to be able to perform the simulations of the bulk systems and systems with a free surface in LAMMPS code. We have performed extensive bulk simulation and studied The empty symbols represent the scanned data points, the solid lines are recreated VFT fits.

Ref. [START_REF] Smith | Relationship between the 𝛼-and 𝛽-relaxation processes in amorphous polymers: insight from atomistic molecular dynamics simulations of 1,4 polybutadiene melts and blends[END_REF] are the data from old simulations using the original CRC model. Refs. [START_REF] Arbe | Study of the dynamic structure factor in the 𝛽 relaxation regime of polybutadiene[END_REF][START_REF] Arbe | Merging of the 𝛼 and 𝛽 relaxations in polybutadiene: a neutron spin echo and dielectric study[END_REF] are the dielectric spectroscopy data. Ref. [START_REF] Berry | The viscosity of polymers and their concentrated solutions[END_REF] are the data based on the measurements of the monomeric friction coefficient. Ref. [START_REF] Solar | Relaxation processes and glass transition of confined polymer melts: a molecular dynamics simulation of 1,4-polybutadiene between graphite walls[END_REF] are dielecric relaxation data obtained from the recent simulations using the original CRC model. Ref. [START_REF] Schmidtke | Temperature dependence of the segmental relaxation time of polymers revisited[END_REF] are the combined data from dielectric spectroscopy (low 𝑇) and NMR (high 𝑇) for PBD of 𝑀 = 1450 g/mol and 𝑀 = 2020 g/mol (in our simulations 𝑀 = 1569 g/mol).

the structure and dynamics of PBD over a wide range of temperatures 213 K < 𝑇 < 353 K.

Our results fit within the experimental references for PBD of similar composition and molecular weight.

We have also studied the influence of the dihedral barriers on the properties of PBD by performing the simulations with the reduced (CRC2 and CRC4 models) and completely disabled barriers (FRC model). We observed minor changes in the structure of the PBD melt with modified/disabled dihedrals, which could not be seen in the similar studies carried out earlier [START_REF] Krushev | The role of internal rotational barriers in polymer melt chain dynamics[END_REF][START_REF] Krushev | Intramolecular caging in polybutadiene due to rotational barriers[END_REF]. By contrast, the dynamics of PBD is orders of magnitude faster for these models, which also leads to significant shifts of the glass transition temperature 𝑇 g to lower values. By the analysis of the local relaxation of the united atoms, we attribute these differences in the dynamics to the fact, that the trapping of the UAs between the dihedral barriers shifts the onset of the structural caging to higher temperatures.

Our results challenge the theories that try to predict the dynamics of the glass forming systems from their structure. In particular, the mode coupling theory (MCT) has been extensively discussed in the previous works [START_REF] Krushev | The role of internal rotational barriers in polymer melt chain dynamics[END_REF][START_REF] Krushev | Intramolecular caging in polybutadiene due to rotational barriers[END_REF][START_REF] Paul | Glass transition in 1,4-polybutadiene: mode-coupling theory analysis of molecular dynamics simulations using a chemically realistic model[END_REF]. Below we discuss the more recent elastically collective nonlinear Langevin equation (ECNLE) theory [START_REF] Mirigian | Elastically cooperative activated barrier hopping theory of relaxation in viscous fluids: II. Thermal liquids[END_REF][START_REF] Mirigian | Dynamical theory of segmental relaxation and emergent elasticity in supercooled polymer melts[END_REF][START_REF] Mirigian | Elastically cooperative activated barrier hopping theory of relaxation in viscous fluids. i. general formulation and application to hard sphere fluids[END_REF].

The theory considers a polymer melt as a system of hard spheres at volume fractions determined from the temperature dependence of the dimensionless compressibility. An energetic dynamic barrier can be calculated from the structure of the system, so that the time of crossing describes the local dynamics. In order to account for the collective rearrangements, the additional elastic barrier is introduced (which is also calculated using the dynamic free energy). Since the theory is non-linear, we expect the small differences in structure to cause significant differences in the dynamics. The theory was successfully applied to predict the segmental relaxation, glass transition temperature and fragility of many polymers, including PBD.

We have put a lot of effort into describing the structure of PBD simulated with our models and indeed observed small differences in the structural properties between the models. The two important quantities here are the characteristic ratio 𝐶 ∞ and the dimensionless compressibility 𝑘 B 𝑇 𝜌𝜅 𝑇 . The ECNLE theory assumes 𝐶 ∞ to be temperature independent and uses it to map the compressibility of the melt to the one of the hard spheres [START_REF] Mirigian | Dynamical theory of segmental relaxation and emergent elasticity in supercooled polymer melts[END_REF]. We have implemented the numerical solution of the ECNLE theory following

Refs. [START_REF] Mirigian | Elastically cooperative activated barrier hopping theory of relaxation in viscous fluids: II. Thermal liquids[END_REF][START_REF] Mirigian | Dynamical theory of segmental relaxation and emergent elasticity in supercooled polymer melts[END_REF][START_REF] Mirigian | Elastically cooperative activated barrier hopping theory of relaxation in viscous fluids. i. general formulation and application to hard sphere fluids[END_REF] and performed the calculations of the relaxation times plugging CHAPTER 3. BULK SYSTEMS ANALYSIS into the theory the 𝐶 ∞ at 𝑇 = 353 K from Section 3.1.2 and the linearized in 1/𝑇 dimensionless compressibility using the Long-Lequeux fits from Section 3.3.3 for our CRC and FRC models. Fig. 3.27 present the results of the calculations of the relaxation times, that correspond to the hopping across the dynamic barrier only (without the additional contribution from the elastic barrier). Even though the differences in the structure between the models were taken into account in the calculations, the data overlaps. Interestingly, for both models, the dynamics at high temperature (where the barrier is very small) we observe an agreement with the MD data of the FRC model. However, with decreasing the temperature, even without the elastic contribution the ECNLE theory already overestimates the relaxation times of the FRC model. This means that in the studied range, the mapping of temperature to the volume fraction via the renormalized dimensionless compressibility is not valid for the FRC model. At the moment the reasons for this disagreement are not clear and further research is required.This is the main reason why we decided not to expose the ECNLE theory, its implementation, and the details of its application to our simulation data in this thesis (despite the large amount of work spent on this project). Chapter 4

Film systems analysis

This chapter presents the results obtained from the simulations of the film systems. We begin our analysis with determining the thickness and studying its temperature dependence upon cooling in Section 4.1. Then in Section 4.2 we relate the obtained results to the layer-resolved investigation of the dynamical properties. Finally, we discuss our findings in Section 4.3.

Density profiles and thickness

Let the 𝑧-axis be perpendicular to the film plane (Fig. 2.8). Then the density profile along the 𝑧 direction can be defined as

𝜌(𝑧) = 1 𝐴 ⟨ 𝑁𝑁 c ∑ 𝑖=1 𝛿(𝑧 -𝑧 𝑖 )⟩ , (4.1)
where 𝐴 is the area of the cross-section of the simulation box in the 𝑥𝑦 plane. 1 In practice, such a calculation is performed by introducing finite (yet, small enough in thickness) slabs along the 𝑧-direction and counting how many particles are in each slab volume. Fig. 4.1

presents the density profiles of our films simulated using the CRC and FRC models at different temperatures. By looking at the shape of the curves, we can get an idea on the influence of some surface effects. In particular, the peaks on the left-hand side of the figures suggest that the presence of the wall leads to layering of united atoms next to it, which is typically found in liquids in contact with an impenetrable substrate [START_REF] Baschnagel | Computer simulations of supercooled polymer melts in the bulk and in confined geometry[END_REF].

CHAPTER 4. FILM SYSTEMS ANALYSIS

We compared this part of our density profiles to the results of Ref. [START_REF] Solar | Relaxation processes and glass transition of confined polymer melts: a molecular dynamics simulation of 1,4-polybutadiene between graphite walls[END_REF], where the PBD was simulated in confinement between two graphite walls using the original CRC model (dashed black line in Fig. 4.1a). An excellent agreement was obtained, which is expected considering that our modified CRC model was tuned to reproduce the same density as the original one.

At the free surface, one observes a gradual decrease of the density down to zero, which is commonly attributed to the presence of capillary waves [START_REF] Peter | Structure et relaxation structurale des fondus de polymères vitrifiables en couches minces[END_REF]. The width of the decrease region is determined by the wave amplitudes, which are temperature dependent and limited by the system size. The plateau in the middle of the films corresponds to the value of the bulk density, as marked by the horizontal dotted lines in Fig. 4.1. From the analysis of the bulk systems, we know that they have almost identical densities, so the plateaus for the CRC and the FRC models overlap. Interestingly, they also overlap at the wall and at the free surface, meaning that surface effects on the density appear to be mildly influenced by the dihedrals.

Measuring the density profile of the film allows us to determine the thickness of the films using the Gibbs dividing surface (GDS) method [START_REF] Hansen | Theory of simple liquids[END_REF]. In order to do this, at first, the mean density ρ at the center of the film 𝑧 c is calculated by averaging 𝜌(𝑧) over a wide enough region 𝑧 c -Δ𝑧/2 ≤ 𝑧 ≤ 𝑧 c + Δ𝑧/2. The choice of the center point 𝑧 c and the averaging region is arbitrary as long as the mean density can be well defined. In our case, the procedure was automatic with 𝑧 c being the middle of the non-zero data of the density profile and Δ𝑧 = 30 Å. The positions of the upper 𝑧 G + and lower 𝑧 G -Gibbs dividing surfaces are calculated as

𝑧 G ± = 𝑧 c + 1 ρ ∫ ±∞ 𝑧 c 𝜌(𝑧)d𝑧 . (4.2)
Qualitatively, the integration in Eq. z (Å) The density profiles of the films obtained from the simulations using the CRC (red) and FRC (blue) models at different temperatures. The dashed black line is the data scanned from Ref. [START_REF] Solar | Relaxation processes and glass transition of confined polymer melts: a molecular dynamics simulation of 1,4-polybutadiene between graphite walls[END_REF], where the PBD film confined between graphite was simulated using the original CRC model at 𝑇 = 353 K. The vertical dot-dashed lines are the positions of the graphite layers. The vertical dotted lines are the positions of the Gibbs dividing surfaces for the CRC model. The horizontal dotted lines correspond to the density of the bulk systems at the same temperatures.

ρ(z) (Å -3 ) CRC FRC (d)
principle, the 𝑧 G -also shifts with decreasing the temperature, since the first peak (that contributes the most to the integration in Eq. (4.2)) increases in height. However, such an influence is actually quite small (∼ 0.1 Å), so it is ignored.

Using the GDS method, the thickness can be monitored not only in the equilibrium simulations, but also upon cooling the films. Fig. 4.2 presents the dependence of thickness on temperature ℎ(𝑇 ) obtained from the cooling simulations at the rate Γ = 0.83 K/ns (colored circles) 2 and from the equilibrium simulations (colored crosses). Since the film system is much bigger than the bulk, the averaging statistics is better and small differences 104 CHAPTER 4. FILM SYSTEMS ANALYSIS in density are now more visible: the thickness of the CRC model is systematically slightly lower than for the FRC model at high temperature.

As it is usually observed in ellipsometry experiments [START_REF] Alcoutlabi | Effects of confinement on material behaviour at the nanometre size scale[END_REF], ℎ(𝑇 ) changes slope when the system is undergoing the glass transition. Yet, our cooling rate Γ = 0.83 K/ns is much faster than typically employed in the experiments, which leads to a broad transition region.

We have extracted the 𝑇 g from the cooling data using a similar approach as during the analysis of the bulk densities. The thickness ℎ(𝑇 ) was fitted with:

ℎ(𝑇 ) = 𝑤 ( 𝑀 -𝐺 2 ) ln [cosh ( 𝑇 -𝑇 g 𝑤 )] + (𝑇 -𝑇 g ) ( 𝑀 + 𝐺 2 ) + ℎ g , (4.4)
with 𝑤 being the width of the glass transition region, 𝑀 and 𝐺 the slopes dℎ(𝑇 )/d𝑇 in the melt and glass regions respectively, 𝑇 g the glass transition temperature and ℎ g the thickness of the film at 𝑇 = 𝑇 g . In order to account for the dependence on the size of the fitting region (𝑇 * g -Δ𝑇 , 𝑇 * g + Δ𝑇 ), we have done the fits over a set of Δ𝑇 ranging from tens of Kelvins to the highest possible value. Comparing to the bulk density data, the data for the film thickness have lower quality. Thus, the influence of the fitting region is more pronounced in the values of 𝑀 and 𝐺, while the 𝑇 g and ℎ g are much less sensitive to the choice of Δ𝑇 (Fig. 4.3). The final iteration of the fits with 𝑇 g , ℎ g , 𝑀 and 𝐺 fixed to The temperature dependence of the film thickness obtained by the GDS method from the simulations using the CRC (red) and the FRC (blue) models. The circles represent the data obtained from the cooling runs, and the crosses are the averaged values from the constant temperature simulations. The solid black lines are the fits using Eq. (4.4) with the parameters from Table 4.1.

their values at the highest possible Δ𝑇 leads to the values presented in Table 4.1 and the fitting curves represented by black solid lines in Fig. 4.2. We obtained an exact match of the glass transition temperature 𝑇 g = 43 K for the FRC model with the bulk simulations.

In case of the CRC model, the resulting value 𝑇 film g = 206 K is slightly higher than the bulk value 𝑇 bulk g = 192 K. Still, the film thickness data in the transition region is quite noisy even after smoothing, so it is hard to argue about reasons of possible differences.

Moreover, refitting the CRC data with the 𝑇 g fixed to the bulk value resulted in the fitting curve of the same quality (not shown) with comparable values of the parameters. We consider them in agreement with each other.

Our 𝑇 g values are, of course, the film averaged values, that result from the contribution of the bulk-like region in the middle of the film and the surface effects. The fact that they are basically the same as in the bulk for both models may suggest that either the surface effects are very weak, or the influence of the graphite wall and the free surface average out.

∆T (K) M ,G (Å/K) M G (d)
We have learned from the analysis of the bulk systems, that disabling the dihedrals leads to significant speed up of the dynamics of PBD while preserving its structure, which explains the presence of the shift of the 𝑇 g to lower values. It is interesting now to investigate the local dynamics of PBD across the film in order to understand the influence of the dihedral barriers on the surface effects, and consequently, on the 𝑇 g .

Mean-square displacements and dielectric relaxation

Our study of the dynamics in the films utilizes two methods. The first one is essentially the same as employed for the bulk systems: the simulations allow direct observation of the mean-square displacements of united atoms 𝑔 0 (𝑡) on different time scales. The second mimics the dielectric spectroscopy technique commonly employed in experiments (see [START_REF] Kremer | Recent advances in broadband dielectric spectroscopy[END_REF] and references therein), where the dynamics is analyzed by measuring the dielectric response of the system to the applied electric field alternating with different frequencies. It is inspired by the previous works [START_REF] Solar | The dielectric axation in polymer films: a comparison between experiments and atomistic simulations[END_REF][START_REF] Solar | Dielectric -relaxation of 1,4-polybutadiene confined between graphite walls[END_REF][START_REF] Smith | A molecular-dynamics simulation study of dielectric relaxation in a 1,4-polybutadiene melt[END_REF], where dielectric relaxation of PBD was studied in the simulations using the original CRC model. The fluctuation-dissipation theorem relates the experimentally relevant complex dielectric permittivity 𝜖 * to the correlation of the total dipole moment of the system ⃗ 𝑀, which was shown to be equal to the sum of independent chain contributions, which, in turn, are the sums of the independent dipoles of the 𝑐𝑖𝑠 groups of PBD3 (Fig. 4.4) [START_REF] Smith | A molecular-dynamics simulation study of dielectric relaxation in a 1,4-polybutadiene melt[END_REF][START_REF] Bahar | Segmental motions of cis-polyisoprene in the bulk state: interpretation of dielectric relaxation data[END_REF]. Thus, the dielectric spectroscopy experiments on PBD probe the relaxation on the segmental scale. The studied sample is placed between the electrodes, that act like walls confining the film, so the relevant quantity extracted from the simulation trajectories is the autocorrelation of the perpendicular to the film component of the dipole moment (DACF) of the 𝑐𝑖𝑠 groups of the PBD:

𝐶 𝑧 (𝑡) = ⟨𝑀 𝑧 (𝑡)𝑀 𝑧 (0)⟩ ⟨𝑀 𝑧 (0)𝑀 𝑧 (0)⟩ . (4.5)
This is an important limitation, since the relaxations in different directions may be different and it is the averaged dynamics that matters for the glass transition in the films. Thus, in our simulations we also study the autocorrelation of in-plane components 𝐶 𝑥,𝑦,𝑥𝑦 (𝑡) as well as the total autocorrelation of a dipole moment vector

𝐶(𝑡) = ⟨ ⃗ 𝑀 (𝑡) ⃗ 𝑀 (0)⟩ ⟨ ⃗ 𝑀 2 (0)⟩ . (4.6)
At first, Fig. 4.5 presents the comparison of the total mean-square displacements of united atoms 𝑔 0 (𝑡) and the total DACF 𝐶(𝑡) obtained from the bulk simulations (black solid lines) and from the film simulations (colored solid lines) at selected temperatures averaged over the whole simulation box. The in-plane and perpendicular components will be discussed below. In case of the CRC model, the averaging over the correlations of the 𝑐𝑖𝑠 dipole moments was performed for 𝐶(𝑡). In case of the FRC model, there is no distinction between the 𝑐𝑖𝑠 and 𝑡𝑟𝑎𝑛𝑠 monomers, because all the dihedral potentials are disabled. Therefore, one expects the dipole moment of each monomer to be equal to zero.

However, due to the excluded volume, there is still a preferred range of the dihedral angles (see Section 3.1.1) which leads to a non-zero value of ⃗ 𝑀. Thus, the DACF for the FRC model was obtained by averaging the correlations of the dipole moments of all monomers.

We use the time domain representation in order to be able to relate the relaxation regimes with the regimes of the MSD, and because performing the Fourier-transform of the data would require either a different approach to the spacing of the simulation trajectories, or enormous interpolation efforts.

The film averaged data of the MSD are slightly faster than the bulk curves at the corresponding temperatures. The DACF curves appear to decay slightly slower at long times in the films. Refs. [START_REF] Solar | The dielectric axation in polymer films: a comparison between experiments and atomistic simulations[END_REF][START_REF] Yelash | Three-step decay of time correlations at polymersolid interfaces[END_REF] attributed such additional slowing down to the adsorption-desorption kinetics of the chains at the graphite wall.

We are interested in the time scales of the structural relaxation. For the bulk systems, they were extracted as points, where a united atom has displaced over a distance of its own size 𝑔(𝑡 = 𝜏 ) = σ 2 . While for the MSD this definition is consistent between the CRC and FRC models, the same relaxation times correspond to different regimes in the dielectric relaxation between the models. Omitting the short time region of the DACF 𝑡 < 100 fs, the data for the CRC model are single decays, which allows to make a rough approximation of the structural relaxation times to the decay of 𝐶(𝑡 = 𝜏 ) ≈ 0.2. Meanwhile, the FRC model exhibits two-step relaxations and the structural relaxation times correspond to the second step with 𝐶(𝑡 = 𝜏 ) ≈ 0.01. Of course, the two discussed quantities cannot be directly related: the MSD represents the motion of individual united atoms, whereas the DACF describes the reorientation of a group comprised of several united atoms. In addition, in the simulations using the realistic model each individual dipole persists at all times, so only their orientations contribute the most to the correlation. In case of the FRC model, the non-zero dipole moment is the result of the average excluded volume effect, so the reorientation of the united atoms within the monomer are of importance, instead of the reorientation of the monomer as a whole.

Having compared the overall film dynamics to the bulk, it is interesting now to look at the relaxations in the different regions of the films. To do so, we have divided the films into slabs and calculated the MSD and the DACF in each of them. The layout of the slabs is presented in the Fig. there. An excellent agreement is obtained between the 𝐶 𝑧 (𝑡) in the middle of the film and between our 2nd thin slab from the wall and the 1st thick slab from the work [START_REF] Solar | Dielectric -relaxation of 1,4-polybutadiene confined between graphite walls[END_REF].

As it has been already mentioned, in the dielectric spectroscopy experiments, the electrodes act as walls for the film, so it is usually the perpendicular component of the dielectric relaxation that is measured. It is interesting to study if there are preferred orientations of the dipoles in the film slabs at the wall and at the free surface, since the experiments would not see the contributions from the dipoles that are oriented in-plane. decreasing the temperature, while the data for the FRC model appears to be independent of temperature. At high temperatures, all projections for each model are identical within the error bar of the measurement. There are small hints of a polarization in 𝑧 direction with preservation of the total magnitude at the wall, and an increase of the magnitude without pronounced polarization at the free surface. Unfortunately, with lowering the temperature our data suffers from insufficient statistics (or equilibration). While in the middle of the film all three projections appear to be identical (again within the error bar of the calculation), it is not possible to interpret physically the influence of the confinement.

Similarly to the bulk analysis, we have extracted the structural relaxation times from the MSD data in each slab in the parallel and perpendicular to the film directions as points where 𝑔 𝑥𝑦 0 (𝑡 = 𝜏 ) = 2 Comparing the data at the same temperatures between the CRC and FRC models we observe a similar behavior to what we have seen in the bulk: the dihedrals set the total order of magnitude of the relaxation time, which is then getting modified by the surface effects. The modifications become dramatic with decreasing the temperature. Yet, visually they penetrate no more than ∼ 20 Å inside the film. The slowing down at the wall and t (fs) t (fs) the speed up at the free surface must be averaging out, so the overall film dynamics (as we reported above) is quite similar to the bulk. Which also explains similar glass transition temperatures to the ones obtained for the bulk systems.

g xy 0 (Å 2 ) (a)
g z 0 (Å 2 ) (a)
Let us focus on the relaxation time gradients at the wall and at the free surface of the films. Figs. 4.11 and 4.12 show these data separating between relaxations in 𝑥𝑦 and 𝑧 directions and Fig. 4.13 presents a superposition of the relaxation time profiles and the density profiles. Additionally, we have fitted the data with a commonly used empirical "double-exponential" relation (see [START_REF] Schweizer | Progress towards a phenomenological picture and theoretical understanding of glassy dynamics and vitrification near interfaces and under nanoconfinement[END_REF][START_REF] Phan | Influence of longer range transfer of vapor interface modified caging constraints on the spatially heterogeneous dynamics of glass-forming liquids[END_REF] and references therein)

𝜏 (𝑧, 𝑇 ) = 𝜏 b (𝑇 ) exp [-𝐴(𝑇 ) exp (- 𝑧 𝜉 )] , (4.7)
with 𝜏 b (𝑇 ) being the relaxation time in the middle of the film at temperature 𝑇, 5 and 𝐴, 𝜉 are adjustable parameters. The results of the fits at selected temperatures are presented in Fig. 4.13 as colored lines and the values of the fitting parameters are plotted in the should be the same between the two models. Such scaling has already been observed before in [START_REF] Phan | Influence of longer range transfer of vapor interface modified caging constraints on the spatially heterogeneous dynamics of glass-forming liquids[END_REF]. Figs. 4.14a and 4.14b show that this is the case for the relaxation in the 𝑥𝑦 direction. In the 𝑧 direction the dependencies of the fit parameters are slightly different between the models, however these fits are definitely influenced by the closest points to the upper GDS which deviate from the overall behavior of the data (especially for the FRC model; see e.g. Figs. 4.10 and 4.11b) 6 . Furthermore, it is important to note that the gradients at the free surface increase with lowering the temperature as it is demonstrated in Figs. [START_REF] De Gennes | Scaling concepts in polymer physics[END_REF].14e and 4.14f, which is exactly the opposite than the behavior of the density profiles. As it has been already mentioned in Section 4.1 a gradual decrease of the density down to zero at the free surface is commonly attributed to the presence of capillary waves [START_REF] Peter | Structure et relaxation structurale des fondus de polymères vitrifiables en couches minces[END_REF], that emerge from a competition between the thermal energy and the surface tension. Their height, and therefore the width of the region of the decrease of the density is proportional to temperature. Moreover, Fig. 4.13 shows that this width is significantly smaller than the length over which the relaxation times are modified. Similar results are reported from simulations of supported films using a bead-spring polymer model [START_REF] Peter | Thickness-dependent reduction of the glass-transition temperature in thin polymer films with a free surface[END_REF] (Fig. 4.15). All these observations suggest that at the free surface the profile of the dynamics and the density profile cannot be related. Especially considering the fact that the decrease of the density profile next to the upper GDS results from the presence of vacuum in between the waves. The local density nearby an arbitrary particle can be in fact not very different from the bulk. Assuming that this is the case, then there is no change in properties of the melt with respect to the bulk-like middle of the film that could be responsible for the change in the dynamics. The existence of the gradient of the relaxation times can be attributed to the presence of the free surface itself and the fact that the cage that surrounds a given particle is open, which significantly speeds up the relaxation. This effect gradually damps towards the interior of the film.

• The situation at the wall side of the film is ambiguous. At a first glance, it seems from Fig. 4.13 that the relaxation times data follows the envelope of the density profiles.

Since the density profiles are identical between the models at the same temperature, one would expect the same scaling of the fit parameters with temperature. However, that is important to note comes from the comparison of the gradients in the 𝑥𝑦 direction with previously mentioned results from the simulations using a coarse- The profile of the relaxation times obtained from Ref. [START_REF] Peter | Thickness-dependent reduction of the glass-transition temperature in thin polymer films with a free surface[END_REF] where the polymer film supported on a structureless wall was simulated using a bead-spring polymer model.

grained model (Fig. 4.15). In particular, both films exhibit the layering of the particles next to the attractive wall that is visible on the density profiles. However, the data from Ref. [START_REF] Peter | Thickness-dependent reduction of the glass-transition temperature in thin polymer films with a free surface[END_REF] show a speed up of the relaxation at the wall, not the slowing down as in our case for both models. This disagreement can be explained by the differences in the properties of the substrate. In our case, the substrate is modeled by particles in a crystal lattice that make the wall very attractive and rough. In the work [START_REF] Peter | Thickness-dependent reduction of the glass-transition temperature in thin polymer films with a free surface[END_REF] the substrate is structureless (smooth) and is modeled using a weakly attractive potential that acts only in 𝑧 direction. Work [START_REF] Hanakata | Interfacial mobility scale determines the scale of collective motion and relaxation rate in polymer films[END_REF] studied the influence of the wall roughness on the dynamics at the substrate also by means of simulations using a bead-spring polymer model for the film and reported that it is sped up at the smooth walls and slowed down at the rough walls. As for our results, it remains unclear how this effect is influenced by the presence of the dihedrals.

Discussion

We have performed molecular dynamics simulations of the polybutadiene films supported by a graphite wall using the chemically realistic (CRC) united-atom model and the freely rotating (FRC) model, where all the dihedral potentials were disabled. We observed no significant differences in the density profiles of the film simulated using both models.

The glass transition temperatures extracted from the dependence of the film thickness on temperature upon cooling agree with the values obtained from the bulk systems.

The overall film dynamics observed in the mean-square displacements of united atoms and the dielectric relaxation is mildly different from the bulk dynamics. The structural The layer resolved analysis of the structural relaxation times extracted from the MSD data showed significant influence of the interfaces with respect to the middle of the film.

However, the effect of the wall and the free surface qualitatively appear to average out when contributing to the total 𝑇 g of the film. Analysis of the relaxation times gradients at the interfaces showed that at the free surface they nearly identical in the 𝑥𝑦 direction between the CRC and FRC models at the temperatures with the same relaxation times in the bulk-like middle of the film. The numerical results for the 𝑧 direction at the free surface and both directions at the wall are not conclusive due to difficulties of interpreting the behavior of the gradient decay length 𝜉. A visual inspection suggests that these gradients are different from each other.

Our data for the CRC model agree with the results reported in the previous studies of confined PBD [START_REF] Solar | The dielectric axation in polymer films: a comparison between experiments and atomistic simulations[END_REF][START_REF] Solar | Dielectric -relaxation of 1,4-polybutadiene confined between graphite walls[END_REF][START_REF] Solar | Relaxation processes and glass transition of confined polymer melts: a molecular dynamics simulation of 1,4-polybutadiene between graphite walls[END_REF] in the fact, that the dynamics is modified in a region of about 1 nm next to the graphite wall. We report slightly deeper perturbations at the free surface. Refs. [START_REF] Solar | The dielectric axation in polymer films: a comparison between experiments and atomistic simulations[END_REF][START_REF] Solar | Dielectric -relaxation of 1,4-polybutadiene confined between graphite walls[END_REF][START_REF] Solar | Relaxation processes and glass transition of confined polymer melts: a molecular dynamics simulation of 1,4-polybutadiene between graphite walls[END_REF] argued that the perturbation of around 1 nm into the film are the result of the interplay between the dihedral barriers present in the realistic model and the density changes created by the interfaces. They should lead to much smaller shifts in the glass transition temperature in the confined geometry with respect to what is typically observed in the widely used coarse-grained bead-spring polymer models (which do not have the dihedrals). Our FRC model mimics the bead-spring models, but it has a realistic length scale and shorted bond length. The comparison between the CRC and FRC models agrees with this hypothesis for the relaxation in the perpendicular to the film direction (which was the only one discussed in Refs. [START_REF] Solar | The dielectric axation in polymer films: a comparison between experiments and atomistic simulations[END_REF][START_REF] Solar | Dielectric -relaxation of 1,4-polybutadiene confined between graphite walls[END_REF][START_REF] Solar | Relaxation processes and glass transition of confined polymer melts: a molecular dynamics simulation of 1,4-polybutadiene between graphite walls[END_REF] and only for the films confined between walls), since there were no similarities found in the gradients between the models.

However, we report a different behavior of the relaxation time gradients in the parallel to the film direction, where the identical gradients were obtained between the CRC and FRC films, that have the same order of magnitude of the relaxation time in the middle. Chapter 5

Conclusions and Outlook

We now give a summary of the key results of our work and provide an answer to the title of the thesis "Is the interfacial mobility in glass-froming polymer films determined by collective motion or intramolecular energetics?"

In order to study the influence of realistic torsional constraints on the properties of glass forming polymers we have adapted the united atom model of 1,4-polybutadiene (PBD)

𝑐𝑖𝑠-𝑡𝑟𝑎𝑛𝑠 random copolymer (developed previously from quantum chemistry calculations and validated against experiments [START_REF] Smith | Molecular dynamics of a 1,4-polybutadiene melt. comparison of experiment and simulation[END_REF][START_REF] Smith | NMR spin-lattice relaxation and conformational dynamics in a 1,4-polybutadiene melt[END_REF][START_REF] Smith | A molecular-dynamics simulation study of dielectric relaxation in a 1,4-polybutadiene melt[END_REF]) for classical molecular dynamics simulations in the LAMMPS code [START_REF] Allen | Computer simulation of liquids[END_REF][START_REF] Lammps | Large-scale Atomic/Molecular Massively Parallel Simulator[END_REF][START_REF] Plimpton | Fast parallel algorithms for short-range molecular dynamics[END_REF]. Four versions of the model were introduced varying the scale of all torsional (dihedral) potentials in the system:

• Chemically realistic chain (CRC) model with the unmodified dihedrals potentials (which provide the distinction between the 𝑐𝑖𝑠 and 𝑡𝑟𝑎𝑛𝑠 monomers).

• CRC2 and CRC4 models with the dihedral potentials reduced by a factor of 2 and 4 correspondingly.

• Freely rotating chain (FRC) model with all dihedral potentials switched off. Here the distinction between the 𝑐𝑖𝑠 and 𝑡𝑟𝑎𝑛𝑠 groups is lost.

Such an interpolation provides a crossover to the widely used generic bead-spring polymer models, in which the realistic chemical details (including the dihedrals) are coarse-grained. Our FRC model resembles the bead-spring models with an important exception, that it has much shorter bond length with respect to the size of the particles.

The most used bead-spring model of [START_REF] Bennemann | Molecular-dynamics simulation of a glassy polymer melt: rouse model and cage effect[END_REF][START_REF] Peter | Thickness-dependent reduction of the glass-transition temperature in thin polymer films with a free surface[END_REF] type has a single length scale that determines both the size of the beads and the length of the bonds between them.

Our simulations were performed in the bulk using all 4 models and in the confined Note that in order to simulate supported films (including a free surface), it is not possible to use tail corrections for the pair-potential. The parameters of the original model had to be adapted to yield the same equation of state with a finite range pair-potential.

The extensive chapter on bulk results is on the one hand a validation of these modified model parameters, and on the other hand a new look at the PBD results. Indeed, today's increased computer power enables us to obtain much longer trajectories.

Bulk systems

The simulations of the bulk system revealed that reducing/disabling the dihedral potentials leads to little change in the conformational properties of PBD (discussed in Section 3.1.2). The analysis of the characteristic ratio 𝐶 ∞ (Fig. 5.1) agrees with the results of previous works [START_REF] Krushev | The role of internal rotational barriers in polymer melt chain dynamics[END_REF][START_REF] Krushev | Intramolecular caging in polybutadiene due to rotational barriers[END_REF] and its temperature dependence shows a similar (though slightly larger) chain extension upon cooling with respect to the experimental results on PBD of comparable composition [START_REF] Mays | Temperature dependence of unperturbed dimensions for stereoirregular 1,4-polybutadiene and poly(𝛼-methylstyrene)[END_REF]. The fact that the characteristic ratio of the FRC is almost identical to the CRC is certainly specific to PBD and makes this polymer an ideal model system to study the influence of the dihedral barriers on the dynamics, 𝑤𝑖𝑡ℎ𝑜𝑢𝑡 𝑐ℎ𝑎𝑛𝑔𝑖𝑛𝑔 𝑡ℎ𝑒 𝑠𝑡𝑟𝑢𝑐𝑡𝑢𝑟𝑒 𝑜𝑓 𝑡ℎ𝑒 𝑚𝑒𝑙𝑡.

The study of the temperature dependence of the bulk density (Section 3.2) of PBD showed that the dihedrals have a negligible influence in the high temperature region (Fig. 5.2). For a small bulk system the densities at 𝑇 ≳ 213 K are basically indistinguishable between the models. However, they have a significant influence on the glass transition temperature 𝑇 g . The Dalnoki-Veress fits to the data obtained from the cooling runs using all models indicate significant shifts of the 𝑇 g with weakening and disabling the dihedrals: The results for 𝑘 B 𝑇 𝜌𝜅 𝑇 , found for all models are indicated by triangles. The orange crosses present the MD results from Ref. [START_REF] Tsolou | Temperature and pressure effects on local structure and chain packing in cis-1,4-polybutadiene from detailed molecular dynamics simulations[END_REF] for cis-1,4-PBD at 1 atm (systems with 128 C-atoms per chain and 32 chains), which were calculated by scanning the data from Fig. 4 and Fig. 7 of [START_REF] Tsolou | Temperature and pressure effects on local structure and chain packing in cis-1,4-polybutadiene from detailed molecular dynamics simulations[END_REF]. The bottom right inset shows the static structure factors 𝑆(𝑞), that overlap for CRC and FRC models.

leads to indistinguishable compressibilities between the simulated models. An alternative approach using the Long-Lequeux theory [START_REF] Long | Heterogeneous dynamics at the glass transition in van der waals liquids, in the bulk and in thin films[END_REF] combines the more precise data of the density and compressibility obtained from the volume fluctuations extracted from long NPT simulations. This demonstrates some influence of the dihedrals, but the differences are still extremely small (Fig. 5.3). Our results are comparable to the experimental [START_REF] Schweizer | Theory of dynamic barriers, activated hopping, and the glass transition in polymer melts[END_REF] and simulation [START_REF] Tsolou | Temperature and pressure effects on local structure and chain packing in cis-1,4-polybutadiene from detailed molecular dynamics simulations[END_REF] data on PBD from the literature.

As an aside, the gathered data on the density and compressibility allowed to verify the internal consistency of our simulations by comparing the heat capacity at constant volume 𝐶 𝑉 and constant pressure 𝐶 𝑝 (Section 3.4). Focusing on the results obtained from

𝐶 𝑝 = 𝐶 𝑉 + 𝑇 𝑉 𝛼 2 𝜅 𝑇 = 𝐶 𝑉 + 𝑁 c 𝑁 𝑇 𝛼 2 𝜌𝜅 𝑇 , (5.1)
where 𝜌 = 𝑁 c 𝑁 /𝑉 is the numeric density of the system (𝑁 c and 𝑁 are the number of chains and the number if united atoms per chain correspondingly, 𝑉 is the volume), 𝛼 = 1 𝑉 𝜕𝑉 𝜕𝑇 | 𝑝,𝑁 c 𝑁 is the thermal expansion coefficient. However, the data showed a qualitative difference with respect to typical experimental results: while the simulated 𝐶 𝑝 increases upon cooling, the experimental 𝐶 𝑝 decreases [START_REF] Grebowicz | Heat capacities of 1,4-polybutadienes[END_REF]. We showed that this disagreement can be solved by considering the quantum mechanical corrections of the vibrational contribution to the 𝐶 𝑝 , which are absent in our classical modelling approach. Our findings are similar to those published recently by Honguy et al. [START_REF] Gao | Comparing simulated specific heat of liquid polymers and oligomers to experiments[END_REF], where an agreement of the heat capacity extracted from MD simulations with experimental data was achieved after taking into account the corrections for quantum effects.

A consolidated framework for describing the structure and dynamics of polymers is provided by the Rouse model [START_REF] Doi | The theory of polymer dynamics[END_REF][START_REF] Rouse | A theory of the linear viscoelastic properties of dilute solutions of coiling polymers[END_REF]. In particular the static and dynamic correlations of the Rouse modes ⃗ 𝑋 𝑝 (𝑡) are of interest. They were investigated for our models of PBD in Section 3.5. As it is expected from our results above and from previous works [START_REF] Krushev | The role of internal rotational barriers in polymer melt chain dynamics[END_REF][START_REF] Krushev | Intramolecular caging in polybutadiene due to rotational barriers[END_REF],

the influence of the dihedrals on the static Rouse mode correlations 𝐶 𝑝𝑝 = ⟨ ⃗ 𝑋 𝑝 (0) 2 ⟩ is mild. In fact, the 𝐶 𝑝𝑝 are completely determined by the conformational properties of PBD (the parametrization has been proposed in the main text) in the studied range of temperatures carrying weak temperature dependence as well (Fig. 5.4a). Contrary to the statics, the dynamic Rouse mode correlations 𝐶 𝑝𝑝 (𝑡) = ⟨ ⃗ 𝑋 𝑝 (𝑡) ⃗ 𝑋 𝑝 (0)⟩ were found to be dramatically dependent on the dihedrals. When the dihedral barriers become comparable to the thermal energy, they impose constraints on the particle movements shifting the caging effect, and thus the glass transition, at much higher temperatures than expected from the collective packing constraints without the dihedrals. The slowing down of the dynamics when approaching the 𝑇 g from above exposes itself via the characteristic feature of the glassy dynamics: the two step decay of the correlation functions with the second step being stretched with respect to the ordinary exponential decay (Fig. 5.4b).

The onset of caging appearing at higher temperatures was confirmed by the analysis of the mean-square displacements (MSD) of united atoms obtained from the simulations using all our models (Section 3.6). Our studied region of temperature 213 ≤ 𝑇 ≤ 353 K is close to the 𝑇 g of the CRC model, so the caging plateaus are clearly pronounced in between the ballistic and subdiffusive regimes of the MSD with decreasing the temperature.

With weakening the dihedrals, the same region of temperatures is becoming farther from the corresponding values of 𝑇 g for the CRC2, CRC4 and FRC models, where the cage effect does not occur yet. Going to lower temperature revealed the caging behavior for the FRC model (Fig. 5.5a). The relaxation map of the structural relaxation times extracted from the MSD data shows a super-Arrhenius increase with decreasing temperature for all models, which is common to all glass-forming systems (Fig. 5.5b). It also illustrates the significant slowing down of the structural relaxation in PBD, that depends on the height of the dihedral barriers.

Having summarized our findings above, we now revisit some critical points of our study of the bulk PBD.

• First of all most of our analysis of the structural properties is biased by the fact that we considered the density to be indistinguishable between the models. We have neglected the minor differences in the pair correlations that arise from the dihedrals (which is reasonable for our small bulk system). It would be interesting to study a bigger system to improve the statistics on subtle packing effects and the resolution The relaxation map of PBD simulated using all our models (colored symbols). The solid colored lines are the extrapolations to lower temperatures using a Vogel-Fulcher-Tammann fit, Eq. (3.84).

of static structure factors at small scattering vectors.

• Another interesting point is that in our crossover from the realistic model to the generic models by disabling all the dihedrals, the distinction between the 𝑐𝑖𝑠 and 𝑡𝑟𝑎𝑛𝑠 monomers of PBD is lost. In Section 3.4 it was shown that dihedrals that keep such a distinction (associated with the rotation along the double-bond axis of the PBD monomer) are the "hard" variables and the dynamics is determined by the other "soft" dihedrals. A deeper insight on the dynamical properties of PBD can be obtained keeping the "hard" dihedrals unmodified and disabling only the "soft" ones2 .

• In a different direction, the speed up of the dynamics using weaker dihedrals while preserving the structure of the polymer melt can be turned into an advantage as a new equilibration method. In principle, it should be possible to use the weaker dihedral configuration to equilibrate the large-scale structure of the system faster.

Then, as a second step one can reintroduce the realistic potential to equilibrate the local structure. Using this trick, it should be possible to obtain equilibrated configurations of the CRC model in the supercooled regime around 𝑇 𝑔 which would be impossible to obtain by brute force dynamics.

• The last, but not least remark can be made when relating our results to the common theories that describe the glassy dynamics, such as the mode-coupling theory [START_REF] Krushev | The role of internal rotational barriers in polymer melt chain dynamics[END_REF][START_REF] Krushev | Intramolecular caging in polybutadiene due to rotational barriers[END_REF][START_REF] Paul | Glass transition in 1,4-polybutadiene: mode-coupling theory analysis of molecular dynamics simulations using a chemically realistic model[END_REF] and the more recent elastically collective nonlinear Langevin equation (ECNLE) theory [START_REF] Mirigian | Elastically cooperative activated barrier hopping theory of relaxation in viscous fluids: II. Thermal liquids[END_REF][START_REF] Mirigian | Dynamical theory of segmental relaxation and emergent elasticity in supercooled polymer melts[END_REF][START_REF] Mirigian | Elastically cooperative activated barrier hopping theory of relaxation in viscous fluids. i. general formulation and application to hard sphere fluids[END_REF]. Being inspired from and tested on the simple glass forming liquids, they attempt to relate the dynamics with the structure of the system. The structure of the liquids is described by the spatial correlations, which are typically treated using central symmetry and averaged isotropically. While this is completely reasonable for the liquids, where the pair correlations are the most important, it is not clear whether this averaging can be applied for the systems with three-and/or four-point potentials, such as the dihedral potentials. The qualitative description, of course, should hold, since the purely structural constraints in liquids and the energetic constraints in polymers effectively lead to the same caging behavior. Yet, the different energy scales between the dihedral barriers and pair interaction barriers result in the cage effect appearing at higher temperatures for polymers. Therefore, it appears that an important challenge is to incorporate intrachain barriers into a liquid-state-theory based approach.

Film systems

The simulations of the PBD films supported on a graphite wall showed that disabling the dihedral potential does not lead to visible changes in the density profile of the films.

The profiles overlap not only in the bulk-like middle of the film, but also in the region next to the wall and at the free surface (Fig. 5.6a). Extracting the thickness of the films using the Gibbs dividing surface (GDS) method revealed a tiny difference: the thickness of the films simulated using the CRC model are a fraction of an Ångström lower than the ones for the FRC model (Fig. 5.6b). We attribute this to the better averaging statistics in the films (they are ∼ 18 times bigger than the bulk systems), so the influence of the dihedrals on the density should become more visible (though it was not detected by a visual inspection of the density profiles).

Upon cooling the films, the temperature dependence of the thickness show a change of slope (Fig. 5.6b) which suggests the glass transition at different temperatures depending on the presence of the dihedral potentials. Extracting the values of 𝑇 g for both models 5.2. FILM SYSTEMS 133 using the Dalnoki-Veress approach revealed that they agree with the bulk values within the precision of our small systems and fast cooling rates.

In order to study the reasons for such an agreement of the 𝑇 g the dynamical properties across the film were studied in a layer-resolved fashion. Two characteristic quantities were chosen: the mean-square displacements (MSD), that describe the dynamics of the united atoms, and the autocorrelation of the dipole moment vector (DACF) associated to each monomer (with the perpendicular to the film component being experimentally relevant).

A consistent picture in the time window of the structural relaxation was observed between both quantities in case of the CRC model: the dynamics is slowed down at the graphite wall (which agrees with the previously reported results of Refs. [START_REF] Solar | The dielectric axation in polymer films: a comparison between experiments and atomistic simulations[END_REF][START_REF] Solar | Dielectric -relaxation of 1,4-polybutadiene confined between graphite walls[END_REF]) and sped up at the free surface in both parallel (𝑥𝑦) and perpendicular (𝑧) to the film directions. While the same behavior is observed in case if the FRC model for the 𝑥𝑦 direction, the structural relaxation times in 𝑧 direction extracted from the MSD correspond to a different regime in the DACF, which is slower at both interfaces (the hints of a similar slowing down was also observed for the layers at the free surface for the CRC model, but at longer times). The early appearance of such a relaxation step was attributed to the fact, that in the FRC model there are no dihedral potentials to preserve the configuration of the monomer. So the dipole moment vector autocorrelation (calculated in the same way as for the CRC The temperature dependence of the film thickness obtained by the GDS method from the simulations using the CRC (red) and the FRC (blue) models. The solid black lines are the fits using Eq. (4.4) with the parameters from Table 4.1. model) is sensitive not only to the reorientation of the monomer (as it is for the CRC model), but also to the rearrangements of the united atoms within the monomer.

Further analysis of the structural relaxation times extracted from the MSD for both models is summarized in Fig. 5.7. It appears, that the slowing down effect of the wall and the speeding up at the free surface average out, which explains the agreement of the values of the total 𝑇 g with the the bulk. Moreover, the surface effects are mildly affected by the presence of the dihedrals themselves. They act with respect to the total order of magnitude of the dynamics, which is set by the dihedrals across the whole system. The quantitative analysis of the relaxation time gradients by fitting the data using the widely assumed empirical double-exponential form [START_REF] Schweizer | Progress towards a phenomenological picture and theoretical understanding of glassy dynamics and vitrification near interfaces and under nanoconfinement[END_REF][START_REF] Phan | Influence of longer range transfer of vapor interface modified caging constraints on the spatially heterogeneous dynamics of glass-forming liquids[END_REF]

𝜏 (𝑧, 𝑇 ) = 𝜏 b (𝑇 ) exp [-𝐴(𝑇 ) exp (- 𝑧 𝜉 )] , (5.2)
showed a similarity for the relaxation in the parallel to the film direction between the CRC and FRC models at the temperatures that have the same order of magnitude of the relaxation in the bulk-like region in the middle of the film. Being at different temperatures (densities), these results suggest that at the free surface the gradients of the dynamics in the 𝑥𝑦 direction can be related to the fact that at the free surface the cage surrounding each particle is open, which facilitates the relaxation. For the relaxation in the perpendicular to the film direction, the gradients are different between the models and temperatures suggesting a complex interplay between the confinement effects, the density and the dihedrals.

After having summarized our results for the films, let us also review the critical points in them.

• In general, a significant amount of our data suffers from insufficient statistics and lack of time for more detailed analysis. Moreover, the slowing down of the dynamics with decreasing the temperature barely allows to reach the structural relaxation time within the time of the simulations. This is the most crucial at the wall, where the dynamics is even slower. Therefore, more simulation time would be required to improve our analysis.

• As a second point, we would like to come back to the discussion of the overlapping density profiles. In principle, it is possible to express the density profile using pair distribution functions similarly to Eq. (3.9) 3 , by introducing the external potentials, that act on the system mimicking the surface effects. Such a description would, however, require knowledge of the radial distribution function in confinement. It would be interesting to study its behavior at the interfaces and how it is influenced by the dihedrals.

• Another remark concerns the dielectric relaxation for the FRC model. We have already noted, that disabling the dihedrals leads to loosing the distinction between the 𝑐𝑖𝑠 and 𝑡𝑟𝑎𝑛𝑠 monomers of PBD. It would be interesting to study the dipole autocorrelation in the model where these "hard" dihedrals are preserved and only the "soft" ones are disabled. Alternatively, it is common in the simulations using the generic models to consider the dipole moment to be perpendicular to the bonds in the chain backbone. Applying this approach for the double bonds of PBD may help to understand better the segmental dynamics and the influence of the dihedrals.

• It is interesting to speculate what happens if the simulations could be pushed to temperatures near the experimental 𝑇 g . While the range of the surface-induced changes of the dynamics is weak in the temperature interval of the current simulations, covering only a few nanometers, it might be expected to increase upon cooling. This low-𝑇, strongly supercooled regime is the focus of the ECNLE theory which then predicts much more long-ranged gradients. If this was true, then in thin supported films the gradients from the substrate and the free interface should interfere, erasing the bulk-like behavior in the center of the film and leading, perhaps, to a continuous decrease of the local relaxation time from high values at the substrate to low values at the free interface. It would be interesting to check whether such a highly inhomogeneous situation could be observed in simulations, e.g. by studying thinner films or by harnessing a more efficient equilibration method based on reduced dihedral barriers.

• A complex interplay between the temperature dependence of the film thickness, the depth of the perturbation of the dynamics by the surface effects, as well as their direction and amplitude determines the film-averaged glass transition temperature.

Thus, a comparison between the realistic and the coarse-grained approaches to the dynamics of polymers in confinement should be reviewed with a well defined mapping between the length and time scales, and the temperature.

Finally, we would like to address the question raised in the title of the thesis: Is the interfacial mobility in glass forming polymer film determined by collective motion or intramolecular energetics? Such a general question is very difficult to answer: a vast 5.2. FILM SYSTEMS 137 amount of experimental and simulation data (including ours) suggest that both factors are intimately intertwined. Moreover, the terms "intramolecular energetics" and "collective motion" themselves cover a wide range of different effects, that vary even depending on the studied polymer. Thus, the better question to ask would be how their interplay affects the dynamics? In our work we have focused on one of the intramolecular contributionsthe energetic constraints on the torsional (dihedral) angle -for a polybutadiene polymer.

We showed that their presence affects the collective caging effect shifting it to higher temperature, which leads to strongly non-Arrhenius glassy dynamics. While this coupling determines the structural relaxation in the system, it has an extremely small effect on the density (for PBD!). At the interfaces of the films, there is an additional factor, that comes into play -the surface effects that are determined by the type and the properties of the confinement. They are intertwined with the energetics and the collective motion influencing each other and, thus, the interfacial mobility in the films.

Systèmes en phase volumique

Les simulations du système volumique ont révélé que la réduction ou désactivation des potentiels dièdres entraîne peu de changements dans les propriétés conformationnelles du PBD (discuté en Section 3.1.2). L'analyse du rapport caractéristique 𝐶 ∞ (Fig. à l'aide de fonctions de distribution par paire de manière similaire à l'Éq. (3.9), 6 en introduisant des potentiels externes, qui agissent sur le système en imitant les effets de surface. Une telle description nécessiterait cependant la connaissance de la 6 En fait, l'Eq. (8.1.34) de la Réf. [START_REF] Winkelmann | Statistical mechanics of phases, interfaces, and thin films[END_REF] propose une telle expression.

Fedir DEMYDIUK

Mobilité interfaciale dans les films polymères vitrifiables : est-elle déterminée par le mouvement collectif ou l'énergie intramoleculaire ? Résumé L'influence des barrières énergétiques de torsion intramoléculaires sur la structure et la dynamique du copolymère aléatoire cis-trans du 1,4-polybutadiène (PBD) a été étudiée au moyen de simulations classiques de dynamique moléculaire. Le PBD en phase volumique et un film de ~10 nm d'épaisseur supporté par une paroi de graphite ont été simulés à l'aide d'un modèle chimiquement réaliste d'atomes unifiés, qui comporte des potentiels de torsion (dièdre). Les simulations ont également été effectuées avec des potentiels dièdres réduits, ainsi que sans potentiel dièdre. Il est démontré que la présence des dièdres renforce l'effet de cage dans la masse, ce qui ralentit la dynamique de plusieurs ordres de grandeur, déplaçant ainsi la température de transition vitreuse vers des valeurs plus élevées par rapport au modèle sans dièdres. En revanche, la densité et la structure du PBD en phase volumique ne sont pratiquement pas affectées par les dièdres. De la même manière, les dièdres déterminent l'ordre de grandeur de la dynamique dans les films. Leur désactivation n'entraîne pas de changements notables du profil de densité des films, non seulement dans la région centrale, semblable au volume, mais aussi aux interfaces. La dynamique a été étudiée par une analyse à résolution de couche des déplacements carrés moyens (MSD) et de l'autocorrélation du moment dipolaire des monomères (relaxation diélectrique). L'extraction des temps de relaxation structurelle a montré que le confinement imposé par la paroi de graphite et la surface libre est le facteur dominant dans la modification des gradients de mobilité interfaciale par rapport à l'ordre de grandeur global de la dynamique fixé par les barrières de torsion à travers tout le film PBD. Mots clés : films polymères supportés, polybutadiène, simulations de dynamique moléculaire, transition vitreuse, relaxation diélectrique, gradients de relaxation, barrières intramoléculaires.

Résumé en anglais

The influence of the intramolecular torsional energetic barriers on the structure and dynamics of 1,4-polybutadiene (PBD) cis-trans random copolymer was studied by means of classical molecular dynamics simulations. Bulk PBD and a ~10 nm thick film supported on a graphite wall were simulated using a chemically realistic united-atom model, which has torsional (dihedral) potentials associated to each bond in the polymer chain. The simulations were also performed with reduced, as well as without dihedral potentials. It is demonstrated that the presence of the dihedrals enhances the cage effect in the bulk, which slows down the dynamics by several orders of magnitude, thus shifting the glass transition temperature to higher values with respect to the model without the dihedrals. Meanwhile the density and the structure of bulk PBD are essentially unaffected by dihedrals. In the same way, the dihedrals determine the order of magnitude of the dynamics in the films. Disabling them does not lead to noticeable changes of the density profile of the films not only in the middle, bulk-like region, but also at the interfaces. The dynamics was studied by layer-resolved analysis of the mean-square displacements (MSD) of united atoms and the monomer dipole moment autocorrelation (dielectric relaxation). Extracting the structural relaxation times showed that the confinement imposed by the graphite wall and the free surface are the dominating factor in modifying the interfacial mobility gradients with respect to the overall order of magnitude of the dynamics set by the torsional barriers across the whole PBD film. Keywords: supported polymer films, polybutadiene, molecular dynamics simulations, glass transition, dielectric relaxation, relaxation gradients, intramolecular barriers.

2. 2 . 15 2. 2

 2152 THE UNITED-ATOM MODEL OF POLYBUTADIENE The united-atom model of polybutadiene A polybutadiene polymer chain has 2 types of carbon-hydrogen groups (Fig. 2.1a): CH and CH 2 , that are represented by united atoms (UAs) with the masses 𝑚 CH 2 = 14.027 u and 𝑚 CH = 13.019 u. A monomer consists of 4 UAs that form either a 𝑐𝑖𝑠 or a 𝑡𝑟𝑎𝑛𝑠 configuration (Fig. 2.1b). A polymer chain contains 𝑁 m = 29 monomers or 𝑁 = 4𝑁 m =

3

  types of bonds in the system: CH-CH, CH 2 -CH and CH 2 -CH 2 . They are plotted in (a) (b)

Figure 2 . 1 :

 21 Figure 2.1: (a) Chemical notation of the 1,4-polybutadiene. (b) Schematic illustration of a 1,4-polybutadiene chain, which contains 116 CH 2 and CH united atoms combined into 29 monomers.

2 CH 2

 22 -CH 2

Figure 2 . 2 :

 22 Figure 2.2: The bond stretching potential Eq. (2.1).

CHAPTER 2 .Figure 2 . 4 :

 224 Figure 2.4: (a) Schematic illustration of the torsion (or dihedral) angles in a fragment of a PBD chain. (b-d) Torsional potential 𝑈 tor (𝜙) (in kcal/mol) for the different torsion angles 𝜙 illustrated in panel (a).

Figure 2 . 5 :

 25 Figure 2.5: The Lennard-Jones potential Eq. (2.5).

  provide a detailed map of the conducted simulations in Figs. A.1 and A.2 in Appendix A. The simulations were conducted on the High Performance Computing (HPC) Center of the University of Strasbourg [59] and the Jean Zay (HPE SGI 8600) supercomputer at the Institute for Development and Resources in Intensive Scientific Computing (IDRIS, Orsay) [60].

Fig. A. 1

 1 Fig. A.1 presents a detailed map of the simulations using the original model. Aiming

Figure 2 . 6 :

 26 Figure 2.6: The bulk system in a periodic simulation box.

Figure 2 . 7 :

 27 Figure2.7: Sketch of the bulk simulation protocol consisting of a cooling run according to Eq. (2.9) followed by the equilibration and production runs. Starting from the melt configuration stored during the cooling run at a given 𝑇 (here 𝑇 = 273 K), equilibration comprises three steps: 1) NPT relaxation over 100 ns with determination of the equilibrium volume 𝑉 (𝑇 ) as the time average over the final part of the time series of the volume, 2) the instantaneous volume of the final configuration of the NPT run is rapidly deformed over 10 ps toward 𝑉 (𝑇 ), 3) NVT relaxation at 𝑉 (𝑇 ) over another 100 ns. Equilibration is followed by a production run which lasts at least 1 𝜇s depending on temperature.

Figure 2 . 8 :

 28 Figure 2.8: The film system in a simulation box.

  have also inherited the values for the timestep 𝛿𝑡 = 1 fs and the thermostat damping time Tdamp = 1000 fs. For the damping time of the barostat, we chose a ten times larger value Pdamp = 10000 fs as it is suggested by the LAMMPS documentation [43]. The default (in LAMMPS) number of chained thermostats and barostats was used: Tchain = 3, Pchain = 3.

Fig. 2 .KFigure 2 . 10 :

 2210 Fig. 2.10 presents the bond length distribution of all bond types obtained in our LAMMPS simulations at 𝑇 = 353 K using two values of the force constant in the harmonic bond potential: 2500 kcal/(molÅ 2 ) (our final value) and 1000 kcal/(molÅ 2 ) (in between

Fig. 2 .

 2 Fig. 2.11 presents bending angle distributions measured in our LAMMPS and GRO-MACS simulations along with the Boltzmann prediction 𝑝 bend (𝜃) = exp(-𝑈 bend (𝜃)/𝑘 B 𝑇 ) ∫ 𝜋 0 d𝜃 exp(-𝑈 bend (𝜃)/𝑘 B 𝑇 ) .(2.13)

Fig. 2 .

 2 Fig. 2.12 presents the complete torsion angle distributions measured in our LAMMPS and constrained GROMACS simulations along with the Boltzmann prediction 𝑝 tor (𝜙) = exp(-𝑈 tor (𝜙)/𝑘 B 𝑇 ) ∫ 2𝜋 0 d𝜙 exp(-𝑈 tor (𝜙)/𝑘 B 𝑇 ) .(2.15)

Figure 2 . 13 :

 213 Figure 2.13: Comparison of the melt (a) and the single chain (b) structure factors obtained form our LAMMPS simulations (red open squares) and constrained GROMACS simulation (blue filled circles) at 𝑇 = 273 K. The black line is the data scanned from Ref. [41].

Figure 2 . 14 :

 214 Figure 2.14: Comparison of the mean-square displacement of united atoms of PBD obtained form our LAMMPS simulations (red open squares) and constrained GROMACS simulation (blue filled circles) at 𝑇 = 273 K and 𝑇 = 353 K. The black line is the data scanned from Ref. [41].

Figure 2 .

 2 Figure 2.15: (a) Comparison of the PBD density 𝜌(𝑇 ) obtained from our LAMMPS simulations upon cooling the system in NPT conditions using a truncated LJ potential with tail corrections with the data from literature. Five independent cooling runs were averaged to obtain the density of PBD during our simulations (pink circles). (b) The dependence of PBD density 𝜌(𝑇 ) at 𝑇 = 353 K on the cutoff distance 𝑟 c and tail corrections.

3 )

 3 original model, Γ = 3.33 K/ns original model, Γ = 0.83 K/ns new model (CRC), Γ = 3.33 K/ns new model (CRC), Γ = 0.83 K/ns original model without dihedrals, Γ = 3.33 K/ns

Figure 2 .

 2 Figure2.16: The density of PBD obtained from the cooling runs using the original and the modified pair interactions. The purple line represents the cooling at the rate Γ = 3.33 K/ns using the original model. This is the same data as in the Fig.2.15. The green line is the data obtained from a slower cooling rate of Γ = 0.83 K/ns. The black and the red lines are from the simulations using the modified model. The data described above were obtained by averaging over 5 independent cooling runs (the simulation details are given in Appendix A). The blue line corresponds to a single cooling run using the original model without the dihedrals.

Figure 2 . 17 :

 217 Figure 2.17: The details of the linear (a) and logarithmic (b) spacing of the trajectory files written during our simulations using LAMMPS.

2 . 3 . 1 ) 42 CHAPTER 3 Figure 3 . 1 :

 23142331 Fig. 3.1a, the distribution of the double bond trans dihedral angle from the MD for the

  3.1b to 3.1d. As in the case of the double bonds, we see that the probability of adopting torsional angles is enhanced near 180 ∘ and supressed near 𝜙 = 0 ∘ or 360 ∘ for the FRC model. Deviations between the MD results for the dihedral distributions and the predictions from Eq. (3.1) are also visible for the CRC, CRC2 and CRC4 models. While Figs. 3.1b and 3.1c show that the deviations are small for the 𝛽 and 𝛼 trans bonds, they are large for the 𝛼 cis bond (Fig.3.1d), in particular for torsional states near 𝜙 = 0 ∘ or 360 ∘ which are suppressed in the MD simulation, contrary to the prediction of Eq. (3.1). Fig.3.2 shows a configuration of the 𝑐𝑖𝑠 monomer and the

Figure 3 . 2 :

 32 Figure 3.2: Fragment of PBD chain comprising 6 UAs. 4 inner UAs make a monomer of PBD. The bond lengths and bending angles are fixed to their equilibrium values according to Tables 2.1 and 2.2. All dihedral angles were fixed to 𝜙 = 0 ∘ representing the ideal 𝑐𝑖𝑠 configuration.

3. 1 .

 1 CONFORMATIONAL PROPERTIES 213 K ≤ 𝑇 ≤ 353 K, Fig. 3.3 shows the results of the CRC model at 𝑇 = 353 K. If bond correlations were only caused by the bond angle, one would expect from the freely rotating chain model theory an exponential decay 𝑃 1 (𝑠) = |⟨cos 𝜃⟩| 𝑠 [2]. Clearly, this prediction (dashed line in Fig. 3.3) must underestimate 𝑃 1 (𝑠) for PBD. It only accounts for correlations between consecutive pairs of bonds [𝑃 1 (𝑠 = 1)], but not for correlations between more distant pairs along the chain backbone [𝑃 1 (𝑠 > 1)], resulting from dihedral and interchain pair interactions. To fit these more distant correlations we use an exponential function 𝑃 1 (𝑠) = 𝐴 exp(-𝑠/𝑠) with 𝐴 and 𝑠 as adjustable parameters. This gives a good description of the MD data for 2 ≤ 𝑠 ≲ 25 (full line in Fig. 3.3). Fig. 3.3 shows that 𝑃 1 (𝑠 ≈ 25) ∼ 10 -3 , which is the noise level of our MD results. Since the statistical accuracy does not allow one to explore bond correlations for large 𝑠 → 𝑁 -1, the best adapted parameterization of the present data is given by

Figure 3 . 3 :

 33 Figure 3.3: Log-log plot of 𝑃 1 (𝑠) versus the number of bonds 𝑠 (1 ≤ 𝑠 ≤ 𝑁 -1) for the CRC model (red circles) at 𝑇 = 353 K. The dashed line indicates the exponential decay based on the freely rotating chain model theory, |⟨cos 𝜃⟩| 𝑠 [2]. The solid line represents the fit result to 𝐴 exp(-𝑠/𝑠) for 2 ≤ 𝑠 ≤ 20, yielding 𝐴 = 0.7507 and 𝑠 = 3.625. 𝑠 is indicated by a vertical dotted line. For 𝑠 > 25 the statistical accuracy of the data deteriorates: 𝑃 1 (𝑠) can become negative, explaining the gap in the log-log plot for 30 ≲ 𝑠 ≲ 100. The associated 'noise' does not allow an analysis of 𝑃 1 (𝑠) for 𝑠 ≪ 𝑠 ≤ 𝑁 -1.

Figure 3 . 4 :

 34 Figure 3.4: (a) Plot of the mean-square internal distance 𝑅 2 e (𝑠) divided by the number of bonds 𝑠 versus 𝑠 (1 ≤ 𝑠 ≤ 𝑁 -1) for the CRC (red circles), CRC2 (green triangles), CRC4 (purple pluses) and FRC models (blue crosses). All data refer to 𝑇 = 353 K. The (black) dashed line shows a fit to Eq. (3.5) for the FRC model (not shown for other models). The horizontal dotted lines show respectively the mean-square bond length 𝑙 2 = 2.158 Å 2 and for the FRC model the square effective bond length 𝑏 2 e = 11.58 Å 2

4 . 3

 43 in case of the original PBD model. Our CRC model essentially reproduces the same density, which fits within the experimental and simulation data from the literature (see Figs. 2.15 and 2.16 and Appendix C for an overview). In this section we complement the discussion by comparing the density obtained from cooling the PBD down to 3 K using all our models, comment on the equilibrium values at the working temperatures and analyze the glass transition.

Fig. 3 .FRCFigure 3 . 5 :

 335 Fig. 3.5 presents the mass densities 𝜌 m (𝑇 ) = 𝑚/𝑉 (𝑇 ) (𝑚 is the mass of PBD in the simulation box of volume 𝑉) of PBD obtained in our work. Different symbol shapes denote the runs, from which the data were extracted, and the colors mark the models.The empty circles represent the data obtained from 5 independent cooling runs at the

Figure 3 . 6 :

 36 Figure 3.6: (a) Contributions to the pressure arising from the different potential terms calculated in LAMMPS during the simulations at 𝑇 = 353 K. The data were averaged over 10 ns of a NVT run. (b-d) The total (a), intermolecular (c) and intramolecular (d) radial distribution functions obtained from the NVT simulations of 1 𝜇s in length at 𝑇 = 353 K. The colours denote different models: CRC -red, CRC2 -green, CRC4 -purple, FRCblue.

  Fig. 3.5 bends over about 𝑤 ≈ 30-50 K, while the typical experimental values are 𝑤 ≈ 2 K [100]. The situation is difficult for the FRC model, where it is not clear if the glassy linear regime is pronounced. While accounting for the width of the transition explicitly,

54 CHAPTER 3 Figure 3 . 7 :

 54337 Figure 3.7: The values of the fitting parameters of Eq. (3.10) depending on the fititting interval Δ𝑇. The colours denote different models: CRC -red, CRC2 -green, CRC4purple, FRC -blue. The dashed lines represent the final values of the parameters.

  consider a binary mixture consisting of 𝑁 particles with 𝑁 A particles of type A and 𝑁 B (= 𝑁 -𝑁 A ) particles of type B. The composition of the mixture is characterized by 𝑥 A = 𝑁 A /𝑁, the number fraction of A particles, and 𝑥 B = 𝑁 B /𝑁 = 1 -𝑥 A , the number fraction of B particles.

  Fig. 3.8 depicts 𝑆(𝑞) and its decomposition 𝑆 AA (𝑞) + 𝑆 BB (𝑞) and 2𝑆 AB (𝑞) according to Eq. (3.24). Panel (a) shows the results for the CRC model at 𝑇 = 353 K, while panel (b) compares the CRC and FRC models. Additionally, panel (b) plots the data versus 𝑞𝜎 AA , i.e. in LJ units for comparison with simple liquids. From panel (b) we see that the

Figure 3 . 8 :

 38 Figure 3.8: Panel (a): Collective static structure factor at 𝑇 = 353 K for the CRC model and its decomposition into partial structure factors according to Eq. (3.24): 𝑆 AA (𝑞)+𝑆 BB (𝑞) and 2𝑆 AB (𝑞). Panel (b): Same as above but for both the CRC model (lines) and the FRC model (circles). In addition, the 𝑞 axis is shown as 𝑞𝜎 AA to facilitate comparison with models where results are given in LJ units.

  Figure 3.9: Bhatia-Thornton structure factors at 𝑇 = 353 K for the CRC model (symbols with dotted lines) and the FRC model (filled circles): 𝑆 nn (𝑞) = 𝑆(𝑞), 𝑆 cc (𝑞) and 𝑆 nc (𝑞). The horizontal dashed lines represent the large-𝑞 limits: 𝑆 nn (𝑞 → ∞) = 1, 𝑆 cc (𝑞 → ∞) = 𝑥 A 𝑥 B and 𝑆 nc (𝑞 → ∞) = 0. The horizontal full line indicates 𝑘 B 𝑇 𝜌𝜅 𝑇 obtained from Fig. 3.11 by an extrapolation to 𝑞 → 0.

Figure 3 . 10 :

 310 Figure 3.10: Partial molar volumes (times density) at 𝑇 = 353 K for the CRC model (open symbols) and the FRC model (filled circles). The wave vector dependence of 𝜌𝑣 A (𝑞) and 𝜌𝑣 B (𝑞) is defined by inserting the 𝑞 dependent partial structure factors into Eq. (3.32) and Eq. (3.33). The dashed lines indicate the fits using Eq. (3.35) with estimated 𝑞 → 0 values, leading to the volumetric factor 𝛿, defined in Eq. (3.30), for the CRC and FRC models.

Figure 3 . 11 :

 311 Figure 3.11: Panel (a): Zoom of the small 𝑞 region for 𝑆(𝑞) (circles) and 𝑘 B 𝑇 𝜌𝜅 𝑇 (𝑞) (triangles) at 𝑇 = 353 K for the CRC model (red) and FRC model (blue). The wave vector dependence of 𝜅 𝑇 (𝑞) is defined by using the 𝑞 dependent partial structure factors in Eq. (3.34). The dashed line shows a polynomial fit, denoted by [𝑘 B 𝑇 𝜌𝜅 𝑇 (𝑞)] fit , in order to estimate the compressibility 𝜅 𝑇 by extrapolation to 𝑞 → 0. The resulting value 𝑘 B 𝑇 𝜌𝜅 𝑇 = 0.188 is indicated as a horizontal dotted line. Inserting [𝑘 B 𝑇 𝜌𝜅 𝑇 (𝑞)] fit , 𝛿 CRC = 1.25 from Fig. 3.10 and the fit result 𝑆 fit cc (𝑞) from panel (b) into Eq. (3.27) gives the solid (black) line. Panel (b): Zoom of the small 𝑞 region for 𝑆 cc (circles) and 𝑆 nc (triangles) at 𝑇 = 353 K for the CRC model (red) and FRC model (blue). The solid (black) line shows a polynomial fit to 𝑆 cc . The dashed lines show Eq. (3.29) when using 𝑆 fit cc and the volumetric factors 𝛿 CRC and 𝛿 FRC obtained from Fig. 3.10.

Fig. 3 .

 3 Fig.3.12 plots the compressibility as 1/√𝑘 B 𝑇 𝜌𝜅 𝑇 versus 1/𝑇. This plot format for the data was suggested by Schweizer and Saltzmann based on the Polymer-Reference-Interaction-Site theory, yielding[START_REF] Schweizer | Theory of dynamic barriers, activated hopping, and the glass transition in polymer melts[END_REF]

Figure 3 . 12 :

 312 Figure 3.12: Temperature dependence of the compressibility. The data are plotted as suggested by Eq. (3.39). The orange squares are experimental results for PBD scanned from Fig. 4 of Ref. [118]. The dashed line shows Eq. (3.39) with parameters 𝐴 = -0.435 and 𝐵 = 987 K taken from Table I of Ref. [118]. The circles show the results for the CRC model obtained from volume fluctuations [Eq. (3.38)]. As the experimental results, they are linearized and can be fit to Eq. (3.39) (not shown). The solid line indicates the fit result to the Long-Lequeux theory obtained from a joint adjustment of the CRC data for the monomer density (inset) to Eq. (3.40) and for the compressibility to Eq. (3.41).The results for 𝑘 B 𝑇 𝜌𝜅 𝑇 , found for all models are indicated by triangles. The orange crosses present the MD results from Ref.[START_REF] Tsolou | Temperature and pressure effects on local structure and chain packing in cis-1,4-polybutadiene from detailed molecular dynamics simulations[END_REF] for cis-1,4-PBD at 1 atm (systems with 128 C-atoms per chain and 32 chains), which were calculated by scanning the data from Fig.4and Fig.7of[START_REF] Tsolou | Temperature and pressure effects on local structure and chain packing in cis-1,4-polybutadiene from detailed molecular dynamics simulations[END_REF].

. 41 ) 68 CHAPTER 3 .

 41683 Here 𝜌 0 is the maximum close-packing density and 𝑇 c the temperature at which air would BULK SYSTEMS ANALYSIS become a good solvent for the polymer. Clearly, this (hypothetical) temperature should be very high so that 𝑇 g ≤ 𝑇 ≪ 𝑇 c is the temperature regime where the theory is expected to hold. The other two parameters, 𝛽 and 𝛾, are numbers introduced to account for the incompressibility of the monomers. A low-𝑇 expansion of Eqs. (3.40) and (3.41) gives back Eq. (3.39) with 𝐵 = 4𝑇 c /𝛾 and 𝐴 = (5 -𝛽)/2𝛾. Moreover, Ref. [123] showed that Eqs. (3.40) and (3.41) provide a consistent description of 𝑇 dependence of 𝜌 and 𝜅 𝑇 for supercooled polymer melts down to the glass transition. Therefore, we fit Eqs. (3.40) and (3.41) to 𝜌 and 𝜅 𝑇 of all our models.

  3.2, the simulations are carried out in two steps: First, NPT relaxation to get the average volume 𝑉 = 𝑉 (𝑇 , 𝑝, 𝑁 c 𝑁 ) corresponding to temperature 𝑇 and pressure 𝑝 (= 0). Second, NVT production runs by imposing the volume 𝑉 (𝑇 , 𝑝, 𝑁 c 𝑁 ). During these canonical production runs the internal energy 𝑈 is determined by 3.4. ENERGIES AND HEAT CAPACITY 69 𝑈 = ⟨ℋ(𝑥)⟩. However, since 𝑉 = 𝑉 (𝑇 , 𝑝, 𝑁 c 𝑁 ), we have (for simplicity the dependence on 𝑁 c 𝑁 is suppressed in the equations below) internal energy 𝑈 = 𝑈 (𝑇 , 𝑉 (𝑇 , 𝑝)) = 𝑈 (𝑇 , 𝑝) , (3.43) enthalpy 𝐻 = 𝑈 (𝑇 , 𝑝) + 𝑝𝑉 (𝑇 , 𝑝) . (3.44)

  large because 𝛽𝐾 bond ℓ 2 0 ≃ 6500 even for the smallest bond length, ℓ 0 = 1.34 Å, and the CHAPTER 3. BULK SYSTEMS ANALYSIS highest temperature studied, 𝑇 = 353 K. From Eq. (3.46) we get the average energy per bond

Fig. 3 .

 3 Fig. 3.13 shows the 𝑇 dependence of the intensive bond energy 𝑒 bond = 𝐸 bond /(𝑁 c 𝑁 ).Two data sets are compared: results from cooling runs with rate = 0.83 K/ns and results for 𝑇 ≥ 213 K after NPT and NVT equilibration. We see that both data sets agree well for 𝑇 ≥ 213 K. Moreover, the main panel of Fig.3.13 suggests that Eq.(3.49) gives an

Figure 3 . 13 :

 313 Figure 3.13: Intensive bond energy 𝑒 bond (= 𝐸 bond /(𝑁 c 𝑁 )) versus temperature 𝑇 for the CRC and FRC models. Results from continuous cooling runs with rate Γ = 0.83 K/ns (colored circles) are compared to equilibrium data obtained by NPT and NVT equilibration for 213 K ≤ 𝑇 ≤ 353 K (colored crosses). The colors mark CRC (red) and FRC (blue) models. The black dashed line depicts Eq. (3.49). The vertical dotted line indicates 𝑇 g = 192 K of the CRC model as determined by fitting Eq. (3.10) to the density. Inset: Plot of the ratio 𝑒 bond /[(1 -1/𝑁 )𝑘 B 𝑇 /2] versus 𝑇 using the same data as in the main panel.

Figure 3 . 14 :

 314 Figure 3.14: Intensive bond angle energy 𝑒 bend (= 𝐸 bend /(𝑁 c 𝑁 )) versus temperature 𝑇 for the CRC and FRC models. Results from continuous cooling runs with rate Γ = 0.83 K/ns (colored circles) are compared to equilibrium data obtained by NPT and NVT equilibration for 213 K ≤ 𝑇 ≤ 353 K (colored crosses). The colors mark CRC (red) and FRC (blue) models. The black dashed line depicts Eq. (3.54). The vertical dotted line indicates 𝑇 g = 192 K of the CRC model as determined by fitting Eq. (3.10) to the density. Inset: Plot of the ratio 𝑒 bend /[(1 -2/𝑁 )𝑘 B 𝑇 /2] versus 𝑇 using the same data as in the main panel. For 𝑇 → 0 the ratio increases to about 1.6 for the FRC model and to about 5.7 for the CRC model.

Figure 3 . 15 : 2 .

 3152 Figure 3.15: Intensive energy of the dihedral angles 𝑒 tor (= 𝐸 tor /(𝑁 c 𝑁 )) versus 𝑇 for the CRC model. Results from continuous cooling runs with rate Γ = 0.83 K/ns (red circles) are compared to equilibrium data obtained by NPT and NVT equilibration for 213 K ≤ 𝑇 ≤ 353 K (red crosses). The black dashed line depicts a fit to the equilibrium data with Eq. (3.55) (fit interval: 213 K ≤ 𝑇 ≤ 353 K), yielding 𝐴 𝜙 = -1.231379 kcal/mol, 𝐵 𝜙 = 1.393 × 10 -3 kcal/(mol K) and 𝐶 𝜙 = -4.644 × 10 -7 kcal/(mol K) 2 . In the glassy phase the black solid line indicates the fit result to 𝑒 tor (𝑇 ) = 𝐴 g 𝜙 + 𝐵 g 𝜙 𝑇 with 𝐴 g 𝜙 = -1.104315 kcal/mol, 𝐵 g 𝜙 = 6.408 × 10 -4 kcal/(mol K) (fit interval: 3.5 K ≤ 𝑇 ≤ 120 K). The vertical dotted line indicates 𝑇 g = 192 K as determined by fitting Eq. (3.10) to the density. Inset: Plot of (𝑒 tor -𝐴 𝜙 )/(𝐵 𝜙 𝑇 + 𝐶 𝜙 𝑇 2 ) versus 𝑇 to examine the quality of the fit with Eq. (3.55).

Fig. 3 .

 3 Fig. 3.15 shows that including only the first anharmonic correction (𝐶 𝜙 ) in Eq. (3.55) suffices to give an excellent fit to the simulation data in the polymer liquid (cf. inset of Fig. 3.15). Numerically, |𝐶 𝜙 𝑇 | is by about a factor of 10 smaller than 𝐵 𝜙 , implying that the anharmonic correction is relatively weak. The fit result for 𝐵 𝜙 (= 1.393 × 10 -3 kcal/(mol K)) has a reasonable value. To see this, let us calculate the energy in

Fig. 3 .

 3 Fig.3.[START_REF] Kremer | Glassy dynamics and glass transition in nanometric layers and films: a silver lining on the horizon[END_REF] shows that the description is indeed excellent. The fit to Eq. (3.56) was carried out for the equilibrated polymer liquid in the interval 213 K ≤ 𝑇 ≤ 353 K (dashed line), where the pair energies of the CRC and FRC models nearly coincide. Interestingly, the extrapolation of the fit to low 𝑇 (solid line) agrees with the FRC results from the cooling run down to the temperature regime where the glass transition of the FRC model occurs (see also the inset in Fig.3.16 supporting the quality of the fit).

Figure 3 . 16 : 5 .

 3165 Figure 3.16: Intensive LJ energy 𝑒 pair (= 𝐸 pair /(𝑁 c 𝑁 )) versus 𝑇 for the CRC and FRC models. Results from continuous cooling runs with rate Γ = 0.83 K/ns (colored circles) are compared to equilibrium data obtained by NPT and NVT equilibration for 213 K ≤ 𝑇 ≤ 353 K (colored crosses). The green dashed line depicts the result of a fit of Eq. (3.56) to the equilibrium data in the interval 213 K ≤ 𝑇 ≤ 353 K. The fit yields 𝐴 RT = -1.8340 kcal/mol and 𝐵 RT = 0.0223 (kcal/mol K) 3/5 . The thin full line shows the extrapolation of this fit to lower 𝑇. In the glassy phase of the CRC model the dash-dotted line indicates the fit result to 𝑒 pair (𝑇 ) = 𝐴 g pair + 𝐵 g pair 𝑇 with 𝐴 g pair = -1.409614 kcal/mol, 𝐵 g pair = 5.286 × 10 -4 (kcal/mol K) (fit interval: 3.5 K ≤ 𝑇 ≤ 120 K). The vertical dotted line indicates 𝑇 g = 192 K for the CRC model as determined by fitting Eq. (3.10) to the density. Inset: Plot of the ratio 𝐸 pair divided by Eq. (3.56) versus 𝑇 to examine the quality of the fit of the equilibrium data by Eq. (3.56).

( 3 .c

 3 49)] and bond angle [Eq. (3.54)] gives the following contribution to the specific heat 1 𝑁 c 𝑁 𝜕 𝜕𝑇 (𝐸 kin + 𝐸 bond + 𝐸 ang )| 𝑝,𝑁 c 𝑁 p from enthalpy fluc. c v from energy fluc.

Figure 3 . 17 :

 317 Figure 3.17: Specific heat at constant pressure 𝑐 𝑝 and constant volume 𝑐 𝑉 versus 𝑇 from equilibrium simulations of the CRC model. The symbols present the results from enthalpy fluctuations for 𝑐 𝑝 (circles) and energy fluctuations for 𝑐 𝑉 (squares). The solid line indicates 𝑐 𝑝 obtained when inserting Eqs. (3.60) to (3.62) into Eq. (3.58). The dashed line shows 𝑐 𝑉 calculated from Eq. (3.64) as discussed in the text. The horizontal dotted line indicates the sum of Eq. (3.60) and 𝐵 𝜙 from Eq. (3.61), yielding the constant 3.188 𝑘 B . When adding Eq. (3.62) to this constant one obtains the dash-dotted line. The difference between the circles and the dash-dotted line shows the contribution of the 𝑇 dependent part of the dihedrals. Inset: Comparison to experiment. The circles and solid line are the results for 𝑐 𝑝 from the main panel. The red dashed line shows Eq. (3.65) (divided by 4𝑅 with 𝑅 being the gas constant) representing the experimental data from Ref.[START_REF] Grebowicz | Heat capacities of 1,4-polybutadienes[END_REF]. The green squares depict the results from Eq. (3.67) using the experimental data (divided by 4𝑅) for the solid state from Fig.9of[START_REF] Grebowicz | Heat capacities of 1,4-polybutadienes[END_REF] as a proxy to account for contributions from quantum mechanical vibrations in the experimental system.

  . (3.60), Eq. (3.61) and Eq. (3.62) the solid line in Fig. 3.17 is obtained, which agrees with the result from the enthalpy fluctuations. Of course, this agreement is expected for systems at thermal equilibrium. Here it illustrates the internal consistency of the simulation, i.e., of the thermostating/barostating methods and the length of the equilibration and production runs.

3. 4 .

 4 ENERGIES AND HEAT CAPACITY 81 The splitting of the internal energy into its contributions according to Eq. (3.59) allows to obtain further insight into the 𝑇 dependence of 𝑐 𝑝 . With the numerical value for 𝐵 𝜙 from Fig. 3.15 the first term of Eq. (3.61) adds a constant contribution to the specific heat of about 0.7𝑘 B . Together with Eq. (3.60) this implies an overall constant contribution of about 3.188𝑘 B to the specific heat, resulting from the kinetic energy and hard variables (bond length, bond angle, and 𝐵 𝜙 ). This baseline accounts for the lion's share of the specific heat (dotted line in Fig. 3.17), while the 𝑇 dependent parts from Eq. (3.61) and Eq. (3.62) represent less than one 𝑘 B of 𝑐 𝑝 . These parts are dominated by

  .68) with ⃗ 𝑟 𝑛 (𝑡) being the position of 𝑛-th bead at time 𝑡, 𝑁 the number of beads in a chain, 𝑝 the mode index. Eq. (3.68) is a decomposition of the equations of motion of the polymer into the normal coordinates (modes), where the 𝑝-th mode describes the part of the chain that includes (𝑁 -1)/𝑝 segments. The inverse transformation reads ⃗ 𝑟 𝑛 (𝑡) = ⃗ 𝑋 0 (𝑡) + 2

FRCFigure 3 . 18 :

 318 Figure 3.18: Main panel: Log-log plot of the static Rouse modes autocorrelations 𝐶 𝑝𝑝 . The colored empty symbols denote the results obtained from our simulations at different temperatures. The filled black circles represent the data scanned from Fig. 3a in Ref. [40]. The vertical dashed lines correspond to 𝑝 that are multiples of the number of chemical monomers in our PBD chains 𝑁 m = 29. Inset: Zoom-in of the data in the main panel at 𝑝 > 27.

CHAPTER 3 .

 3 BULK SYSTEMS ANALYSIS they correspond exactly to the values of 𝑝 = 29, 58, 87. Those values are in fact multiples of the number of chemical monomers in our PBD chains 𝑁 m = 29 (each of them consists of 4 UAs resulting in 𝑁 = 4𝑁 m = 116 UAs per chain). Thus, the deviations of the values of 𝐶 𝑝𝑝 at 𝑝 = 29, 58, 87 from the overall trend can be explained by the fact that the cosine transform in Eq. (3.68) catches the periodicity in the PBD chain comprised of repeating sets of 4 UAs.The theoretical prediction for the 𝐶 𝑝𝑝 in the ideal Rouse chain reads[START_REF] Rouse | A theory of the linear viscoelastic properties of dilute solutions of coiling polymers[END_REF]

properties in Section 3 . 1 . 2 (

 312 cf. Eqs.(3.4) and (3.5)), that include the correction due to the stiff bending angle potential and an effective angle created by the dihedrals. By expressing the positions of UAs in Eq. (3.68) as the sums of bondvectors and plugging the result into Eq. (3.71) one obtains (after a long calculation that, due to its length, we do not reproduce in this thesis):

FRCFigure 3 . 19 :

 319 Figure 3.19: Main panel: Log-linear plot of the rescaled static Rouse modes autocorrelations 𝐶 𝑝𝑝 . The colored symbols denote the results obtained from our simulations at 𝑇 = 353 K. The colored lines with black dots are the results of Eq. (3.73). Inset: Zoom-in of the data from the main panel at 𝑝 𝑁 → 0. The colored horizontal dashed lines correspond to the values of 𝑏 2 e for each model.

Fig. 3 .

 3 Fig. 3.20 presents the functions 𝐶 𝑝𝑝 (𝑡) for 𝑝 = 1, 10, 20, 50, 80, 100 for our models of PBD at 𝑇 = 213 K. The first modes look like single decays with the relaxation times significantly faster when reducing the strength of the dihedrals (Figs. 3.20a and 3.20b).In case of the CRC2 and CRC models, with increasing the mode number 𝑝 the two step decays gradually develop with the intermediate plateaus becoming more pronounced when decreasing the temperature (Figs. 3.20c to 3.20f). The curves also become more stretched.In case of the CRC4 and FRC modes the relaxations are quite fast, so the exponential relaxations are merged with the initial relaxation.In order to characterize the relaxations quantitatively, at first, we extracted the relaxation times 𝜏 𝑝 of all modes as points, where 𝐶 𝑝𝑝 (𝑡 = 𝜏 𝑝 ) = 1/𝑒. The result is plotted in Fig.3.21. The axes are rescaled to test if the behavior of 𝜏 𝑝 follows the 𝑝 -2 prediction of Eq. (3.74). As it can be seen from the figure, only the first few modes, that correspond to large length scale relaxations have 𝑝 -2 scaling. For the smaller length scale modes deviations are observed, which can be explained by the fact that in our simulations the chains are more complicated (stiffness, dihedrals, excluded volume etc.) than the Rouse model considers. Another observation is, that depending on the model, the times are shifted up systematically for all modes with respect to the data for the FRC model.With lowering the temperature the shifts become dramatic and span over many orders of magnitude. For the lowest studied temperatures 𝑇 = 225, 213 K, the first modes of the CRC model have not even reached 1/𝑒 over the simulation time of several 𝜇s. Such

Figure 3 . 20 :

 320 Figure 3.20: Dynamic Rouse modes autocorrelations 𝐶 𝑝𝑝 (𝑡) for 𝑝 = 1, 10, 20, 50, 80, 100.The colored symbols represent the data obtained from the simulations at 𝑇 = 213 K. The colored lines with black dots are the fits using Eq.(3.76). The horizontal dashed lines correspond to the value 1/𝑒.

FRCFigure 3 . 21 :

 321 Figure 3.21: Relaxation times 𝜏 𝑝 of the dynamic Rouse modes correlations extracted as points, where 𝐶 𝑝𝑝 (𝑡 = 𝜏 𝑝 ) = 1/𝑒. The colors mark the simulation models and the symbols denote different temperatures.

FigFigure 3 . 22 :

 322 Fig. 3.22a).

  appears to become progressively violated with decreasing 𝑁 /𝑝. With 𝑠 = 𝑁 /𝑝, 𝑁 /𝑝 ≲ 10 corresponds to 𝑠 ≲ 16. According to Fig. 3.4, the regime 𝑠 ≲ 16 corresponds to chain segments where deviations from the Gaussian chain model, 𝑅 2 e (𝑠) = 𝑏 2 e 𝑠, become very prominent.

Fig. 3 . 4 𝑁- 1 ∑ 4 Figure 3 . 23 :

 3414323 Fig. 3.23 presents the 𝑔 0-4 (𝑡) obtained from the simulations using the CRC model at 𝑇 = 353 K. The united atoms at the ends of a PBD chain (orange diamonds in the figure)

96 CHAPTER 3

 963

Figure 3 . 24 :FRCFigure 3 . 25 :

 324325 Figure 3.24: Log-log plots of the mean-square displacement 𝑔 0 (𝑡) of united atoms for the CRC [panel (a), red], CRC2 [panel (b), green], CRC4 [panel (c), purple] and FRC [panel (d), blue] models . The colored symbols represent the data obtained from the simulations at different temperatures. The legend in the plot for the FRC model is also valid for other models. The colored solid lines with black dots are the predictions using the Rouse modes Eqs. (3.80), (3.82) and (3.83). The horizontal dashed line corresponds to the value 𝜎 2 .

Figure 3 . 26 :

 326 Figure 3.26: Comparison of the relaxation map of PBD simulated using the CRC and FRC models (filled symbols) with experimental and simulation references from the literature.The empty symbols represent the scanned data points, the solid lines are recreated VFT fits. Ref.[START_REF] Smith | Relationship between the 𝛼-and 𝛽-relaxation processes in amorphous polymers: insight from atomistic molecular dynamics simulations of 1,4 polybutadiene melts and blends[END_REF] are the data from old simulations using the original CRC model. Refs.[START_REF] Arbe | Study of the dynamic structure factor in the 𝛽 relaxation regime of polybutadiene[END_REF][START_REF] Arbe | Merging of the 𝛼 and 𝛽 relaxations in polybutadiene: a neutron spin echo and dielectric study[END_REF] are the dielectric spectroscopy data. Ref.[START_REF] Berry | The viscosity of polymers and their concentrated solutions[END_REF] are the data based on the measurements of the monomeric friction coefficient. Ref.[START_REF] Solar | Relaxation processes and glass transition of confined polymer melts: a molecular dynamics simulation of 1,4-polybutadiene between graphite walls[END_REF] are dielecric relaxation data obtained from the recent simulations using the original CRC model. Ref.[START_REF] Schmidtke | Temperature dependence of the segmental relaxation time of polymers revisited[END_REF] are the combined data from dielectric spectroscopy (low 𝑇) and NMR (high 𝑇) for PBD of 𝑀 = 1450 g/mol and 𝑀 = 2020 g/mol (in our simulations 𝑀 = 1569 g/mol).

Figure 3 . 27 :

 327 Figure 3.27: The relaxation times obtained in our work (filled symbols with solid lines) compared to the predictions of the ECNLE theory (empty symbols) and the experimental reference (crosses).

( 4 . 2 )

 42 means that starting from a middle point 𝑧 c the position of the GDS is propagated to the left (or to the right) proportionally to the fraction of the density profile at position 𝑧 to the mean value. The positions 𝑧 G + and 𝑧 G -are marked in the Fig. 4.1 as vertical dotted lines. The thickness of the film is then obtained as ℎ = 𝑧 G + -𝑧 G -. (4.3) Since the density of PBD increases upon cooling, the thickness decreases, as it is observed from lower positions of 𝑧 G + with decreasing the temperature in the Fig. 4.1. In 4

Figure 4 . 1 :

 41 Figure 4.1:The density profiles of the films obtained from the simulations using the CRC (red) and FRC (blue) models at different temperatures. The dashed black line is the data scanned from Ref.[START_REF] Solar | Relaxation processes and glass transition of confined polymer melts: a molecular dynamics simulation of 1,4-polybutadiene between graphite walls[END_REF], where the PBD film confined between graphite was simulated using the original CRC model at 𝑇 = 353 K. The vertical dot-dashed lines are the positions of the graphite layers. The vertical dotted lines are the positions of the Gibbs dividing surfaces for the CRC model. The horizontal dotted lines correspond to the density of the bulk systems at the same temperatures.

Figure 4 . 2 :

 42 Figure 4.2:The temperature dependence of the film thickness obtained by the GDS method from the simulations using the CRC (red) and the FRC (blue) models. The circles represent the data obtained from the cooling runs, and the crosses are the averaged values from the constant temperature simulations. The solid black lines are the fits using Eq. (4.4) with the parameters fromTable 4.1.

Figure 4 . 3 : 106 CHAPTER 4 .

 431064 Figure 4.3: The values of the fitting parameters of Eq. (4.4) depending on the fititting interval Δ𝑇. The colours denote different models: CRC -red, FRC -blue. The dashed lines represent the final values of the parameters.

Figure 4 . 4 :

 44 Figure 4.4: Schematic representation of the dipole moment of the 𝑐𝑖𝑠 group of PBD. The CH 2 united atoms are treated together. The hydrogen atoms from the CH groups are separated and placed in their equilibrium positions. The charges are given in the figure in the units of electronic charge 1 e = 1.60217733 × 10 19 C along with the equilibrium bond lengths. Such a configuration results in a small non-zero dipole moment depicted as an orange arrow. Since the dihedral angle for the 𝑐𝑖𝑠 monomer is strongly localized, the dipole is considered in plane of the monomer. Due to the stiff bonds and bending angles it can also be considered perpendicular to the double bond.

4 . 6 .

 46 The bottom of the first slab coincides with the topmost layer 4.2. MEAN-SQUARE DISPLACEMENTS AND DIELECTRIC RELAXATION 109 of the graphite. The slabs have width equal to 2 σ and overlap with the neighboring slabs over a distance σ for better statistics in the data. A united atom or a dipole belongs to a slab if its position (center of mass of a corresponding monomer in case of a dipole) is

Figure 4 . 5 :

 45 Figure 4.5: Panels (a-b):The mean-square displacements 𝑔 0 (𝑡) obtained from the bulk simulations (black solid lines) and film simulations (colored solid lines) using the CRC (red) and FRC (blue) models at different temperatures averaged over all united atoms in the simulation box. Panels (c-d): The dielectric relaxation 𝐶(𝑡) obtained from the bulk simulations (black solid lines) and film simulations (colored solid lines) using the CRC (red) and FRC (blue) models at different temperatures averaged over all dipoles in the simulation box. The vertical dotted lines correspond to the structural relaxation times extracted from MSD curves. The horizontal dotted lines correspond to the values 𝐶 = 0.2 (CRC) and 𝐶 = 0.01 (FRC), obtained as averaged values of the DACF at the times where 𝑔 0 = σ 2 .

Figure 4 . 6 :

 46 Figure 4.6: Schematic illustration of the geometry of the slabs in the PBD film.

Fig. 4 .

 4 Fig. 4.9 presents a layer-resolved comparison of the projections of the dipole moment vector onto different directions. The limits of the 𝑦 axis are fixed for each model to highlight the temperature dependence of the surface effects. One immediately notices the difference in the magnitude of the dipole moment between the models. Our calculations show on average | ⃗ 𝑀 | ≈ 0.33 D (1 Debye = 0.2081943 eÅ) for the CRC model and | ⃗ 𝑀 | ≈ 0.23 D for the FRC model. The values for the CRC model also show a slight increase with

3 σ 2 3 σ 2

 3232 and 𝑔 𝑧 0 (𝑡 = 𝜏 ) = 1 respectively (the factors 2/3 and 1/3 are introduced to account for different number of degrees of freedom with respect to the bulk data). They are plotted in Fig.4.10 as a function of the position of the slab.

Figure 4 . 7 :

 47 Figure 4.7: The parallel to the film components of the MSD (a-b) and DACF (c-d) obtained from the simulations using the CRC and FRC model at 𝑇 = 240 K. The colored lines correspond to the fim slabs with blue side of the spectrum marking the slabs at the wall and the red side -at the free surface. The black long dashed line represents the bulk data and the black short dashed line is the in-plane component of the film-averaged data. In case of the MSD, the bulk curve is multiplied by 2/3 to account for different number of degrees of freedom. The horizontal dotted line represents the value2

Figure 4 . 8 :

 48 Figure 4.8: The perpendicular to the film components of the MSD (a-b) and DACF (c-d) obtained from the simulations using the CRC and FRC model at 𝑇 = 240 K. The colored lines correspond to the fim slabs with blue side of the spectrum marking the slabs at the wall and the red side -at the free surface. The black long dashed line represents the bulk data and the black short dashed line is the perpendicular component of the film-averaged data. In case of the MSD, the bulk curve is multiplied by 1/3 to account for different number of degrees of freedom. The horizontal dotted line represents the value1

Figure 4 . 9 :

 49 Figure 4.9: Projections 𝑀 𝑥 , 𝑀 𝑦 , 𝑀 𝑧 of the dipole moment obtained in each slab in the film plotted agains the slab coordinate. The colours denote different models used in the simulations: CRC -red, FRC -blue. The shapes mark different directions: 𝑥 -squares, 𝑦circles, 𝑧 -triangles. The colored solid lines are the rescaled total magnitudes | ⃗ 𝑀| √ 3 . The limits of the 𝑦 axis are fixed for each model to highlight the temperature dependence of the surface effects.

4. 2 .

 2 MEAN-SQUARE DISPLACEMENTS AND DIELECTRIC RELAXATION 115

Fig. 4 . 14 .

 414 Fig. 4.14. Comparing this ample amount of data one can extract the following features: • First of all, aligning the data so that the positions of the upper GDS coincide (Fig. 4.11) leads to an impression that the relaxation time gradients at the free surface are determined only by the time in the middle of the films and are the same for both studied models. This is seen the best comparing the CRC model at 353 K (red squares) and the FRC model at 225 K (blue triangles) in Fig. 4.11. Such an observation suggests that the scaling of the fitting parameters in Eq. (4.7) with 𝜏 b (𝑇 )

Figs. 4 .

 4 Figs.[START_REF] De Gennes | Scaling concepts in polymer physics[END_REF].14g and 4.14h show that this is not the case. While the amplitudes of the relaxation time gradients can be discussed, the values of the decay length are very scattered. Moreover, only the amplitudes of the gradients of the dynamics in 𝑧 direction can probably be averaged to the same constant value for both models at all temperatures. The data in the 𝑥𝑦 direction actually show the same scaling with log 10 (𝜏 b ) as for the free surface (empty symbols in Fig.4.14c). Another observation

Figure 4 . 10 :

 410 Figure 4.10: The natural logarithm of the structural relaxation times extracted from the MSD data at different slabs in the film. The symbols denote different temperatures and the colors mark different models.

Figure 4 . 11 :

 411 Figure 4.11: Subset of the data from Fig. 4.10 for 𝑧 > 70 Å (next to the free surface of the film). Panel (a) presents the relaxation times in 𝑥𝑦 direction; panel (b) presents the relaxation times in 𝑧 direction. The legend is the same as in Fig. 4.10. The data were shifted so that the positions of the Gibbs dividing surfaces at the free surface of the film coincide at 𝑧 = 110 Å for all models and temperatures (vertical dotted line).

Figure 4 . 12 : 120 CHAPTER 4

 4121204 Figure 4.12: Subset of the data from Fig. 4.10 for 𝑧 < 50 Å (next to the wall). Panel (a) presents the relaxation times in 𝑥𝑦 direction; panel (b) presents the relaxation times in 𝑧 direction. The legend is the same as in Fig. 4.10.

Figure 4 . 13 :

 413 Figure 4.13: Structural relaxation times extracted from MSD data in film slabs. Filled symbols correspond to the perpendicular to the film direction, empty symbols correspond to the parallel to the film component. Colors mark different models. Colored lines are the results of the fits using Eq. (4.7). The density profiles are given for comparison as black solid lines. Vertical black dotted lines are the positions of the Gibbs dividing surfaces.

Figure 4 . 14 :

 414 Figure 4.14: The values of the fitting coefficients in Eq. (4.7) for the relxation time gradients at the free surface (a-b,e-f) and the wall (c-d,g-h) of the PBD supported films simulated using the CRC (red) and FRC (blue) models. The empty symbols represent the in-plane relaxations, the same filled symbols are the relaxation in the perpendicular to the film direction.

Figure 4 .

 4 Figure 4.15:The profile of the relaxation times obtained from Ref.[START_REF] Peter | Thickness-dependent reduction of the glass-transition temperature in thin polymer films with a free surface[END_REF] where the polymer film supported on a structureless wall was simulated using a bead-spring polymer model.

  from the MSD data correspond to different regimes in the dielectric relaxation between the CRC and the FRC models. Qualitatively different behavior of the dielectric relaxation between the parallel and perpendicular to the film components was observed: the in-plane relaxation is slower at the wall and faster at the free surface, while the perpendicular to the film relaxation is slowed down at both interfaces.

124 CHAPTER 4 . 48 Figure 4 . 16 :

 124448416 Figure 4.16: The natural logarithm of the structural relaxation times extracted from the MSD data at different slabs in the film. The symbols denote different temperatures and the colors mark different models.

126 CHAPTER 5 .

 1265 CONCLUSIONS AND OUTLOOKgeometry (supported film on graphite) using the CRC and the FRC models in the temperature range 213 ≤ 𝑇 ≤ 353 K.

1 κFRCFigure 5 . 1 :FRCFigure 5 . 2 :

 15152 Figure 5.1: Temperature dependence of the characteristic ratio 𝐶 ∞ from Eq.(3.6) for the CRC (red circles), CRC2 (green triangles), CRC4 (purple pluses) and FRC models (blue crosses). The plot format, ln 𝐶 ∞ versus 𝑇, is motivated by Eq. (3.7) when assuming 𝜅 to be constant. The full line presents a linear fit to the FRC data, yielding 𝜅 = -0.52 K -1 . The dashed line shows the fit result to the CRC model, leading to 𝜅 = -0.64 K -1 . The dotted line indicates the experimental result, 𝜅 = -0.10 K -1 (obtained for 298 K ≤ 𝑇 ≲ 373 K), from Table VIII of Ref.[START_REF] Mays | Temperature dependence of unperturbed dimensions for stereoirregular 1,4-polybutadiene and poly(𝛼-methylstyrene)[END_REF].

128 CHAPTER 5 Figure 5 . 3 :

 128553 Figure 5.3: Temperature dependence of the compressibility. The data are rectified as suggested by Eq. (3.39). The orange squares are experimental results for PBD scanned from Fig. 4 of Ref. [118]. The dashed line shows Eq. (3.39) with parameters 𝐴 = -0.435 and 𝐵 = 987 K taken from TableIof Ref.[START_REF] Schweizer | Theory of dynamic barriers, activated hopping, and the glass transition in polymer melts[END_REF]. The circles show the results for the CRC model obtained from volume fluctuations [Eq.(3.38)]. As the experimental results, they are linearized and can be fit to Eq. (3.39) (not shown). The solid line indicates the fit result to the Long-Lequeux theory obtained from a joint adjustment of the CRC data for the monomer density (top left inset) to Eq.(3.40) and for the compressibility to Eq.(3.41). The results for 𝑘 B 𝑇 𝜌𝜅 𝑇 , found for all models are indicated by triangles. The orange crosses present the MD results from Ref.[START_REF] Tsolou | Temperature and pressure effects on local structure and chain packing in cis-1,4-polybutadiene from detailed molecular dynamics simulations[END_REF] for cis-1,4-PBD at 1 atm (systems with 128 C-atoms per chain and 32 chains), which were calculated by scanning the data from Fig.4and Fig.7of[START_REF] Tsolou | Temperature and pressure effects on local structure and chain packing in cis-1,4-polybutadiene from detailed molecular dynamics simulations[END_REF]. The bottom right inset shows the static structure factors 𝑆(𝑞), that overlap for CRC and FRC models.

5. 1 .

 1 BULK SYSTEMS 129 the simulations using the CRC model, it was shown that the thermodynamic relation between them holds:

Figure 5 . 4 :

 54 Figure 5.4: Panel (a): Log-log plot of the static Rouse modes autocorrelations 𝐶 𝑝𝑝 . The colored empty symbols denote the results obtained from our simulations at different temperatures. The filled black circles represent the data scanned from Fig. 3a in Ref. [40]. Panel (b): Dynamic Rouse modes autocorrelations 𝐶 𝑝𝑝 (𝑡) for 𝑝 = 20 at 𝑇 = 213 K. The colored symbols represent the data obtained from the simulations. The colored lines with black dots are th stretched exponential fits using Eq. (3.76).

Figure 5 . 5 :

 55 Figure 5.5: Panel (a): Log-log plot of the mean-square displacement 𝑔 0 (𝑡) of united atoms. The colored symbols represent the data obtained from the simulations for different models. The colored solid lines with black dots are the predictions using the Rouse modes Eqs. (3.80), (3.82) and (3.83). Panel (b):The relaxation map of PBD simulated using all our models (colored symbols). The solid colored lines are the extrapolations to lower temperatures using a Vogel-Fulcher-Tammann fit, Eq. (3.84).

Figure 5 . 6 :

 56 Figure 5.6: Panel (a): The density profiles of the films obtained from the simulations using the CRC (red) and FRC (blue) models at different temperatures. The vertical dot-dashed lines are the positions of the graphite layers. The vertical dotted lines are the positions of the Gibbs dividing surfaces for the CRC model. The horizontal dotted lines correspond to the density of the bulk systems at the same temperatures. Panel (b):The temperature dependence of the film thickness obtained by the GDS method from the simulations using the CRC (red) and the FRC (blue) models. The solid black lines are the fits using Eq. (4.4) with the parameters fromTable 4.1.

Figure 5 . 7 :

 57 Figure 5.7: The natural logarithm of the structural relaxation times extracted from the MSD data at different slabs in the film. The symbols denote different temperatures and the colors mark different models.

1 κFRCFigure 5 . 8 :FRCFigure 5 . 9 :

 15859 Figure 5.8: Dépendance en température du rapport caractéristique 𝐶 ∞ , Eq. (3.6), pour les modèles CRC (cercles rouges), CRC2 (triangles verts), CRC4 (plus violet) et FRC (croix bleues). Le format du graphique, ln 𝐶 ∞ en fonction de 𝑇, est motivé par l'Éq. (3.7)en supposant que 𝜅 est constant. La ligne pleine présente un ajustement linéaire aux données de la FRC, donnant 𝜅 = -0, 52 K -1 . La ligne en pointillé présente le résultat de l'ajustement au modèle CRC, conduisant à 𝜅 = -0, 64 K -1 . La ligne en pointillés indique le résultat expérimental, 𝜅 = -0, 10 K -1 (obtenu pour 298 K ≤ 𝑇 ≲ 373 K), d'après le Tableau 8 de la Réf.[START_REF] Mays | Temperature dependence of unperturbed dimensions for stereoirregular 1,4-polybutadiene and poly(𝛼-methylstyrene)[END_REF].

Figure 5 . 11 :Figure 5 . 12 : 5 •Figure 5 . 13 :

 5115125513 Figure 5.11: Panneau (a) : Tracé log-log des autocorrélations des modes de Rouse statiques 𝐶 𝑝𝑝 . Les symboles vides colorés indiquent les résultats obtenus à partir de nos simulations à différentes températures. Les cercles noirs remplis représentent les données scannées de la Fig. 3a dans la Réf. [40]. Panneau (b) : Autocorrélations des modes de Rouse dynamiques 𝐶 𝑝𝑝 (𝑡) pour 𝑝 = 20 à 𝑇 = 213 K. Les symboles colorés représentent les données obtenues à partir des simulations. Les lignes colorées avec des points noirs sont les ajustements exponentiels étirés en utilisant Éq. (3.76).

Figure 5 . 14 :

 514 Figure 5.14: Le logarithme naturel des temps de relaxation structurelle extraits des données MSD à différentes profondeures du film. Les symboles indiquent les différentes températures et les couleurs les différents modèles.

  fonction de distribution radiale sous confinement. Il serait intéressant d'étudier son comportement aux interfaces, et comment il est influencé par les dièdres.• Une autre remarque concerne la relaxation diélectrique pour le modèle FRC. Nous avons déjà noté que la désactivation des dièdres conduit à la perte de la distinction entre les monomères 𝑐𝑖𝑠 et 𝑡𝑟𝑎𝑛𝑠 du PBD. Il serait intéressant d'étudier l'autocorrélation des dipôles dans le modèle où ces dièdres "durs" sont préservés et où seuls les dièdres "mous" sont désactivés. Alternativement, il est courant dans les simulations utilisant les modèles génériques de considérer que le moment dipolaire est perpendiculaire aux liaisons du squelette de la chaîne. L'application de cette approche pour les doubles liaisons du PBD peut aider à mieux comprendre la dynamique segmentale et l'influence des dièdres.• Il est intéressant de spéculer sur ce qui se passerait si les simulations pouvaient être poussées à des températures plus proches des 𝑇 g expérimentale. Bien que la portée des changements de la dynamique induits par la surface soit faible dans l'intervalle de température des simulations actuelles, couvrant seulement quelques nanomètres, on pourrait s'attendre à ce qu'elle augmente lors du refroidissement. Ce régime à faible 𝑇, fortement surfondu, est au centre de la théorie ECNLE qui prédit alors des gradients beaucoup plus étendus. Si cela était vrai, alors dans les films minces supportés, les gradients du substrat et de l'interface libre devraient interférer, effaçant le comportement de masse au centre du film et conduisant, peut-être, à une diminution continue du temps de relaxation local des valeurs élevées au substrat aux valeurs basses à l'interface libre. Il serait intéressant de vérifier si une telle situation très inhomogène peut être observée dans les simulations, par exemple en étudiant des films plus fins ou en exploitant une méthode d'équilibration plus efficace basée sur des barrières dièdres réduites.• Une interaction complexe entre la dépendance de la température de l'épaisseur du film, la profondeur de la perturbation de la dynamique par les effets de surface, ainsi que leur direction et leur amplitude, détermine la température de transition vitreuse moyenne du film. Ainsi, une comparaison entre les approches réalistes et à gros grains de la dynamique des polymères en confinement devrait être entreprise à nouveau avec une correspondance bien définie entre les échelles de longueur et de temps, et de la température.Enfin, nous aimerions aborder la question soulevée dans le titre de la thèse : La mobilité interfaciale dans un film polymère vitrifiable est-elle déterminée par le mouvement collectif ou par l'énergie intramoléculaire ? Il est très difficile de répondre de manière franche à une question aussi générale : un grand nombre de données expérimentales et de simulation (y compris les nôtres) suggèrent que les deux facteurs sont intimement liés. De plus, les termes "énergie intramoléculaire" et "mouvement collectif" couvrent eux-mêmes un large éventail d'effets différents, qui varient même en fonction du polymère étudié. Par conséquent, la meilleure question à poser serait de savoir comment leur interaction affecte la dynamique. Dans notre travail, nous nous sommes concentrés sur l'une des contributions intramoléculaires -les contraintes énergétiques sur l'angle de torsion (dièdre) -pour un polymère de polybutadiène. Nous avons montré que leur présence affecte l'effet de cage collectif en le déplaçant vers une température plus élevée, ce qui conduit à une dynamique vitreuse fortement non-Arrhenius. Alors que ce couplage détermine la relaxation structurelle du système, il a un effet extrêmement faible sur la densité (pour le PBD !). Aux interfaces des films, un facteur supplémentaire entre en jeu : les effets de surface qui sont déterminés par le type et les propriétés du confinement. Ils sont intimement liés à l'énergétique et au mouvement collectif qui s'influencent mutuellement et, par conséquent, à la mobilité interfaciale dans les films.

Table 2 .

 2 1: Parameters for the bond potential of Eq. (2.1).

	Bond length type CH-CH CH 2 -CH 2 CH-CH 2	𝐾 bond (kcal/molÅ 2 ) 2500 1.34 ℓ 0 (Å) 2500 1.53 2500 1.50

Table 2 .

 2

	2: Parameters for the bending po-tential of Eq. (2.2). They are identical to the original model [26].
	Bond angle type CH 2 -CH-CH CH-CH 2 -CH 2	𝐾 bend (kcal/mol) 68.1475 125.896 𝜃 0 (deg) 66.5925 111.652

6

∑ 𝑛=1

Table 2 .

 2

	Torsion angle 𝛼 (𝑡𝑟𝑎𝑛𝑠) 𝛼 (𝑐𝑖𝑠) 𝛽 db (𝑡𝑟𝑎𝑛𝑠) db (𝑐𝑖𝑠)	𝐴 1 0.385179 3.466610 1.050490 𝐴 2 𝐴 3 1.160540 0.360163 0.680313 -3.261500 -4.942280 0.340155 24.011000 5.316470 0.100537 27.54850 -0.160949 -8.914360 𝐴 4 𝐴 5 𝐴 6 -3.86178 -0.320148 -0.720332 1.48069 -0.640294 -3.041420 5.58258 0.280129 2.000930 24.011000 -5.316470 0.100537 -27.54850 -0.160949 8.914360

3: Coefficients 𝐴 𝑛 (in kcal/mol) of the torsional potential of Eq. (2.3) in the CRC model for different dihedral angles of PBD indicated in Fig. 2.4. realistic chain (CRC) model of PBD. As it was mentioned in Section 2.1, we also study model variants with reduced torsional potential by replacing 𝑈 tor (𝜙) for all dihedral angles 2.2. THE UNITED-ATOM MODEL OF POLYBUTADIENE 19

Table 2

 2

.

[START_REF] De Gennes | Scaling concepts in polymer physics[END_REF]

show that PBD deviates from the standard Lorentz-Berthelot mixing rule

[START_REF] Hansen | Theory of simple liquids[END_REF]

.

By

contrast to the Lorentz rule, PBD is a nonadditive mixture, since 𝜎 𝑖𝑗 ≠ (𝜎 𝑖𝑖 +𝜎 𝑗𝑗 )/2, and has enhanced mixing ability relative to the Berthelot rule because 𝜖 𝑖𝑗 > √ 𝜖 𝑖𝑖 𝜖 𝑗𝑗 . In simple liquids and metallic alloys both features are favorable 2.3. THE SIMULATION METHODOLOGY 21

Table 2 . 4 :

 24 Parameters of the Lennard-Jones potential Eq. (2.5) for the different types of interaction sites: 𝑖, 𝑗 = CH 2 , CH, GRA.

	Pair type CH ↔ CH CH 2 ↔ CH 2 CH 2 ↔ CH GRA ↔ GRA	𝜖 𝑖𝑗 (kcal/mol) 0.114999 3.38542 𝜎 𝑖𝑗 (Å) 0.107639 4.00904 0.116723 3.79256 0.055640 3.40000

Table 2

 2

	.2, whereas

𝜃 (𝜃 -𝜃 0 ) 2 . (2.11) Assuming 𝜃 to be close to the equilibrium bond angle 𝜃 0 the bending stiffnesses 𝐾 bend and CHAPTER 2. SIMULATION MODEL AND COMPUTATIONAL DETAILS LAMMPS GROMACS Figure 2.11: Bending angle distributions measured during the LAMMPS (red open squares) and GROMACS (blue filled circles) simulations. The vertical lines mark the average equilibrium values of 𝜃 0 . The black solid lines are the sum of Eq. (2.13) calculated for each bending angle type and normalized by the number of corresponding angles in the system.

𝑘 𝜃 are related to one another by

4

Table 2

 2 𝑡𝑟𝑎𝑛𝑠 configurations. The deviations from the planar configurations are possible, but quickly become penalized energetically due to large values of 𝑈 tor (𝜙) relative to the thermal energy that is 𝑘 B 𝑇 ≃ 0.7 kcal/mol at 𝑇 = 353 K 5 . The potential energy of the inflection points at 𝜙 = 90 ∘ , 270 ∘ corresponds to about 12 000 K and that of the maxima at 𝜙 = 180 ∘ for cis or 𝜙 = 0 ∘ , 360 ∘ for trans to about 24 000 K. These temperatures

	.3 were obtained by fitting Eq. (2.3) to Eq. (2.14) for
	each type of dihedral angle (the details are described in Appendix B).
	The double-bond dihedral potentials are much stronger than the 𝛼 and 𝛽 dihedral
	potentials. This ensures that the double-bond torsion angles stay close to the planar average
	equilibrium state during the simulations. It also implies that a monomer initially prepared
	in a 𝑐𝑖𝑠 or 𝑡𝑟𝑎𝑛𝑠 state will remain in this state for all times of current simulations with a
	length of few microseconds. In the initial configurations, that we received from the work
	[36], the potentials were randomly assigned to the monomers in a chain with the ratio of
	45% 𝑐𝑖𝑠 and 55% are much larger than the temperatures studied (𝑇 ≤ 353 K) so that the corresponding
	angles are never adopted in the simulation 6 . We can estimate typical fluctuations around
	this state by a harmonic expansion of 𝑈 tor (𝜙) around 𝜙 min . This expansion provides
	an excellent approximation up to 𝑈 tor (𝜙 ≈ 𝜙 min ± 20 ∘ ) ≈ 2.7 kcal/mol with a force
	constant 𝑘 𝜙 = 0.0133 kcal/(mol deg 2 ), implying that √⟨(𝜙 -𝜙 min ) 2 ⟩ = √𝑘 B 𝑇 /𝑘 𝜙 ≈ 7 ∘
	for 𝑇 = 353 K. Deviations from the equilibrium state are therefore small, though larger
	than for ℓ and 𝜃. The double bonds, together with the bond lengths and bond angles, are
	therefore stiff degrees of freedom.
	The main source of chain flexibility and conformational rearrangements in PBD stems
	from torsional rotations around the 𝛽 bond and the 𝛼 cis and 𝛼 trans bonds. This is the
	case because the maximum potential energy of these dihedrals is an order of magnitude
	smaller than for the double bonds, as seen from Fig. 2.4. The energy difference between
	trans and gauche states of the 𝛽 potential is about 0.41 kcal/mol ≈ 0.58𝑘 B 𝑇 at 353 K and
	the barrier at 𝜙 ≃ 180 ∘ ± 60 ∘ is about 3.84 kcal/mol ≈ 5.48𝑘 B 𝑇 at 353 K. The barriers

Table 2 .

 2 5: Lennard-Jones energy parameters for the different types of united atoms 𝑖, 𝑗 = CH 2 , CH. Second column: 𝜖 𝑖𝑗 of the original model[START_REF] Smith | United atom force field for molecular dynamics simulations of 1,4polybutadiene based on quantum chemistry calculations on model molecules[END_REF][START_REF] Yelash | Three-step decay of time correlations at polymersolid interfaces[END_REF]. Third column: 𝜖 𝑖𝑗 for the shifted and smooth model obtained from the original model by multiplying 𝜖 𝑖𝑗 by 𝜆 = 1.149462. Fourth column: 𝜎 𝑖𝑗 , which is the same for both models. 𝑗 = CH 2 , CH), that is, we replace 𝜖 𝑖𝑗 of the original model by 𝜆𝜖 𝑖𝑗 . The resulting values of 𝜖 𝑖𝑗 are given in the third column of Table

	Pair type CH ↔ CH CH 2 ↔ CH 2 CH 2 ↔ CH	original 𝜖 𝑖𝑗 (kcal/mol) 0.1000460 0.0936436 0.1015460	𝜖 * 𝑖𝑗 = 𝜆𝜖 (kcal/mol) 0.114999 3.38542 𝜎 𝑖𝑗 (Å) 0.107639 4.00904 0.116723 3.79256
	𝜖 as for the original model, and determine again the average pair energy 𝐸 pair (𝑇 ). The
	NVT simulations were performed in LAMMPS only and consisted of an equilibration run
	of 10 ns followed by a production run of another 10 ns to measure the pair energy and
	pressure. At both 𝑇 = 273 K and 353 K the run length of 10 ns allows the UAs to displace
	on average over a distance of at least 2 particle diameters. As expected due to the shorter
	range of the modified model, we find |𝐸 pair | < |𝐸 ref pair |. Since 𝑈 pair (𝑟) is proportional to 𝜖,
	we define 𝜆 by	𝜆 =	𝐸 ref pair 𝐸 pair	.	(2.20)
	At 𝑇 = 273 K this gives 𝜆 = 1.149462. We utilize this value for all pairs 𝑖𝑗 of united
	atoms (𝑖,

Table 3 .

 3 1: Characteristic ratio 𝐶 ∞ from Eq. (3.6) for the various models studied [Section 2.2] at temperature 𝑇.

	𝑇 (K) CRC CRC2 CRC4 FRC 213 6.2 6.0 5.9 5.8 225 5.8 5.9 5.9 5.7 240 5.8 5.8 5.9 5.7 253 6.1 5.9 5.7 5.6 273 5.7 5.7 5.7 5.5 293 5.7 5.7 5.6 5.5 353 5.5 5.5 5.5 5.4

Table 3 .

 3 ⟨cos 𝜙⟩ 𝐶 ∞ ⟨cos 𝜙⟩ 𝐶 ∞ ⟨cos 𝜙⟩ 𝐶 ∞ ⟨cos 𝜙⟩ 𝐶 ∞

	CRC 𝑇 (K) 213 -0.323 5.6 -0.275 5.0 -0.253 4.8 -0.246 4.7 CRC2 CRC4 FRC 225 -0.312 5.4 -0.270 4.9 -0.249 4.7 -0.244 4.7 240 -0

2: 𝐶 ∞ at temperature 𝑇 from the hindered rotation model [Eq. (3.8)]. The values for ⟨cos 𝜙⟩ were determined from the simulation. .305 5.3 -0.264 4.9 -0.246 4.7 -0.242 4.7 253 -0.302 5.3 -0.261 4.9 -0.243 4.7 -0.240 4.6 273 -0.292 5.2 -0.255 4.8 -0.242 4.6 -0.237 4.6 293 -0.285 5.1 -0.250 4.8 -0.235 4.6 -0.234 4.6 353 -0.267 4.9 -0.238 4.6 -0.226 4.5 -0.227 4.5 3.2. DENSITY 49

Table 3 .

 3

	3: The final values of the fit parameters of Eq. (3.10)
	model 𝑇 g (K) 𝜌 g (g/cm 3 ) 𝑀 × 10 -4 (g/cm 3 𝐾) 𝐺 × 10 -4 (g/cm 3 𝐾) 𝑤 (K) CRC 192 0.967 -6.6 -2.0 50 CRC2 143 1.004 -6.9 -2.2 45 CRC4 108 1.032 -7.0 -2.5 34 FRC 43 1.089 -8.5 -4.0

Structure factors, compressibility 3.3.1 Lennard-Jones interactions: Comparison of polybutadiene with the Kob-Andersen binary mixture

	From the viewpoint of the Lennard-Jones (LJ) interactions 1,4-cis-trans polybutadiene
	(PBD) is a binary mixture. Thus, it is interesting to compare PBD to models of glass-
	forming binary mixtures, an extensively studied example being the Kob-Andersen (KA)
	mixture [58, 102-104].

  The KA mixture is produced by blending A and B particles which by themselves crystallize at low temperature. The nonadditive cross interactions in the KA system are specifically designed to suppress demixing and ensuing crystallization. Good glass-forming properties are thus expected. Nonetheless, composition fluctuations at various length scales are unavoidable in any mixed system. In the cold liquid these fluctuations can lead to regions containing just a single component. If such a region is large enough and long-lived enough, it will nucleate a crystal of this component. Exactly this ordering mechanism has recently been described for the KA mixture

	𝑓 PBD 𝜖 𝜖 𝑓 KA	≃ 6.674 ≃ 2.243 . mol kcal	⇒ 𝑓 PBD 𝜖	𝜖 AA ≃ 0.718 ,	(3.18)

is convenient to describe deviations from this rule by an equation similar to Eq. (3.14), i.e.

𝜖 𝑖𝑗 = √ 𝜖 𝑖𝑖 𝜖 𝑗𝑗 [1 + 𝑓 𝜖 |𝜖 𝑖𝑖 -𝜖 𝑗𝑗 |] .

(3.17)

As before, the modified mixing rule does not change the energy parameter of identical particles, but implies stronger attractions between unlike particles relative to the Berthelot rule, if 𝑓 𝜖 > 0. Stronger attractions implies enhanced mixing ability

[START_REF] Ingebrigtsen | Crystallization instability in glass-forming mixtures[END_REF][START_REF] Binder | Polymer + solvent systems: phase diagrams, interface free energies, and nucleation[END_REF]

. Applying Eq. (3.17) to PBD and the KA mixture we get Thus, mixing is enhanced for both models. However, as for the LJ diameters, the deviations from the standard mixing rule (Berthelot rule) is stronger for the KA system than for PBD-by a factor of about 3. Again, this difference just reflects the similarity of the A 3.3. STRUCTURE FACTORS, COMPRESSIBILITY 57 and B particles in the case of PBD.

  CHAPTER 3. BULK SYSTEMS ANALYSIS 109, 114], the structure factors (appear to) vanish for PBD. Vanishing of 𝑆 cc (𝑞 → 0) implies that there are no (thermodynamic) composition fluctuations at large scales. Such fluctuations are prevented by chain connectivity: At large scales, the character of the binary mixture disappears for PBD; PBD is a homopolymer. What is the consequence of 𝑆 cc (𝑞 → 0) = 0 in view of the general relations of Eqs. (3.27) to (3.29)? From Eq. (3.27) we conclude that 𝑆(𝑞 → 0) = 𝑘 B 𝑇 𝜌𝜅 𝑇 and from Eq. (3.29) that 𝑆 nc (𝑞 → 0) = 0, provided the volumetric factor 𝛿 [cf. Eq. (3.31)] does not diverge. Although it is physically clear that 𝛿 cannot diverge, Fig. 3.10 tests this point by parametrizing the 𝑞 dependence of 𝜌𝑣 A (𝑞) and 𝜌𝑣 B (𝑞) by a polynomial fit function:

	similar oscillations are also found for simple binary mixtures [107, 109, 114], reflecting
	the (nonuniversal) local composition-composition and number-composition correlations in
	these systems.

has two CH 2 UAs and two CH UAs per monomer) we calculate 𝑆 cc (𝑞) and 𝑆 nc (𝑞) from Eqs.

(3.25)

and

(3.26)

. Fig.

3

.9 shows these structure factors together with 𝑆(𝑞) at 𝑇 = 353 K for the CRC and FRC models. Due to the small values of 𝑆 cc (𝑞) and 𝑆 nc (𝑞) the figure zooms on ordinate values below 1, thereby truncating the maxima of 𝑆(𝑞). We find that 𝑆 cc (𝑞) is positive for all 𝑞. For 𝑞 > 𝑞 * , 𝑆 cc (𝑞) oscillates around 𝑥 A 𝑥 B (= 0.25)-the value expected for an ideal (equimolar) mixture-while 𝑆 nc (𝑞) oscillates around 0 and is roughly in antiphase to 𝑆 cc (𝑞) for large 𝑞. Both 𝑆 cc (𝑞) and 𝑆 nc (𝑞) exhibit a maximum for 𝑞 near 𝑞 * . Qualitatively However, qualitative differences between simple binary mixtures and PBD occur in the limit 𝑞 → 0. While 𝑆 cc (𝑞) and 𝑆 nc (𝑞) tend to finite values for simple liquids

[107,

  [𝑘 B 𝑇 𝜌𝜅 𝑇 (𝑞)] fit = 𝑘 B 𝑇 𝜌𝜅 𝑇 + 𝐴 𝜅 𝑞 2 + 𝐵 𝜅 𝑞 4 , (3.37) where 𝐴 cc , 𝐵 cc , 𝑘 B 𝑇 𝜌𝜅 𝑇 , 𝐴 𝜅 and 𝐵 𝜅 are adjustable parameters. The simulation results for 𝑘 B 𝑇 𝜌𝜅 𝑇 (𝑞) are determined by inserting the 𝑞 dependent partial structure factors into Eq. (3.34).Fig. 3.11 shows the results of this analysis. From panel (b) we see that 𝑆 cc (𝑞) for the CRC and FRC models agree with each other and can be well described by Eq. (3.36).Using then 𝑆 fit cc (𝑞) and the results for 𝛿 from Fig.3.10, the difference for 𝑆 cc (𝑞) between the CRC and FRC models can be understood [cf. dashed lines in Fig.3.11b]. Similarly,

	1.688
	1.625
	0.375
	0.312

Table 3 .

 3 4:The values of the volumetric factor 𝛿 obtained from the partial molar volumes. Fig.3.11a shows that the 𝑞 dependent generalization of the compressibility is well described by Eq.(3.37). Within the error bars the simulation results for 𝑘 B 𝑇 𝜌𝜅 𝑇 (𝑞) for both models are the same so that no separate fit was attempted. Using [𝑘 B 𝑇 𝜌𝜅 𝑇 (𝑞)] fit , 𝛿 CRC and 𝑆 fit cc (𝑞) in Eq. (3.27), the solid line in Fig. 3.11a reveals that a good description of the simulation data for 𝑆(𝑞) is obtained for the CRC model. Due to the extra term 𝛿 2 CRC 𝑆 fit cc (𝑞) in Eq. (3.27), the compressibility plateau is less visible in 𝑆(𝑞) than in 𝑘 B 𝑇 𝜌𝜅 𝑇 (𝑞), for the simulated range of 𝑞 values. Apparently, 𝑘 B 𝑇 𝜌𝜅 𝑇 (𝑞) together with Eq. (3.37) is a more convenient method (with respect to straightforward extrapolation 𝑆(𝑞 → 0)) to determine the (reduced) compressibility 𝑘 B 𝑇 𝜌𝜅 𝑇 from canonical simulations of binary-mixture-like polymer systems. Here this analysis gives 𝑘 B 𝑇 𝜌𝜅 𝑇 = 0.188 at 𝑇 = 353 K for the CRC model (and other models considering a single fit in Fig. 3.11a).

	𝑇 (K) 𝛿 CRC 𝛿 CRC2 𝛿 CRC4 𝛿 FRC 213 0.829 1.001 1.025 1.178 225 0.905 1.012 1.045 1.176 240 0.960 1.021 1.058 1.224 253 0.865 1.017 1.099 1.217 273 0.999 1.061 1.106 1.228 293 1.029 1.078 1.102 1.260 353 1.250 1.187 1.228 1.377

  function (cf. Eqs.(3.35) to(3.37)). This function was inspired by the Taylor expansions of the structure factors up to 4th order in 𝑞 in Eqs. (3.28) and (3.32) to(3.34). While visually the fits look good, the post-factum look on the fitting parameters revealed that the coefficients at the 𝑞 4 term are either of the same order or even larger than the ones at

the 𝑞 2 term (the values are given in Appendix F). Clearly, this suggests that the Taylor expansion does not work in this region of 𝑞. There was no clear temperature dependence of the corresponding length scales observed as well.

Table 3 .

 3 5:The values of the fit parameters of the Long-Lequeux theory for the density and compressibility.

	model	𝑇 c (K) 𝜌 0 (g/cm 3 )	𝛽	𝜌 0 /𝛾 2
	CRC	920.751	1.104 1.251 3.902896e+27
	CRC2 1054.873	1.109 1.692 2.773495e+27
	CRC4	997.521	1.108 1.493 3.190466e+27
	FRC	1087.512	1.108 1.851 2.516277e+27

Table 2

 2 .2. As for the bond length, the agreement of the bond angle distribution obtained from the simulations (see Section 2.4.2) with the ones for intependent angles (i.e. with the distribution determined only by Eq. (3.50))

Table 2 .

 2 2 one finds erf([1 ± 𝑥 0 ]𝑞) ≃ 1 and so

	𝑍 bend (𝛽) ≃ √	𝜋 𝛽𝐾 bend	.	(3.52)

  Since the experimental 𝑐 𝑝 is reported per monomer, the dashed line in the inset shows Eq. (3.65) divided by 4 𝑅 with 𝑅 being the gas constant. We see that the order of magnitude is the same for both simulation and experiment, but the 𝑇 dependence is different: The simulated 𝑐 𝑝 increases upon cooling, whereas the experimental 𝑐 𝑝 decreases.

	52.63 + 0.178 𝑇 [J/(mol K)] .	(3.65)

  𝑐 vib (𝑇 )) to replace the terms 𝑐 vib + 𝑇 𝛼 2

	specific heat then reads	𝜌𝜅 𝑇	. The resulting prediction for the experimental
		𝑐 𝑝 ≈ 𝑐 vib (𝑇 ) + [	3 5	𝐵 RT 𝑇 2/5 -2|𝐶 𝜙 |𝑇] ,	(3.67)

where the term […] stems from Eqs. (3.61) and

(3.62)

. This term replaces 𝐶 conf in Eq.

(3.66)

Table 3 .

 3

	6: The values of the fitting parameters of Eq. (3.84) and the glass transition temperatures 𝑇 VFT g obtained by extrapolating the fits to 100 sec. The values obtained from the cooling runs are also given for reference.
	model log 𝜏 ∞ (ps) 𝐵 (K) 𝑇 0 (K) 𝑇 VFT g CRC -1.2 641.14 119.2 161.3 (K) 𝑇 cool g CRC2 -0.3 316.11 116.8 139.0 CRC4 -0.3 351.70 72.6 97.1 FRC -0.1 293.62 41.0 61.8	(K) 192 143 108

Table 4 .

 4 1: The final values of the fit parameters of Eq. (4.4) model 𝑇 g (K) ℎ g (Å) 𝑀 (Å/K) 𝐺 (Å/K) 𝑤 (K)

	CRC FRC	206.9 43.0	94.0 82.8	0.0766 0.0704	0.0132 0.0295	71.1 7.7

140

  RÉSUMÉ EN FRANÇAISbilles-ressorts, largement utilisés, dans lesquels les détails chimiques réalistes (y compris les dièdres) sont confondus dans une seule interaction à "gros grains". Notre modèle FRC ressemble à un modèle billes-ressorts, à une exception près : la longueur des liaisons est beaucoup plus courte par rapport à la taille des particules. Le modèle billes-ressort le plus utilisé de type Kremer-Grest[START_REF] Kremer | Dynamics of entangled linear polymer melts: a molecular-dynamics simulation[END_REF][START_REF] Bennemann | Molecular-dynamics simulation of a glassy polymer melt: rouse model and cage effect[END_REF][START_REF] Peter | Thickness-dependent reduction of the glass-transition temperature in thin polymer films with a free surface[END_REF] a une seule échelle de longueur qui détermine à la fois la taille des billes et la longueur des liaisons entre elles.

	Nos simulations ont été effectuées en volume en utilisant les 4 modèles et dans une
	géométrie confinée (film supporté sur graphite) en utilisant les modèles CRC et FRC dans
	la gamme de températures 213 ≤ 𝑇 ≤ 353 K.
	Il est à noter que pour simuler des films supportés (incluant donc une surface libre),

il n'est pas possible d'utiliser des corrections de queue ("tail corrections") pour prendre en compte la portée infinie du potentiel de paire. Les paramètres du modèle original ont dû être adaptés pour donner la même équation d'état avec un potentiel de paire à portée finie. Le chapitre sur les résultats en phase volumique est d'une part une validation de ces paramètres de modèle modifiés, et d'autre part un nouveau regard sur les résultats PBD. En effet, la puissance informatique accrue d'aujourd'hui nous permet d'obtenir des trajectoires beaucoup plus longues.

The high value of the force constant is a compromise between the experimentally observed bond vibrations and the constrained bonds of the original model, see also Section

2.4.

This is not yet our FRC model, since the pair interactions are still described by the truncated LJ potential with the tail corrections.

To get Eq. (2.12) we insert the identity cos 𝜃 -cos 𝜃 0 = -2 sin[(𝜃 + 𝜃 0 )/2] sin[(𝜃 -𝜃 0 )/2] into Eq. (2.2) and assume deviations of 𝜃 from 𝜃 0 to be small.

We utilize the value of Boltzmann constant defined in the LAMMPS code, 𝑘 B = 0.0019872067 kcal/(mol K).

When testing the reduced torsions, we found some flipping events at 1/10 of the initial height of the dihedral barriers.

The main argument of the derivation is the fact, that the angles do not change when rescaling the volume of the system. Therefore, the derivatives of the angular potentials with respect to the volume vanish. Our simulations employed an isotropic barostat, so this statement holds.

The enthalpy 𝐻(𝑥) = ℋ(𝑥) + 𝑝𝑉 (𝑥) depends on the microstate 𝑥 of the system. The specific heat is given by 𝑐 𝑝 = (⟨𝐻 2 ⟩ -⟨𝐻⟩ 2 )/(𝑁 c 𝑁𝑘 B 𝑇 2 ). However, the NPT simulations were carried out at 𝑝 = 0 pressure so that 𝐻(𝑥) = ℋ(𝑥) and 𝑐 𝑝 can be calculated from the fluctuations of the hamiltonian ℋ only.

For the CRC(2,4) models Δ𝐶 𝑝𝑝 (𝑡) = 𝐶 𝑝𝑝 (𝑡) -𝐶 ref 𝑝𝑝 models the data of the FRC model were taken as a reference. For the FRC model, the difference between the CRC4 and FRC was used.

This discussion is also motivated by the fact, that when we tried to discard the upper bound of 𝐴 𝑝 during the fits, we obtained the curves with 𝐴 𝑝 > 1, which has no physical meaning.

The lower limit of 𝑥 = 0.5 is in fact a result of the ideal Rouse chains. However by introducing chain stiffness, the value of 𝑥 increases with the 𝑥 = 0.75 being the limit of the stiff chains[2].

The periodic boundary conditions are employed in the simulations in 𝑥, 𝑦 directions.

On the contrary to the bulk systems, it was not possible to perform multiple cooling runs for the same film systems and average over them. Only the smoothing was applied in order to obtain a data point per Kelvin.

We refer the reader to the works[35,

[START_REF] Smith | A molecular-dynamics simulation study of dielectric relaxation in a 1,4-polybutadiene melt[END_REF] for the theoretical and implementation details, since (as we will show later) our calculations are identical.

Considering the width of the slabs 2 σ , this method yields nearly identical results for the structural relaxation to checking the positions in the beginning and the end of the correlation time window, or even at each time step within the time window. For large time windows the additional checks deteriorate the averaging statistics.

σ 2 .

σ 2 . The purple lines show the results of the original work[START_REF] Solar | Dielectric -relaxation of 1,4-polybutadiene confined between graphite walls[END_REF], where wider slabs were studied.

We used the value averaged over 5 slabs in the center of each film as 𝜏 b (𝑇 ). The bulk relaxation times sometimes slightly disagree with the values in the middle of the films, which, when used in the fits, lower their quality. We believe the disagreement comes from the statistical errors and insufficient equilibration rather than having a physical origin, since the films should be thick enough to have a bulk-like region in the middle.

The origin of these artifacts is not clear. Excluding these points from the fits for all data resulted in overestimating the gradients for all models and temperatures.

It is probably redundant to reduce only the "soft" dihedrals. Our results show that when reducing all of them by factor a of 4, the double bond dihedrals preserve their property of being "hard" variables effectively making the change from 𝑐𝑖𝑠 to 𝑡𝑟𝑎𝑛𝑠 impossible.

In fact, Eq. (8.1.34) of Ref.[START_REF] Winkelmann | Statistical mechanics of phases, interfaces, and thin films[END_REF] proposes such an expression.

Notre modèle de PBD possède deux types d'atomes unifiés, ce qui permet d'appliquer la même approche.

Il est probablement redondant de réduire uniquement les dièdres "mous". Nos résultats montrent qu'en les réduisant tous par un facteur 4, les dièdres de double liaison conservent leur propriété de variables "dures" rendant effectivement impossible le passage d'un état 𝑐𝑖𝑠 à un état 𝑡𝑟𝑎𝑛𝑠.

https://docs.lammps.org/dihedral_nharmonic.html.

Acknowledgements

The authors would like to acknowledge the High Performance Computing Center of the University of Strasbourg for supporting this work by providing scientific support and access to computing resources. Part of the computing resources were funded by the Equipex Equip@Meso project (Programme Investissements d'Avenir) and the CPER Alsacalcul/Big Data. This work was granted access to the HPC resources of IDRIS under the allocation 2020-101572 made by GENCI. Preparing this thesis could not have been possible without the support of many people. First of all, I would like to express my gratitude to the thesis supervisors Dr. Hendrik Meyer and Prof. Wolfgang Paul for provided guidance and advices over this period. Infinite thanks go to the members of the TSP group at the Institut Charles Sadron. In particular Prof. Joerg Baschnagel for numerous very inspiring scientific discussions, helping with the preparation of the manuscript and significant participation in the project. Also Dr. Olivier Benzerara for providing all necessary technical and every day life support and advisory. Special thanks go to Dr. Mathieu Solar for providing the initial configuration for the simulations, helping with the GROMACS tools and useful discussions. I would also like to express gratitude to the permanent members and students in the Theoretical Polymer Physics team at the Institute of Physics in Halle for the support during the visits. Additionally, I would like to mention the support of many friends that I made during these three years at the Institut Charles Sadron and the Europeran Doctoral College, as well as my family and my life partner Anastasiia Shpiruk. I would not have made it without them.

Finally, I would like to acknowledge the IdEx Inernational Doctoral Program of the University of Strasbourg for funding the project and providing exciting networking ACKNOWLEDGEMENTS opportunities with many young researchers from all areas of science.

Résumé en français

Mobilité interfaciale dans les films polymères vitrifiables : est-elle déterminée par le mouvement collectif ou l'énergie intramoleculaire ?

Nous proposons une étude systématique de l'influence des barrières intramoléculaires et de l'interaction avec un substrat sur la dynamique vitreuse, et ceci par des simulations numériques chimiquement réalistes. Notre travail se concentre sur un type particulier d'interactions intramoléculaires: les barrières de torsion (dièdres), qui contrôlent majoritairement les conformations d'une chaîne polymère, contrairement les énergies de liasons et de valence.

Afin d'étudier l'influence de contraintes de torsion réalistes sur les propriétés des polymères vitreux, nous avons adapté un modèle atomes unifiés du copolymère aléatoire 𝑐𝑖𝑠-𝑡𝑟𝑎𝑛𝑠 du 1,4-polybutadiène (PBD) pour des simulations de dynamique moléculaire classique avec le code LAMMPS [START_REF] Allen | Computer simulation of liquids[END_REF][START_REF] Lammps | Large-scale Atomic/Molecular Massively Parallel Simulator[END_REF][START_REF] Plimpton | Fast parallel algorithms for short-range molecular dynamics[END_REF]. Ce modèle a été développé précédemment à partir de calculs numériques de chimie quantique et validé par rapport à des expériences [START_REF] Smith | Molecular dynamics of a 1,4-polybutadiene melt. comparison of experiment and simulation[END_REF][START_REF] Smith | NMR spin-lattice relaxation and conformational dynamics in a 1,4-polybutadiene melt[END_REF][START_REF] Smith | A molecular-dynamics simulation study of dielectric relaxation in a 1,4-polybutadiene melt[END_REF]. Ce travail de doctorat fait suite aux travaux récents étudiants le PBD confiné entre deux parois de graphite. Le but est d'étudier le comportement avec une surface libre, ainsi que l'influence des potentiels dièdres. Notre modèle se décline dans en quatre versions , où nous faisons varier l'échelle de tous les potentiels de torsion (dièdres) dans le système :

• Modèle de chaîne chimiquement réaliste (CRC) avec les potentiels dièdres non modifiés (qui fournissent la distinction entre les monomères 𝑐𝑖𝑠 et 𝑡𝑟𝑎𝑛𝑠).

• Modèles CRC2 et CRC4 avec les potentiels dièdres réduits d'un facteur 2 et 4 respectivement.

• Modèle FRC (Freely rotating chain) avec tous les potentiels dièdres désactivés. Ici, la distinction entre les groupes 𝑐𝑖𝑠 et 𝑡𝑟𝑎𝑛𝑠 n'a pas lieu.

Une telle interpolation permet de faire le lien avec les modèles génériques de type

Appendix A

Simulation map

See next page.

This appendix gives a "genealogical tree" to trace the history of the simulated systems,

Appendix B Torsional potential

The original model of Ref. [START_REF] Smith | United atom force field for molecular dynamics simulations of 1,4polybutadiene based on quantum chemistry calculations on model molecules[END_REF] represents the torsional potential 𝑈 tor (𝜙) by the function

where the torsion angle lies in 0 ∘ ≤ 𝜙 ≤ 360 ∘ and the energy coefficients 𝑘 𝑛 are given in

Refs. [START_REF] Smith | United atom force field for molecular dynamics simulations of 1,4polybutadiene based on quantum chemistry calculations on model molecules[END_REF][START_REF] Smith | Molecular dynamics of a 1,4-polybutadiene melt. comparison of experiment and simulation[END_REF][START_REF] Yelash | How well can coarse-grained models of real polymers describe their structure? the case of polybutadiene[END_REF]. When determining the torsion angles from the Cartesian coordinates of the united atoms one gets the cosine of 𝜙. It is therefore computationally expedient to expand Eq. (B.1) in powers of cos 𝜙. This gives [START_REF] Krushev | Computersimulationen zur dynamik und statik von polybutadienschmelzen[END_REF] 𝑈

The GROMACS code provides such a representation via the Ryckaert-Bellemans function [START_REF] Ryckaert | Molecular dynamics of liquid 𝑛-butane near its boiling point[END_REF] (which we write here in the form employed by LAMMPS 1 )

The values for 𝐴 𝑛 resulting from 𝑘 𝑛 of Ref. [START_REF] Yelash | How well can coarse-grained models of real polymers describe their structure? the case of polybutadiene[END_REF] are given in Table B.1.

Since 𝑘 6 vanishes for the double bonds in the original model [START_REF] Smith | United atom force field for molecular dynamics simulations of 1,4polybutadiene based on quantum chemistry calculations on model molecules[END_REF][START_REF] Smith | Molecular dynamics of a 1,4-polybutadiene melt. comparison of experiment and simulation[END_REF][START_REF] Yelash | How well can coarse-grained models of real polymers describe their structure? the case of polybutadiene[END_REF], the choice When comparing the double bond potentials, we noticed, that in the recent works using GROMACS [START_REF] Solar | The dielectric axation in polymer films: a comparison between experiments and atomistic simulations[END_REF][START_REF] Solar | Dielectric -relaxation of 1,4-polybutadiene confined between graphite walls[END_REF][START_REF] Solar | Relaxation processes and glass transition of confined polymer melts: a molecular dynamics simulation of 1,4-polybutadiene between graphite walls[END_REF][START_REF] Yelash | Three-step decay of time correlations at polymersolid interfaces[END_REF][START_REF] Solar | Dipolar correlations in 1,4-polybutadiene across the timescales: a numerical molecular dynamics simulation investigation[END_REF] there are different sets of parameters {𝐴 𝑛 } for the cis and the trans double bonds, which does not agree with the ∘ corresponds to about 12 000 K and is much higher than the temperatures of a few hundred Kelvin typically studied. In practice, this implies that a monomer initially prepared in a cis or trans state will remain in this state for all times of current simulations with length of a few microseconds. We believe that the authors of the works [START_REF] Solar | The dielectric axation in polymer films: a comparison between experiments and atomistic simulations[END_REF][START_REF] Solar | Dielectric -relaxation of 1,4-polybutadiene confined between graphite walls[END_REF][START_REF] Solar | Relaxation processes and glass transition of confined polymer melts: a molecular dynamics simulation of 1,4-polybutadiene between graphite walls[END_REF][START_REF] Yelash | Three-step decay of time correlations at polymersolid interfaces[END_REF][START_REF] Solar | Dipolar correlations in 1,4-polybutadiene across the timescales: a numerical molecular dynamics simulation investigation[END_REF] have preferred to distinguish between the cis and trans states by turning the potentials to having single minima at 𝜙 = 0 ∘ (or 360 ∘ ) and 𝜙 = 180 ∘ , rather then two. In this way, even if a possibility of a rare escape from the barrier occurs, the cis monomer will still remain cis, and the trans monomer will still remain trans. In practice, this was achieve

Appendix C PBD density references

See next page.

APPENDIX C. PBD DENSITY REFERENCES

Appendix D

LAMMPS input scripts

Bulk system. To run: $LAMMPS -in inputScript.lmp -var sID packageId

Appendix E

Bond correlations and internal distances

The mean-square internal end-to-end distance and the bond correlation function are defined by

where ⃗ 𝑟 𝑛 is the position vector of the 𝑛th united atom (𝑛 = 1, … , 𝑁), ⃗ 𝑙 𝑛 = ⃗ 𝑟 𝑛+1 -⃗ 𝑟 𝑛 the bond vector and ⟨…⟩ the average over all configurations and chains of the system and over all possible pairs 𝑛, 𝑚 = 𝑛 + 𝑠 of a chain. A viable parameterization of MD results for

where 𝛼 = -⟨cos 𝜃⟩ (> 0). Using

where 𝛼 = e -1/𝑠 (< 1) and 𝐶 ∞ is the characteristic ratio of the infinitely long chain,

APPENDIX E. BOND CORRELATIONS AND INTERNAL DISTANCES

The characteristic ratio determines three often discussed length scales of a polymer: the effective bond length 𝑏 e [START_REF] Doi | The theory of polymer dynamics[END_REF],

Appendix F

Partial structure factor fits