Mechanical and viscoelastic properties of glass-forming polymers in the bulk and thin films : molecular dynamics study of model systems

Ivan Kriuchevskyi

Résumé Introduction

Les fondus de polymères sont des liquides composés de chaînes macromoléculaires [START_REF] Doi | The Theory of Polymer Dynamics[END_REF][START_REF] Rubinstein | Polymer Physics[END_REF].

Lors d'un refroidissement à basse température, le fondu de polymères liquide se transforme en un solide qui peut être amorphe (verre) ou semi-cristallin [START_REF] Strobl | The Physics of Polymers: Concepts for Understanding their Structures and Behavior[END_REF]. La capacité à former des plans cristallins dépend du détail chimique de chaques polymères. Seul les chaînes avec des configurations stéréo-régulières (par exemple, orientations isotactiques ou syndiotactiques des groupements chimiques latéraux [START_REF] Rubinstein | Polymer Physics[END_REF]) ou des chaînes sans groupes latéraux (par exemple le poly-éthylène) peuvent s'aligner parallèlement les uns aux autres pour former des lamelles cristallines. Malgré tout, même dans ces conditions favorables, la cristallisation totale du fondu est difficile à obtenir (voir par exemple, [START_REF] Muthukumar | Nucleation in polymer crystallization[END_REF]). La difficulté intrinsèque de la cristallisation implique que les fondus de polymères sont en générale de bons candidats à la formation de verres [START_REF] Mckenna | Comprehensive Polymer Science[END_REF][START_REF] Angell | Relaxation in glassforming liquids and amorphous solids[END_REF][START_REF] Donth | The Glass Transition[END_REF].

Étant l'une des constantes élastiques, le module de cisaillement à l'équilibre (dans l'approximation de la réponse linéaire) G eq est une mesure directe de la transition du liquide (G eq = 0) vers l'état solide (G eq > 0). Cette observable décrit la réponse macroscopique du matériau aux faibles contraintes extérieures de cisaillement ; elle contient des informations à propos du paysage énergétique du système à l'équilibre. Dans l'article de Squire, Holt et Hoover une expression de cette constante élastique pour un système isotherme est dérivée à partir du potentiel d'interaction de paire de l'hamiltonien du modèle [START_REF] Squire | Isothermal elastic constants for argon: Theory and monte carlo calculations[END_REF]. Ils ont étendu la théorie classique de Born [START_REF] Born | Dynamical Theory of Crystal Lattices[END_REF] aux températures finies et ont trouvé un terme de correction à l'expression classique de Born pour les constantes élastiques, cette correction implique la moyenne des carrés des fluctuations de la contrainte (du "stress"). Ils ont nommé cette correction le "stress fluctuation formalism" et l'ont proposé comme une méthode adaptée au calcul des propriétés élastiques par les méthodes de simulation numérique sur ordinateur. Il était généralement admis que cette v méthode n'était uniquement valable que pour les solides élastiques et pour un champ de déplacement bien défini [START_REF] Squire | Isothermal elastic constants for argon: Theory and monte carlo calculations[END_REF][START_REF] Lutsko | Generalized expressions for the calculation of elastic constants by computer simulation[END_REF][START_REF] Frenkel | Understanding Molecular Simulation[END_REF], mais il a été démontré que cette approche pouvait aussi être appliquée pour toutes températures en passant par la transition vitreuse et que la valeur correcte de G eq = 0 pouvait être obtenue pour l'état liquide [START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary conditions, temperature, and sampling time[END_REF]. Pour les solides cristallins G eq disparait de manière discontine au point de fusion pour des températures T croissantes. Pour les solides amorphes et les verres, le comportement de G eq à la proximité de la température de transition vitreuse est matière à débat. Deux scénarios qualitativement différents ont été mis en avant, l'un postulant la singularité de ce saut discontinu à l'approche de la transition vitreuse [START_REF] Götze | Complex Dynamics of Glass-Forming Liquids: A Mode-Coupling Theory[END_REF][START_REF] Szamel | Emergence of long-range correlations and rigidity at the dynamic glass transition[END_REF][START_REF] Klix | Elastic properties of 2d amorphous solids[END_REF][START_REF] Klix | Discontinuous shear modulus determines the glass transition temperature[END_REF] ou l'autre postulant une transition continue "cusp-like" à T g [START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary conditions, temperature, and sampling time[END_REF][START_REF] Barrat | Elastic response of a simple amorphous binary alloy near the glass transition[END_REF][START_REF] Yoshino | Emergence of rigidity at the structural glass transition: A first-principles computation[END_REF][START_REF] Zaccone | Disorder-assisted melting and the glass transition in amorphous solids[END_REF].

Une autre situation intéressante émerge quand on confine les matériaux vitreux, par exemple, comme par confinement géometrique pour le cas des couches minces. Il est communément accepté que les propriétés des couches minces de polymères peuvent dévier substantiellement de celles du volume [START_REF] Alcoutlabi | Effects of confinement on material behaviour at nanometer size scale[END_REF][START_REF] Forrest | The glass transition in thin polymer films[END_REF]. Pourtant, la plupart des travaux expérimentaux se sont focalisés sur l'élucidation de la dépendance des propriétés thermiques des couches minces en fonction de l'épaisseur de la couche, plus précisement de la dépendance de la température de transition vitreuse (T g ) en fonction des interactions interfaciales et de la géométrie de la couche (e.g. couche mince en contact avec l'air ou déposé sur substrat). Souvent les résultats obtenus se sont montrés conflictuels [START_REF] Alcoutlabi | Effects of confinement on material behaviour at nanometer size scale[END_REF]. Du point de vue de la nanofabriquation et de sa conception, il serait utile de pouvoir déterminer comment les propriétés mecaniques sont affectées comparativement à la phase volumique.

Bien que de nombreuses méthodes expérimentales soient disponibles pour mesurer les propriétés mécaniques des couches minces de polymères, incluant les méthodes basées sur l'indentation, la dilatométrie, les ondes sonores de surface, l'angle de courbure sous contrainte, l'adaptation des ces méthodes de mesure expérimentale à la mesure des propriétés mécaniques des couches minces reste un défi technique. Pour parfaitement comprendre comment le confinement influe sur les propriétés mécaniques et la réponse à la contrainte pour les matériaux polymères, il est critique qu'un certain nombre de méthodes indépendantes aussi bien experimentales que numériques et théoriques soient mises en oeuvre dans la détermination des propriétés mécaniques des couches minces en fonction de l'epaisseur de la couche. Pour comprendre comment le confinement affecte les propriétés mécaniques et la réponse à la déformation d'une couche polymère ultra-mince, il est crucial qu'un certain nombre d'expériences soit réalisé (comme la déformation par création de nano-bulles [START_REF] O'connell | Novel nanobubble inflation method for determining the viscoelastic properties of ultrathin polymer films[END_REF]) et des simulations numériques pour mesurer les propriétés mécaniques vi en fonction de l'épaisseur du film soient réalisées. Nous présentons dans ce travail de thèse cette approche numérique en nous fondant sur les résultats de la phase volumique.

En se focalisant sur des modèles génériques -qui ne contiennent que les plus élémentaires fonctionalités d'un polymère (connectivité, volume exclu, etc...) -ceci donne la possibilité de capturer le comportement universel de la transition liquide-verre [START_REF] Schnell | Simulated glassforming polymer melts: Glass transition temperature and elastic constants of the glassy state[END_REF][START_REF] Yoshimoto | Mechanical heterogeneities in model polymer glasses at small length scales[END_REF][START_REF] Schulmann | Strictly two-dimensional self-avoiding walks: Thermodynamic properties revisited[END_REF]. Un exemple classique et largement utilisé de ceci est le modèle de collier de perles ("flexible bead-spring model" [START_REF] Baschnagel | Polymer Glasses, ch. Glass Transition and Relaxation Behavior of Supercooled Polymer Melts: An Introduction to Modeling Approaches by Molecular Dynamics Simulations[END_REF]). Dans ce modèle les liaisons sont représentées par un potentiel harmonique associé à une distance de liaison l [START_REF] Schnell | Simulated glassforming polymer melts: Glass transition temperature and elastic constants of the glassy state[END_REF]:

U bond (l) = 1 2 k bond (l -l 0 ) 2 , (1)
où k bond = 1110ǫ LJ /σ 2 LJ et l 0 = 0.967σ LJ sont la constante de force et la longueur d'équilibre de la liaison. Tous les autres monomères qui ne sont pas directement connectés interagissent à travers un potentiel Lennard-Jones (LJ) représentant les forces de Van der Waals:

U LJ (r) =    4ǫ LJ σ LJ r 12 -σ LJ r 6 + C, r < r cut 0 else, (2)
avec ǫ LJ et σ LJ étant la profondeur du puit de potentiel (minimum d'énergie) et le diamètre du monomère, respectivement [START_REF] Schnell | Simulated glassforming polymer melts: Glass transition temperature and elastic constants of the glassy state[END_REF].

Dans notre travail nous avons utilisé ce modèle et étudié l'influence du confinement sur les propriétés des fondus de polymères. Nous nous sommes concentré sur deux aspects des matériaux amorphes : la température de transition vitreuse et les propriétés rhéologiques.

Dans les deux cas nous avons commencé motre analyse par l'étude de la phase volumique.

Nous avons continué par les couches minces et nous nous sommes servis des résultats de la phase volumique comme point de référence. Pour ces deux systèmes nous avons aussi addressé la question de savoir si oui ou non G eq (T ) est discontinu à la température de transition vitreuse.

Température de transition vitreuse

En démarrant avec la phase volumique, nous mesurons la température de transition vitreuse T g via la dilatometrie, mais en utilisant la densité comme abscisse au lieu du volume.

Nous avons confirmé ensuite la dépendance en taille de chaîne T g (N ) = T ∞ La ligne pointillée rouge est une interpolation de Eq. ( 3) avec h 0 = 1.69. Médaillon:

Résolution couche par couche de la température de transition vitreuse T g (z * ) pour une couche mince avec h = 21.33, où z * est la position renormalisé dans l'épaisseur z * = z/h de refroidissement Γ [START_REF] Baschnagel | Polymer Glasses, ch. Glass Transition and Relaxation Behavior of Supercooled Polymer Melts: An Introduction to Modeling Approaches by Molecular Dynamics Simulations[END_REF]. En faisant cela nous avons obtenu un point de référence bien characterisé pour nos couches minces. Au regard du profil de densité des monomères ρ(z)

dans une couche mince nous avons approximé celui-ci analytiquement avec la fonction erf pour obtenir l'épaisseur du film h et le plateau de densité ρ 0 , qui correspond à la densité au milieu de la couche. Bien que T g pour les couches minces est habituellement déterminée à partir de la dépendance de l'épaisseur de la couche h en fonction de la température, nous avons utilisé le plateau de densité plutôt que h de telle manière à être cohérent avec l'analyse réalisé en phase volumique. Nous avons montré que T g dans les couches minces est plus bas que dans la phase volumique, et que cet abaissement est plus prononcé pour des valeurs de faible h. La dépendance de T g en fonction de l'épaisseur du film est bien décrit par [START_REF] Herminghaus | The glass transition of thin polymer films: some questions, and a possible answer[END_REF] (voir Fig. 1)

T g (h) = T bulk g (1 + h 0 /h) (3)
En enregistrant les profiles de densité monomèrique ρ(z) lors des refroidissements, nous avons aussi obtenu la résolution couche par couche "pseudo-thermodynamique" de T g (z).

En accord avec des expériences et simulations récentes [START_REF] Peter | Structure et relaxation structurale des fondus de polymères vitrifiables en couches minces[END_REF][START_REF] Baschnagel | Computer simulation of supercooled polymer melts in the bulk and in confined geometry[END_REF][START_REF] Mirigian | Theory of activated glassy relaxation, mobility gradients, surface diffusion, and vitrification in free standing thin films[END_REF] nous avons trouvé que T g (z)

décroit fortement en fonction de la distance à l'interface. Ceci prouve que la réduction de viii T g avec les valeurs décroissantes de h pour les couches minces à interface libre est liée à l'interface libre.

Propriétés mécaniques

En utilisant le "stress fluctuation formalism" [START_REF] Squire | Isothermal elastic constants for argon: Theory and monte carlo calculations[END_REF] nous nous sommes focalisé sur deux réponses : celle aux temps courts (ou fréquence infinie) du module de cisaillement µ a et celle aux temps longs ("equilibrium") du module de cisaillement G eq = µ a -µ f , où µ f mesure les fluctuations des contraintes.

Dépendance en fonction de la température

Nous avons démarré notre analyse à partir du cas de référence N = 4 dans la phase volumique. Nous avons trouvé que la fréquence d'échantillonnage à un effet faible sur la valeur de µ a mais devient crucial dans le cas de µ f . Nous avons fait le choix du temps d'échantillonnage ∆t = t max /10, pour un temps d'acquisition total de t max = 10 5 τ LJ Si nous regardons la dépendance de la température, nous obtenons que la réponse à temps court, µ a (T ), peut être décrit de manière correcte en fonction de la densité comme µ a (T ) = Aρ 3 (T ) + Bρ(T ). Le potentiel LJ est identifié comme le responsable de la dépendance en ρ 3 , quand la partie proportionnelle à ρ est associé avec le potentiel de liaison. Une dépendance similaire avec la densité a aussi été trouvée pour le module élastique d'un liquide de Lennard-Jones au niveau de la coexistence solide-liquide [START_REF] Heyes | Scaling of lennardjones liquid elastic moduli, viscoelasticity and other properties along fluidsolid coexistence[END_REF].

Pour la réponse à temps long, le module de cisaillement à l'équilibre G eq est identifié comme étant une mesure sans ambiguité de la transition à partir du liquide. Celui-ci montre une transition franche (sigmoidale), de l'état liquide G eq = 0 vers l'état solide avec G eq > 0. La forme de cette transition est dépendante du temps d'échantillonnage. Nous insistons sur le fait que la transition devient plus franche quand le temp d'échantillonnage augmente, mais reste continu. En plus, nous avons montré qu'il est possible de déterminer G eq à partir des déplacements carrés moyens des monomères ou de correlation de déplacement de ces mêmes monomères [START_REF] Klix | Elastic properties of 2d amorphous solids[END_REF][START_REF] Klix | Discontinuous shear modulus determines the glass transition temperature[END_REF][START_REF] Schnell | Simulated glassforming polymer melts: Glass transition temperature and elastic constants of the glassy state[END_REF][START_REF] Fuchs | Aspects of the dynamics of colloidal suspensions: Further results of the mode-coupling theory of structural relaxation[END_REF]. données pour les cinq tailles de chaînes en traçant le module de cisaillement à l'équilibre en unité réduite Geq = G eq (T )/ρk B T g comme une fonction de la température en unité réduite T = T /T g (voir Fig. 2).

Dépendance en longueur de chaîne

x

Dépendence en épaisseur du film

En ce qui concerne les films, nous avons commencé l'analyse de µ a (T ) pour quatre épaisseurs h(T g ) = 21.33, 8.50, 6.61, 4.82 avec la longueur de chaîne N = 16. Nous avons trouvé que le module de cisaillement à fréquence infinie ne dépend pas de h dans l'état vitreux, en accord avec [START_REF] Forrest | Brillouin light scattering studies of the mechanical properties of thin freely standing polystyrene films[END_REF]. Par contre, ce module augmente avec h dans l'état liquide.

En déterminant le module local, µ a (z), nous avons démontré que cette augmentation est reliée à la baisse de µ a (z) proche de l'interface. Nous avons démontré que µ a (z) peut être complètement décrit en utilisant la dépendence en densité du module en volume :

µ a (z) = µ bulk a (ρ(z)).
En utilisant µ bulk a (ρ) ∼ ρ 3 + ρ, nous pouvant déduire la dépendance en h de µ a . Nous soulignons que le terme en ρ 3 est responsable de la baisse de µ a en diminuant h. Nous avons aussi constaté que cette baisse est masquée dans l'état vitreux par l'augmentation de la densité du plateau (ceci est en contradiction avec [START_REF] Forrest | Brillouin light scattering studies of the mechanical properties of thin freely standing polystyrene films[END_REF] qui suppose la densité étant indépendant de h).

Nous avons montré ensuite la dépendance en épaisseur du module d'équilibre G eq (T ).

Celui-ci décroît en diminuant h ; ceci est conforme avec des résultats de littérature sur des films polymères amorphes [START_REF] Stafford | Elastic moduli of ultrathin amorphous polymer films[END_REF] ainsi que pour des films cristallins en phase fcc [START_REF] Van Workum | Local elastic constants in thin films of an fcc crystal[END_REF].

Nous avons mis en évidence que cette dépendance en épaisseur provient de la dépendance en pression de G eq (T ). C'est donc une conséquence des détails de la préparation du film, qui dans notre cas (où nous fixons la surface de la boite de simulation) induit des contraintes tangentielles, P T , dans le film. Pour démontrer ceci, nous sommes revenus au système volumique et avons effectué des simulations pour différentes pressions P à T = 0.2 < T g . Nous avons donc obtenu G eq en fonction de la pression qui se trouve être linéaire : G eq (P ) = G eq (0)[1 + AP ], avec A ≈ 0.11. En traçant G eq (T )/(1 + AP T ) vs T /T g pour nos quatre films, il est possible d'obtenir une courbe maitresse avec toutes les données (voir Fig. 3).

Comparaison à la littérature

Enfin, nous avons comparé nos résultats de G eq pour le massif (bulk) et les systèmes de films minces à d'autre systèmes venant de la littérature. Nous avons considéré deux systèmes modèles, le modèle de Kob-Andersen en 3D (3dKA) [START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary conditions, temperature, and sampling time[END_REF] et un modèle en 2D de particules Lennard-Jones polydisperse (2dpLJ) [START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary conditions, temperature, and sampling time[END_REF] ainsi qu'un système expérimental, un mélange 2D binaire de colloïdes [START_REF] Klix | Discontinuous shear modulus determines the glass transition temperature[END_REF]. Nous avons obtenu une similarité remarquable dans le comportement de Geq . En se qui concerne le saut de G eq à T g , nous observons que ce saut est presqu'identique pour tous les systèmes. Néanmoins, le comportement dans l'état de verre est dépendant de la dimensionnalité du système considéré.

xii

Perspectives

Bien que nos travaux apporte des éclaircissements sur l'influence de la surface libre sur les propriétés des films de polymères formant un verre, il reste beaucoup de questions ouvertes. Des études de simulations supplémentaires pourrait traiter les questions suivantes.

Pour le système massif (bulk) :

• Le travail ref. [START_REF] Heyes | Scaling of lennardjones liquid elastic moduli, viscoelasticity and other properties along fluidsolid coexistence[END_REF] traite des lois d'échelle du module de cisaillement. Les auteurs suggèrent que la puissance de la densité est obtenue un utilisant les distances adimensionnées r = rρ -1/3 . Dans la ref. [START_REF] Heyes | Scaling of lennardjones liquid elastic moduli, viscoelasticity and other properties along fluidsolid coexistence[END_REF], ils montrent que la partie repulsive du potentiel est à l'origine de la dépendance en ρ 5 et la partie attractive en une dépendance en ρ 3 de µ a (pour le potentiel 12-6 LJ). Puisque le formalisme des fluctuations de contraintes donne la possibilité de séparer directement µ a en toutes les contributions possibles, il est possible de séparer la partie LJ en parties "attractives" et "répulsives" . Ceci montreras l'origine du facteur d'échelle observé en ρ 3 de µ a , et répondras à la question de savoir si le facteur d'échelle en ρ 5 + ρ 3 suggéré dans la ref. [START_REF] Heyes | Scaling of lennardjones liquid elastic moduli, viscoelasticity and other properties along fluidsolid coexistence[END_REF] est également possible dans notre cas. Pour prolonger cela, il serait intéressant de sonder directement l'influence de la puissance du potentiel (U (r) ∝ 1 r n ) sur la puissance du facteur d'échelle de densité de µ a .

• Pour l'analyse de µ a , il serait utile de séparer G eq en contributions de LJ et contributions des liaisons ou, plus pratique, en contributions intra-et inter-chaînes.

Pour les films minces:

• Il serait intéressant de verifier la relation µ a (z) = µ bulk a (ρ(z)) et si elle s'applique pour les films simulés et pour des polymères "réels".

• G eq (z) résolu localement par couche doit être calculé par le formalisme des contraintesfluctuations. Ceci donnera une compréhension plus profonde sur la dépendance en h. En plus, ceci permettra une étude plus précise de la forme de G eq (T ) au seuil de la transition vitreuse.

• Implémenter la méthode des champs de déplacement pour les films, puisque celleci est reliée à la complaisance de fluage en cisaillement. Implémenter ensuite la déformation biaxiale également pour avoir une connection directe avec l'expérience de "nanobubble inflation" [START_REF] O'connell | Novel nanobubble inflation method for determining the viscoelastic properties of ultrathin polymer films[END_REF].

xiii Chapter 1

Introduction

Polymer melts are liquids composed of macromolecular chains [START_REF] Doi | The Theory of Polymer Dynamics[END_REF][START_REF] Rubinstein | Polymer Physics[END_REF]. A macromolecule can have various microstructures (homopolymer, copolymer, etc.) and architectures (linear, ring, star, etc.). We will only be concerned with the simplest structure, that of a linear homopolymer, where N monomeric repeat units of the same type are connected to form a chain. Moreover, we will assume that the polymers in the melt are monodisperse (all have the same N ), that their chain conformations are flexible, and that the nonbonded interactions are short-range and spatially isotropic. This is a standard chain model pertinent to many synthetic polymers [START_REF] Rubinstein | Polymer Physics[END_REF]. Experimentally, the chain length N can be large, a typical range being 10 3 N 10 5 . Since the average size of a polymer grows with N , these large chain lengths also imply the chain size to be large. A possible measure of the chain size is the (average) end-to-end distance (R) of a polymer [START_REF] Doi | The Theory of Polymer Dynamics[END_REF][START_REF] Rubinstein | Polymer Physics[END_REF]. For 10 3 N 10 5 , the end-to-end distance is typically in the range of R ≈ 100 -1000 Å. The size of a chain thus exceeds that of a monomer (≈ 5 Å ) by several orders of magnitude.

These widely distributed length scales manifest themselves in the structural and dynamic properties of a polymer melt. In the melt, the monomers spatially arrange in densely packed nearest-neighbor shells. This collective structure of the melt, does not change when N increases; it remains liquid-like. By contrast, the dynamics is strongly affected. With increasing N , the monomer relaxation time τ 1 and the chain relaxation time τ N become more and more separated. This opens a time window τ 1 < t < τ N where slow relaxation of the polymers occurs, before viscous flow sets in for t > τ N . Therefore, triggered by the increase of chain length, viscoelasticity emerges in polymer liquids, already at high temperature.

Glass formation

When cooled to low temperatures, the polymer liquid transforms into a solid that can be either amorphous (glassy) or semicrystalline [START_REF] Strobl | The Physics of Polymers: Concepts for Understanding their Structures and Behavior[END_REF]. Semicrystalline polymers have both amorphous and crystalline regions. In many cases, the crystalline regions are formed of lamellar sheets in which the polymers fold back and forth so that chain sections align parallel to each other. The sheets twist and branch as they grow outward from a nucleus, forming spherulitic structures at large scales [START_REF] Strobl | The Physics of Polymers: Concepts for Understanding their Structures and Behavior[END_REF]. The hierarchy of these morphological features range from the lamellar ordering of the chains (10 nm) to the macroscopic packing of the spherulites (100 µm and larger). This reflects the complexity of the underlying crystallization process, which is not fully understood [START_REF] Muthukumar | Nucleation in polymer crystallization[END_REF][START_REF] Luo | Frozen topology: Entanglements control nucleation and crystallization in polymers[END_REF].

The ability to form crystals hinges on the microstructure of the polymer. Only chains with regular configurations (for example, isotatic or syndiotatic orientations of the sidegroups [START_REF] Rubinstein | Polymer Physics[END_REF]) or chains without sidegroups (for example, polyethylene) can align parallel to each other so as to pack into crystalline lamellae. However, even in these favorable cases, full crystallization is hard to achieve (see, for example, [START_REF] Muthukumar | Nucleation in polymer crystallization[END_REF]). This intrinsic difficulty of crystal formation implies that polymer melts are in general good glass formers [START_REF] Mckenna | Comprehensive Polymer Science[END_REF][START_REF] Angell | Relaxation in glassforming liquids and amorphous solids[END_REF][START_REF] Donth | The Glass Transition[END_REF]. Either they can be readily supercooled or, due to irregular chain configurations, a crystalline phase does not exist at all. There are several examples for the latter case, for instance, homopolymers with an atactic orientation of (bulky) sidegroups (atactic polystyrene, etc.) or random copolymers, such as cis-trans polybutadiene, in which monomers, having the same chemical composition, but different microstructures (cis/trans configuration of butadiene), are randomly concatenated.

Properties of glasses

These polymeric glass formers exhibit several properties that are also found in other nonpolymeric (intermediate and fragile [START_REF] Böhmer | Nonexponential relaxations in strong and fragile glass formers[END_REF][START_REF] Berthier | Theoretical perspective on the glass transition and amorphous materials[END_REF]) glass-forming liquids and are thus typical of the glass transition in general [START_REF] Angell | Relaxation in glassforming liquids and amorphous solids[END_REF][START_REF] Donth | The Glass Transition[END_REF][START_REF] Berthier | Theoretical perspective on the glass transition and amorphous materials[END_REF][START_REF] Cavagna | Supercooled liquids for pedestrians[END_REF][START_REF] Biroli | Perspective: The glass transition[END_REF][START_REF] Ediger | Perspective: Supercooled liquids and glasses[END_REF]. As the liquid is cooled toward low temperatures, it becomes kinetically arrested below a characteristic temperature T g and forms an amorphous solid for T ≪ T g . Often, T g is determined by heat capacity measurements and is therefore also referred to as "calorimetric glass transition temperature" (see Panel (a) of Fig. 1.1). But it is not the only method. Several other thermodynamic, physical, mechanical and electrical properties have a distinct change close to T g . For example, we can look at the volume of the system along the cooling. In contrast to crystalliza- Along with that, commonly measured structural properties only change smoothly upon the cooling through the glass transition. For instance, the structure of the glass and the liquid, characterized by the static structure factor S(q), is essentially identical. A glass just appears to be a "frozen liquid", a liquid that has stopped to flow. This suggests to focus on the dynamic behavior and its T dependence as the glass transition is approached from the high-temperature side. Indeed, on cooling toward T g , characteristic dynamic features, sometimes collectively referred to as glassy dynamics, emerge. Certainly, the most prominent feature is the massive increase of all structural relaxation times, over about 14 orders of magnitude, from ≈ 10 -12 s at high T to ≈ 10 2 s at T g (see Panel (c) of

Mechanical response

Being one of the elastic constants, the equilibrium shear modulus G eq describes the macroscopic response of the solid to weak external deformations and encodes information about the potential energy landscape of the system. Establishing the link between the elastic constants and the features of the potential energy landscape is a first step in understanding the mechanical behavior of the system. In a seminal work, Squire, Holt and Hoover derived an expression for the (isothermal) elastic constants in terms of the (pairwise) interaction potential [START_REF] Squire | Isothermal elastic constants for argon: Theory and monte carlo calculations[END_REF]. By extending the classical Born theory [START_REF] Born | Dynamical Theory of Crystal Lattices[END_REF] to finite temperature they found a correction term to the Born expression for the elastic constants, which involves the mean-square fluctuations of the stress. They proposed this "fluctuation formalism" as a convenient route to calculating the elastic properties in computer simulations. Since then, many simulation studies have utilized this approach and it has been generalized

to systems with nonzero initial stress [START_REF] Ray | Molecular dynamics calculation of elastic constants for a crystalline system in equilibrium[END_REF], hard-sphere interactions [START_REF] Farago | Fluctuation formalism for elastic constants in hardspheres-and-tethers systems[END_REF] and arbitrary continuous potentials [START_REF] Lutsko | Generalized expressions for the calculation of elastic constants by computer simulation[END_REF], or to the calculation of local mechanical properties [START_REF] Yoshimoto | Mechanical heterogeneities in model polymer glasses at small length scales[END_REF].

In particular from Ref. [START_REF] Lutsko | Generalized expressions for the calculation of elastic constants by computer simulation[END_REF] it has become clear that the fluctuation correction to the Born term does not necessarily vanish in the zero-temperature limit. This is due to the fact that when a system is subjected to a homogeneous deformation, the ensuing particle displacements need not follow the macroscopic strain affinely; there may be nonaffine displacements. The fluctuation term quantifies the extent of the nonaffine contributions, whereas the Born term reflects the affine part of the particle displacement. How important the nonaffine displacements are, depends on the system under consideration [START_REF] Lutsko | Generalized expressions for the calculation of elastic constants by computer simulation[END_REF][START_REF] Barrat | Elastic response of a simple amorphous binary alloy near the glass transition[END_REF]. While the elastic properties of crystals with one atom per unit cell are given by the Born term only, stress fluctuations are significant for crystals with more complex unit cells [START_REF] Lutsko | Generalized expressions for the calculation of elastic constants by computer simulation[END_REF].

They become particularly pronounced for amorphous solids [START_REF] Maloney | Universal breakdown of elasticity at the onset of material failure[END_REF]. This is revealed by several recent simulation studies on various glass formers, like Lennard-Jones mixtures [START_REF] Tanguy | Continuum limit of amorphous elastic bodies: A finite-size study of low-frequency harmonic vibrations[END_REF][START_REF] Léonforte | Continuum limit of amorphous elastic bodies (III): Three dimensional systems[END_REF], polymer glasses [START_REF] Yoshimoto | Mechanical heterogeneities in model polymer glasses at small length scales[END_REF] or silica melts [START_REF] Léonforte | Inhomogeneous elastic response of silica glass[END_REF].

Many of these simulation studies concentrate on the mechanical properties deep in the glassy state, exploring for instance correlations between the nonaffine displacement field and vibrational anomalies of the glass ("Boson peak") [START_REF] Léonforte | Continuum limit of amorphous elastic bodies (III): Three dimensional systems[END_REF][START_REF] Léonforte | Inhomogeneous elastic response of silica glass[END_REF], the mechanical heterogeneity at the nanoscale [START_REF] Yoshimoto | Mechanical heterogeneities in model polymer glasses at small length scales[END_REF][START_REF] Tsamados | Local elasticity map and plasticity in a model Lennard-Jones glass[END_REF] or the onset of molecular plasticity in the regime where the macroscopic deformation is still elastic [START_REF] Papakonstantopoulos | Molecular plasticity of polymeric glasses in the elastic regime[END_REF].

Since the equilibrium shear modulus G eq (T ) of crystalline solids vanishes discontinuously at the melting point with increasing temperature T , this naturally begs the question of the behavior of G eq (T ) for amorphous solids and glasses in the vicinity of the glass transition temperature T g [START_REF] Barrat | Elastic response of a simple amorphous binary alloy near the glass transition[END_REF]. Two qualitative different scenarios have been put forward being either in favor of a discontinuous jump singularity at the glass transition [START_REF] Götze | Complex Dynamics of Glass-Forming Liquids: A Mode-Coupling Theory[END_REF][START_REF] Szamel | Emergence of long-range correlations and rigidity at the dynamic glass transition[END_REF][START_REF] Klix | Elastic properties of 2d amorphous solids[END_REF][START_REF] Klix | Discontinuous shear modulus determines the glass transition temperature[END_REF] or a continuous (cusp-like) transition [START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary conditions, temperature, and sampling time[END_REF][START_REF] Barrat | Elastic response of a simple amorphous binary alloy near the glass transition[END_REF][START_REF] Yoshino | Emergence of rigidity at the structural glass transition: A first-principles computation[END_REF][START_REF] Zaccone | Disorder-assisted melting and the glass transition in amorphous solids[END_REF].

Thus it is important to study the mechanical behavior of the glasses in the wide temperature range that covers both the deep glassy state and the glass transition region.

Glass-forming polymer films

Another interesting situation arises when one somehow confines the glassy material, for example in a thin film geometry. It is widely accepted that the physical properties of thin glass-forming polymer films can deviate substantially from their bulk counterparts [START_REF] Alcoutlabi | Effects of confinement on material behaviour at nanometer size scale[END_REF][START_REF] Forrest | The glass transition in thin polymer films[END_REF]. Most experimental studies have focused on elucidating the thickness dependence of the thermal properties in thin polymer films, specifically the glass transition temperature (T g ), as a function of interfacial interactions and molecular entanglements, often with conflicting results [START_REF] Alcoutlabi | Effects of confinement on material behaviour at nanometer size scale[END_REF]. However, from an application perspective it would be helpful to directly determine whether the mechanical properties display similar deviations from the bulk response. While a number of experimental methods are available for measuring the mechanical properties of polymer films, including methods based on indentation [50,51],

surface acoustic waves [52], and beam curvature [53], adapting these measurement techniques to ultrathin polymer films remains challenging. To understand how confinement affects the mechanical properties and response of such ultrathin polymeric materials, it is crucial that a range of experimental (like the nanobubble inflation method [START_REF] O'connell | Novel nanobubble inflation method for determining the viscoelastic properties of ultrathin polymer films[END_REF]) and computational methods for measuring the mechanical properties of ultrathin films as a function of film thickness are made available. We present here such a computational approach building on the results obtained for the bulk.

Outline

Focusing on generic models -that only retain the most basic features of a polymer as connectivity and excluded volume -gives a possibility to grasp a universal behavior of the liquid-glass transition [START_REF] Schnell | Simulated glassforming polymer melts: Glass transition temperature and elastic constants of the glassy state[END_REF][START_REF] Yoshimoto | Mechanical heterogeneities in model polymer glasses at small length scales[END_REF][START_REF] Schulmann | Strictly two-dimensional self-avoiding walks: Thermodynamic properties revisited[END_REF]. A classical and widely utilized example is the flexible bead-spring model [START_REF] Baschnagel | Polymer Glasses, ch. Glass Transition and Relaxation Behavior of Supercooled Polymer Melts: An Introduction to Modeling Approaches by Molecular Dynamics Simulations[END_REF].

In our work we have employed this model and studied the influence of confinement on the properties of the polymer melts. We focus on two important aspects of amorphous materials: the glass transition temperature and the mechanical properties. In both cases we start the analysis by looking at the bulk system. We continue with film systems and use bulk results as the reference point.

The manuscript is organized as follows:

• Chapter 2 introduces the simulation methodology and the theoretical background of the analysis. We briefly describe the used simulation technique and the simulation model. Further, we present the simulation procedure with the times and other parameters used for each step. In the end, we discuss different methods to sample viscoelastic properties of the system.

• Chapter 3 is dedicated to the glass transition temperature. In the first part we focus on the chain length and cooling rate dependence of the glass transition in the bulk system. In the second part we analyze the film thickness dependence of the glass transition temperature.

• Chapter 4 presents the result on viscoelastic properties. Starting with the bulk system, we first investigate the temperature dependence of the short and long time response functions. In this context we compare different methods and discuss the effect of the finite time of the simulation. We continue with the bulk system, but shift our focus to the influence of the chain length on the viscoelastic properties. Further, we use the results on the bulk system to explain the dependence of viscoelastic properties on the thickness of the film. In the end, we compare our results for the bulk and the film systems to a different simulated models and recent experiment.

• Chapter 5 gives a summary of our work along with the future perspectives.

List of publications:

Along the thesis I have contributed to several publications (you can find the abstracts in Appendix B).

Three articles below are dedicated to the development and analysis of computational techniques that give the possibility to measure the time-dependent shear relaxation modulus G(t). We also discuss its connection to the equilibrium shear modulus G eq (which is one of the main quantities that we focus on in Chapter 4 of this thesis).

• 2015 Shear-strain and shear-stress fluctuations in generalized Gaussian ensemble simulations of isotropic elastic networks.

m d 2 r i (t) dt 2 = F i (t), (2.1)
where F i denotes the total force on particle i and derives from the potential U N (r 1 , . . . r N ) describing the interections:

F i = - ∂U N (r 1 , . . . r N ) ∂r i (2.2)
If the ergodic principle holds [START_REF] Frenkel | Understanding Molecular Simulation[END_REF][START_REF] Allen | Computer Simulation of Liquids[END_REF], the set of microscopic configurations obtained from MD simulation and compatible with macroscopic constrains, such as fixed volume, temperature, etc., defines the corresponding statistical ensemble. Thus all physical properties of the system could be calculated as time averages that correspond to ensemble averages.

Originally, MD simulation was invented in the microcanonical ensemble that corresponds to the macroscopic constraints of constant number of particles N , constant volume V and constant total energy E. It is possible to extend MD simulation to other ensembles by introducing additional degrees of freedom. This will be discussed in section 2.1.2.

Integration of the equations of motion

To simulate the time evolution of the system, MD numerically integrates Equation (2.1) in an iterative way. Lets assume that we know the positions and velocities for all particles at time t. The positions r i (t) allow us to determine the forces F i (t). Then, a Taylor expansion gives an estimate of the new positions and velocities at a small time increment ∆t later:

r i (t + ∆t) ≈ r i (t) + v i (t)∆t + F i (t) 2m (∆t) 2 (2.3) v i (t + ∆t) ≈ v i (t) + F i (t) m ∆t (2.4)
The time increment ∆t is referred to as the "time step" in MD simulations. From r i (t + ∆t) we can obtain the new forces F i (t + ∆t) (Eq. (2.2)), which completes the input for Equations (2.3) and (2.4) to calculate the positions and velocities at time t+2∆t. Iteration of this procedure therefore furnishes a discretized trajectory of the system, x(t k = k∆t), with k = 0, 1, . . . , N max , starting from the initial configuration (k = 0) up to the final configuration for the maximum number N max of time steps simulated.

While the preceding description illustrates the principle of the MD simulation, Equations (2.3-2.4) are not used in practice due to the following drawback: They are not invariant with respect to time reversal, whereas Equation (2.1) is. Time-reversal symmetry implies that the structure of the equations of motion is preserved, if we move backward in time. Assume that we know the positions and velocities at time t+∆t and we want to calculate r i (t) and v i (t) by inserting the reverse time step -∆t in Equations (2.3) and (2.4).

This would alter the equations, thereby violating time-reversal symmetry.

To remedy this problem, assume again that we know the positions and velocities at time t + ∆t. From the positions, we calculate the new forces F i (t + ∆t). Then, we can apply the Taylor expansion of the velocities in reverse:

v i (t) ≈ v i (t + ∆t) - F i (t + ∆t) m ∆t ⇒ v i (t + ∆t) ≈ v i (t) + F i (t + ∆t) m ∆t
Adding this result to Equation (2.4) gives

v i (t + ∆t) ≈ v i (t) + F i (t) + F i (t + ∆t) 2m ∆t (2.5)
Now Equation (2.5) is time reversible because it is symmetric with respect to t and t + ∆t, contrary to Equation (2.4). Moreover, Equation (2.5) also makes the positions time reversible. To see this, lets apply the reverse time step -∆t from t + ∆t to t in Equation (2.3):

r i (t) ≈ r i (t + ∆t) -v i (t + ∆t)∆t + F i (t + ∆t) 2m (∆t) 2
Inserting Equation (2.5) , we get

r i (t) ≈ r i (t + ∆t) -v i (t)∆t - F i (t) 2m (∆t) 2
which can be rearranged to give back Equation (2.3). The resulting algorithm

r i (t + ∆t) ≈ r i (t) + v i (t)∆t + F i (t) 2m (∆t) 2 v i (t + ∆t) ≈ v i (t) + F i (t) + F i (t + ∆t) 2m ∆t (2.6)
is called "velocity-Verlet algorithm" [START_REF] Frenkel | Understanding Molecular Simulation[END_REF][START_REF] Allen | Computer Simulation of Liquids[END_REF]. It is the most commonly used algorithm in

MD simulations

The preceding description of the MD method shows that the time step ∆t is a crucial operational parameter. Certainly, we would like to take ∆t as large as possible, since this allows to extend the longest simulation time, t max = N max ∆t, at fixed computational effort N max . However, the accuracy of the simulation may then suffer because the Taylor expansion in Equation (2.6) is only appropriate for small ∆t.

Simulating different ensembles: NVE, NVT, NPT

Solving Newtons equations of motion we are bound to energy conservation and thus can explore the phase space of the microcanonical ensemble (NVE). However, experiments typically control the temperature (T ) instead of total energy and external pressure (P ) instead of volume. That is why MD simulation were extended to be able to realize NVT and NPT ensembles, which is usually achieved by introducing new dynamical variables [START_REF] Frenkel | Understanding Molecular Simulation[END_REF][START_REF] Allen | Computer Simulation of Liquids[END_REF].

NVT ensemble

The canonical ensemble represents the system in thermal equilibrium with a larger system (heat bath) at a fixed temperature. The system can exchange energy with the heat bath, thus the energy in the system is not conserved. The heat bath is assumed to be so large that it is not affected by the energy transfer. However, direct introduction of the heat bath in Newton's equations would give rise to an enormous amount of new degrees of freedom, that makes the already complicated N -body problem not feasible.

Nosé and Hoover resolved this problem by collecting all these additional degrees of freedom into one new dynamical variable ξ, giving rise to the equations of motion:

ṙi = p i m (2.7) ṗi = F i -ξp i (2.8) ξ = 3N k B Q T (t) -T ext (2.9)
where p i -the generalized momentum of particle i, T ext -the external temperature of the heat bath, and T is the instantaneous temperature of the system, related to the kinetic energy via the equipartition-theorem-like relation:

E k = p 2 i 2m = 3 2 N k B T (t) (2.10)
The additional term -ξp controls the exchange of the kinetic energy between the system and the bath. If ξ is negative, it accelerates particles. If it is positive, particles are slowed down. The change of ξ is controlled by the difference between the temperature of the system and the temperature of the bath. In another words, if the temperature of the system is higher than the temperature of the bath, ξ increases (slowing down the particles) and vice versa. The strength of the coupling between the bath and the system is determined by the thermal inertia coefficient Q.

NPT ensemble

The control of the pressure in the system can be realized in a similar way as the control of the temperature. For doing so, the volume V of the simulation box is allowed to fluctuate in order to keep the pressure constant. In an isotropic system with central forces the instantaneous pressure can be calculated as:

P (t) = 1 3V i mv 2 i + i<j F ij • r ij , (2.11)
where F ij is the force acting between particles i and j, r ij -the distance between these particles.

Introducing a new dynamical variable η, the equations of motion in the NPT ensemble are given by ṙi =

p i m + η (r i -R cm ) (2.12) ṗi = F i -(ξ + η)p i (2.13) ξ = 3N k B Q T (t) -T ext (2.14) η = V M (P (t) -P ext ) (2.15) V = 3V η (2.16)
with R cm is the system's center of mass.

The time evolution of the volume depends on η, which is coupled to the deviation of the system's instantaneous pressure P (t) from the imposed external pressure P ext . The coupling parameter M determines the frequency of the volume fluctuations. In the case of M → ∞ the canonical ensemble is recovered. The position of center of mass is subtracted in Eq. (2.12) to guarantee that the NPT ensemble is realized [START_REF] Melchionna | [END_REF].

The exact values of Q and M used in our simulation are presented in Sect. 2.3.

Polymer model

Here we use a generic flexible bead-spring model for linear homopolymer chains [START_REF] Schnell | Simulated glassforming polymer melts: Glass transition temperature and elastic constants of the glassy state[END_REF][START_REF] Peter | Modeling dielectric relaxation in polymer glass simulations: Dynamics in the bulk and in supported polymer films[END_REF].

The bonds are represented by a harmonic potential with the bond length l:

U bond (l) = 1 2 k bond (l -l 0 ) 2 , (2.17)
where k bond = 1110ǫ/σ 2 and l 0 = 0.967σ are the force constant and the equilibrium bond length. All other monomers that are not connected by bonds interact by a Lennard-Jones (LJ) potential:

U LJ (r) =    4ǫ σ r 12 -σ r 6 + C, r < r cut 0 else, (2.18)
where ǫ is the depth of the potential minimum and σ the particle diameter. The constants ǫ and σ set respectively a scale for the energy and the length. The potential is truncated to increase numerical efficiency and shifted to make it smooth (continuous with its 1st

derivative, but not with its 2nd which gives additional corrections to some observables [START_REF] Peter | Modeling dielectric relaxation in polymer glass simulations: Dynamics in the bulk and in supported polymer films[END_REF]). The cutoff radius r cut = 2.3σ ≃ 2r min , where r min is the minimum of the LJ potential.

Like other bead spring models [START_REF] Bennemann | Molecular-dynamics simulations of the thermal glass transition in polymer melts -alpha-relaxation behavior[END_REF], the present model effectively eliminates crystalization. It is the consequence of two properties of such models. First, the distance r min of the LJ potential is incommensurate with l 0 of the bond potential. Thus the bond potential locally distorts possible crystalline arrangements of the LJ monomers (fcc or bcc).

Second and very crucial, the polymer chains are very flexible. Only the repulsive part of the LJ potential supresses back-folding of adjacent bonds [START_REF] Baschnagel | Computer simulation of supercooled polymer melts in the bulk and in confined geometry[END_REF].

Lennard-Jones units and approximate mapping to SI units

The constants ǫ and σ in Eq. (2.18) set respectively a scale for the energy and the length.

So, we can introduce nondimensional (reduced) quantities:

U * LJ = U LJ /ǫ, the distance r * = r/σ, time t * = t/τ LJ with τ LJ = mσ 2 ǫ
and temperature T * = T /(ǫ/k B ). It is common practice in simulation to use so called "LJ units" [START_REF] Allen | Computer Simulation of Liquids[END_REF]. In this case all fundamental quantities are set to 1 (ǫ = 1, σ = 1, m = 1 and k B = 1). Then, r * = r, t * = t, etc., and we can drop the "*" to simplify the notation. We use LJ units for all quantities presented in following sections.

The mapping of LJ units to SI units has been discussed in several works [ and polybutadiene with that of a beadspring model [START_REF] Paul | Structure and dynamics of amorphous polymers: computer simulations compared to experiment and theory[END_REF]. The result is τ LJ ≃ 2.1 × 10 -13 s.

However, these values for σ, ǫ and τ LJ depend on the real polymer for which the conversion is carried out, and may vary a lot. This caveat is discussed in [59,[START_REF] Kröger | Simple models for complex nonequilibrium fluids[END_REF][START_REF] Peter | Multiscale simulation of soft matter systems -from the atomistic to the coarse-grained level and back[END_REF]. In their seminal work on reptation dynamics [59] Kremer and Grest reported that an LJ bead corresponds to the range of 1/2 to 5 monomers for real polymers. This variation reflects differences in size and flexibility of the specific monomer, and entail fluctuations in the results for σ, ǫ and τ LJ . Typically, they fall in the range: 5 × 10 -10 m σ 13 × 10 -10 m, 300 K ǫ/k B 500 K and 2 × 10 -12 s τ LJ 2 × 10 -10 s [59].

Due to these uncertainties we will make the following choices when giving an approximate conversion from LJ to physical units: σ = 5 × 10 -10 m, ǫ/k B = 450 K and τ LJ = 2 × 10 -12 s. This implies a monomer mass of m ≃ 60 g/mol. With that we also obtain reference values for other quantities, e.g., a mass density of 1 corresponds to 0.8 g/cm 3 and a pressure of 1 corresponds to 497 bar. Since a pressure of 1 bar is very small in LJ units (P ≃ 0.002), vanishing reduced pressure (P = 0) is a good proxy for ambient pressure conditions.

Simulation aspects

The simulations are performed in two different ensembles: NVT and NPT. For fixing the temperature and/or pressure we use the Nosé-Hoover algorithm provided by LAMMPS ("fix nvt" or "fix npt") [START_REF]Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS)[END_REF]. The damping coefficients for the barostat is P damp = 75 and for the thermostat T damp = 1. These parameters related to Q and M from Sect. 2.1.2 via

Q = 3N k B T × T 2 damp M = N k B T × P 2 damp
where N is the total number of particles in the system .

The equations of motion are integrated using the velocity Verlet algorithm (see Sect. 2.1.1), with time step ∆t = 0.005. This value of the time step is the usual choice when using LJ units, it gives a good balance between the accuracy and the maximum simulation time. All simulations consist of three parts: continuous cooling, "equilibration" and the sampling of the properties.

Bulk

We The number of particles is enough to expect continuum elasticity to be applicable [START_REF] Schnell | Simulated glassforming polymer melts: Glass transition temperature and elastic constants of the glassy state[END_REF].

The initial configurations ave been taken from the PhD thesis of Stephen Frey [START_REF] Frey | Propriétés viscoélastiques des fondus de polymères vitrifiables[END_REF].

They are equilibrated at temperature T i = 0.7 well above the glass transition temperature T g (see below) and at pressure P = 0. As the criterion for the equilibration the orientational correlation function of end-to-end vector φ e (t) were taken

φ e (t) = R e (t) • R e (0) R 2 e (0) (2.19)
where R e (t) is the end-to-end vector. Since φ e (t) measures the slowest relaxation process in terms of the chains [START_REF] Doi | The Theory of Polymer Dynamics[END_REF], the system is considered to be equilibrated when φ e (t) drops to the value of 0.1.

Keeping the pressure constant, the imposed mean temperature is linearly decreased according to:

T (t) = T i -Γt, (2.20)
where Γ is the constant cooling rate and t is the time. Here we used the range of cooling rates Γ from 2 × 10 -6 to 2 × 10 -4 . This range was used to explore the dependence of the glass transition temperature on Γ (see Sect. 3.1.2). For all other purposes we mainly used Γ = 2 × 10 -5 . For each cooling run we store different configurations above and below T g (see Fig. 2.1).

We are particularly interested in the "equilibrium" shear modulus G eq (for a precise definition see Sect. 2.4) as a function of temperature T for the imposed pressure P equal to zero. But the sampling procedure for G eq obliges us to work in the NVT ensemble (see Sect. 2.4.2 ). Therefore we made additional "equilibration" steps before sampling (see Fig. 2.1). Starting from the saved configuration at temperature T i from a cooling run with rate Γ = 2 × 10 -5 , we perform a volume relaxation in the NPT ensemble over a time t relax NPT = 0.75 × 10 5 . After the relaxation we measure the mean volume (t V = 0.25 × 10 5 ), and then quickly deform (t def = 10 3 ) our system to this mean volume and fix it. This typically creates some additional stresses in the system, so we also add an "equilibration" in an NVT ensemble (t eq NVT = 10 5 ) to relax these unwanted stresses. Finally we sample the quantities of interest in the NVT ensemble with the average pressure P = 0 over a maximum time t max = 10 5 . The times of the "equilibration" steps are tuned in such a way that the aging effects are weak on times of the data recording. All times here are presented in τ LJ . The procedure is similar to the one used in [START_REF] Schnell | Simulated glassforming polymer melts: Glass transition temperature and elastic constants of the glassy state[END_REF][START_REF] Mangalara | The relationship between dynamic and pseudo-thermodynamic measures of the glass transition temperature in nanostructured materials[END_REF].

Film

For the films we focused on the one chain length N = 16. Initial configurations are prepared from the bulk configuration at T = 0.7 by removing periodicity in z direction which mimics the free standing film. Fixing the area A of the film (lateral system size of The simulation procedure for the films is essentially the same as for the bulk; yet there are a few differences. In contrast to the bulk we cool the system while keeping the area of the film fixed. We cannot do the constant pressure (isotropically) cooling for the films, because of the presence of the surface tension which cannot be compensated in the liquid state that leads to the permanent decrease of the film area.

L = L x = L y = 23.
Since the system can only adapt the volume in z direction because we keep A fixed, it creates stresses in the film, which are parallel to the surface. To confirm this statement, we use the Irving-Kirkwood to define the perpendicular (or "normal", P N ) and parallel (or "tangential", P T ) components of the local pressure in the system [START_REF] Irving | The statistical mechanical theory of transport processes. iv. the equations of hydrodynamics[END_REF][START_REF] Varnik | Molecular dynamics results on the pressure tensor of polymer films[END_REF]:

P N (z) = ρ(z)k B T - 1 A a<b ∂U ∂r ab (z ab ) 2 r ab × 1 | z ab | Θ z -z a z ab Θ z b -z z ab P T (z) = ρ(z)k B T - 1 2A a<b ∂U ∂r ab (x ab ) 2 + (y ab ) 2 r ab × 1 | z ab | Θ z -z a z ab Θ z b -z z ab (2.21)
where ρ(z) denotes the density at z averaged over tangential coordinates x and y, Θ is the Heaviside step function. where P bulk T (z) is the "bulky" pressure inside the film, P γ T (z) -the surface pressure which is directly connected with surface tension γ:

2γ = ∞ -∞ (P N (z) -P γ T (z)) dz P N ≡ 0 === - ∞ -∞ P γ T (z)dz liq. = -P T h (2.23)
In the liquid the particles can easily rearrange inside the films, relaxing any stress inside the film (P bulk (z) ≡ 0). Thus in the liquid the surface tension can be determined from the total tangential pressure p T (second part of Eq. (2.23)), which is defined as:

P T = 1 h ∞ -∞ P T (z)dz (2.24)
In the glassy phase the particles are known to be trapped in instantaneous configurations depending on the cooling protocol, hence they are struggling to relax the stress inside the film (P bulk (z) < 0 because the film wants to shrink, see red curve in Fig. 2

.2).

Actually P bulk (z) in the glass depends purely on the preparation of the film. Keeping this in mind, we use two different equilibration procedures after cooling ("film preparations").

First, we equilibrate in the NVT ensemble keeping the area of the film fixed. Second, we use a short NPT (P = 0) run to relax the area of the film and then proceed with an NVT equilibration. Both preparations give the same pressure profile in the liquid. As for the glassy phase, additional NPT run creates P bulk (z) > 0 which tries to compensate P γ (z) (see green curve in Fig. 2.2). Showing the temperature dependence of the total tangential and normal pressures (multiplied by film thickness h for demonstration reasons), Figure 2.3 underlines the difference between two preparations: the same normal pressure equal to zero for all temperatures, the same tangential pressure in the liquid phase, different tangential pressures in the glassy state. Comparing these two preparations we will show the influence of the tangential pressure on the properties of the film.

Sampling viscoelastic properties

In this section we answer the following questions: What is the main quantity that we are interested in? Why is it so good? How can we measure it?

Equilibrium shear modulus and explicit deformation method

The equilibrium shear modulus G eq is an important order parameter characterizing the transition from liquid, G eq = 0, to the solid state, G eq > 0 [START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary conditions, temperature, and sampling time[END_REF][START_REF] Wittmer | Shear-stress relaxation and ensemble transformation of shear-stress autocorrelation functions[END_REF]. The shear modulus of an isotropic solid body may be determined from the long time limit of the stress increment δτ (t) = τ (t) -τ (0 -) due to an imposed small step strain γ (see. Fig. 2.4).

G eq = lim t→∞ G(t), (2.25)
where

G(t) = δτ (t)
γ and the instantaneous shear stress τ (t) at time t can be measured experimentally from the forces acting on the walls of the shear cell. At time t = 0 + the imposed strain γ provokes a uniform displacement of all particles. In general, such an affine displacement creates a shear stress τ (0 + ) with corresponding shear modulus

µ a = G(0 + ) = τ (0 + ) -τ (0 -) /γ.
With time the particles can rearrange in order to decrease the free energy. This leads to a further, non-affine, contribution to the shear modulus which we call µ f (where "f" stands for "fluctuation"):

G eq = µ a -µ f , (2.26)
This method is straightforward but the problem is obvious -you can not directly measure inherent G eq in a non-deformed state (γ = 0).

Certainly, the subscript "eq" in the name of G eq means that it must be used for the equilibrium systems, stable at all times. Although glasses are not equilibrium systems, they eventually flow in the long time limit (t → ∞), we will keep this denomination to characterize a long-time, in term of our simulation time, response of the system. Hence, we will additionally discuss the influence of the sampling time on G eq (see Sect. 4.1.4).

Stress fluctuation formalism

In contrast to the described method, where the stress fluctuations were induced by external deformation (γ > 0), one may also obtain G eq from the spontaneous fluctuations of the stress δτ 2 = (τ -τ ) 2 in the non-deformed cell (γ = 0). This is possible because of the linear response. Assuming a system at imposed particle number N , constant volume V and mean temperature T , the "stress fluctuation formalism" [START_REF] Squire | Isothermal elastic constants for argon: Theory and monte carlo calculations[END_REF][START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary conditions, temperature, and sampling time[END_REF][START_REF] Wittmer | Shear-stress relaxation and ensemble transformation of shear-stress autocorrelation functions[END_REF]70] defines the shear modulus G eq as a difference of the shear stress fluctuations in the conjugated ensemble at constant mean shear stress τ = τ (NVτ T-ensemble) and imposed shear strain γ (NVγT-ensemble):

G eq = βV δτ 2 | τ -βV δτ 2 | γ = µ a -µ f (2.27)
here β = 1/k B T is the inverse temperature with k B being Boltzmann's constant. In the limit of small deformations (|γ| ≪ 1) the affine term µ a = βV δτ 2 | τ can be rewritten as "simple average" [START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary conditions, temperature, and sampling time[END_REF][START_REF] Schnell | Simulated glassforming polymer melts: Glass transition temperature and elastic constants of the glassy state[END_REF][START_REF] Wittmer | Shear-stress relaxation and ensemble transformation of shear-stress autocorrelation functions[END_REF][START_REF] Wittmer | Fluctuation-dissipation relation between shear stress relaxation modulus and shear stress autocorrelation function revisited[END_REF] characterizing the second-order energy change:

µ a = μa = 1 V Ĥ′′ (γ) γ=0 (2.28)
with Ĥ being the full Hamiltonian of the system, Ĥ′′ (γ) the second derivative of the Hamiltonian with respect to the shear strain γ. The term "simple average" means that it can be computed in any ensemble assuming that the same state point is sampled.

For pair potentials:

µ a = µ B -P ex + P id (2.

29)

with P id = k B T ρ being the ideal pressure contribution, P ex = P -P id the excess pressure and µ B is a so-called "Born-Lamé coefficient" [START_REF] Born | Dynamical Theory of Crystal Lattices[END_REF]:

µ B = μB = 1 V l ∂ 2 U ∂r 2 l - 1 r l ∂U ∂r l (r α l ) 2 (r β l ) 2 r 2 l (2.30)
and the shear stress τ :

τ = τ = 1 V l ∂U ∂r l r α l r β l r l - N a=1 mv α a v β a (2.31)
where a labels the particles, l labels the interaction U (r l ) between the particles i and j with i < j, r l is the distance between corresponding particles, v a is the velocity of particle a (= 1, . . . , N ) with mass m, α and β define the plane in which the shear stress τ and the shear modulus G eq are measured, and • is the "simple average". In an isotropic body all three planes (xy, xz, yz) are equivalent, therefore to increase accuracy we define G eq as an average over these planes.

For isotropic d-dimensional system (d = 3 in our case) it is common practice to rewrite µ B in terms of the radial pair distribution function g(r) [START_REF] Zwanzig | Highfrequency elastic moduli of simple fluids[END_REF]:

µ B = 2π d(d + 2) ρ 2 ∞ 0 drg(r) d dr r 4 dU dr (2.32)
where the d-dependent prefactor appears from the averaging over angles. This form also emphases that of the instantaneous elastic response (µ a ) depend on structure of the material (g(r)) and the intraction potential (U ).

It is worth to mention that the ideal (kinetic) contribution to affine elasticity µ id a = P id are identical to the ideal part of the fluctuation term [START_REF] Wittmer | Fluctuation-dissipation relation between shear stress relaxation modulus and shear stress autocorrelation function revisited[END_REF]. Therefore all ideal contributions to the shear modulus G eq (Eq. (2.26)) cancel.

µ id f = δτ 2 id | γ with τid = -1 V N a=1 mv α a v β a
Since µ a is independent of the ensemble in which it is determined (Eqs. 2.28 to 2.30) and µ f is measuring the stress fluctuation in constant γ ensemble (µ f = βV δτ 2 | γ ), it is convenient to obtain the shear modulus G eq performing simulations in NVγT ensemble only. By taking into account that γ = 0, the ensemble is equivalent to the classical canonical ensemble (see "NVT sampling" part of Panel (a) on Fig. 2.1).

Displacement correlations

Another finite-T method for computing G eq , which is also of experimental relevance, is based on the connection of G eq with displacement correlations of the particles [START_REF] Klix | Elastic properties of 2d amorphous solids[END_REF][START_REF] Klix | Discontinuous shear modulus determines the glass transition temperature[END_REF].

Here we briefly summarize the theoretical approach. Starting with the displacements u i (t) = r i (t) -ri for each particle i of the system, where the average in time position of the particle (r i ) is taken as reference, the displacement field is defined as:

u(r, t) = i u i (t)δ(r -r i (t)) (2.33)
and its Fourier transform:

u(q, t) = 1 √ N N i=1 u i (t) exp(iqr i ).
(2.34)

The wavevector q must be commensurate to the square simulation box of linear length

L (q x = 2π
Lx n x with n x = 0, ±1, ±2, . . . ). The component of u(q) perpendicular to q corresponds to the transverse component u ⊥ (q) of the Fourier transformed displacement field. Using that according to the equipartition theorem every independent elastic mode corresponds to an average kinetic or potential energy k B T /2, continuum mechanics implies that [15, 73]

q 2 |u ⊥ (q)| 2 → k B T ρ/G eq for q → 0 (2.35)
The beauty of this method is that it can be implemented in the real experiments by recording the trajectories of all the particles for a long time, without knowing the interaction potentials in the system. But the general disadvantage of this method both for the computer simulation and the experiment is quite obvious: at each time t we must store all positions of the particles (N = 12288 for each component in our case). This is a huge amount of data comparing to the "stress fluctuation" method, were we store six quantities (3 for the µ a and 3 shear stress components τ ) for each time step. Nevertheless we will show in Section 4.1.3 that both methods give the same values of G eq .

G eq from the mean square displacement of the particles

Yet another method is to appeal to ideal MCT, where we may expect that the collective property G eq is related to a single-monomer quantity, the (so-called) Lindemann localization length r s (T ) by [START_REF] Schnell | Simulated glassforming polymer melts: Glass transition temperature and elastic constants of the glassy state[END_REF][START_REF] Fuchs | Aspects of the dynamics of colloidal suspensions: Further results of the mode-coupling theory of structural relaxation[END_REF]:

G eq (T ) ρk B T ∝ σ 2 6r 2 s (T ) (2.36)
with σ being the monomer diameter. The Lindemann localization length measures the average displacement of a monomer around its equilibrium position in the glass. The Lindemann criterion of melting says [START_REF] Götze | Complex Dynamics of Glass-Forming Liquids: A Mode-Coupling Theory[END_REF]: the polymer glass should melt if the displacement attains 10% of the monomer diameter. Following [START_REF] Schnell | Simulated glassforming polymer melts: Glass transition temperature and elastic constants of the glassy state[END_REF] we defined r s (T ) from the the meansquare displacement g 0 (t, T ) of a monomer through the equation g 0 (t s , T ) = 6r 2 s (T ). Here t s is determined at T g as a time where the typical displacement of a monomer is about 0.1σ and is kept constant for all other T . In Section 4.1.3 we compare the results of this method with the results of two methods presented above. One of the advantages of this method is that it gives an easy way to determine a layer resolved G eq (z) in the film (see Sect. 4.3.2).

Chapter 3

Glass transition temperature

Bulk properties

We begin this chapter with a discussion of the glass transition temperature T g in the bulk, which will serve as an important reference point for the analysis in the films. First we introduce the method to determine T g , which is commonly used in experiment. Further we investigated the influence of the chain length N and cooling rate Γ on T g .

T g determination and its chain length dependence

One of the classical experimental methods to determine T g is dilatometry [START_REF] Sperling | Introduction to physical polymer science[END_REF]. In this method T g is determined as the point in the volume versus temperature curve where the change of the slope occurs. This method is based on the assumption that the fluid and the glassy branch have the thermal expansion coefficient α = ∂ ln 1/ρ ∂T which is independent of the temperature far from the transition. Technically it is applied by finding an intersection of the straight-line extrapolations of the liquid and glassy parts of the volume-temperature curve. One can also find T g from the total energy versus temperature curve E tot (T ). Since volume is an extensive variable, i.e. it depends on the number of particles in the system, it is more convenient to use the number density ρ = N/V , which gives a possibility to easily compare our results to other studies.

The dependence of the density ρ on the temperature is illustrated in Fig. 3.1(a) where N = 4 and Γ = 2 × 10 -5 . It is seen that ρ has a bend on cooling, where the thermal expansion coefficient α jumps from the high liquid to a small glassy state value. Using the above method for determining T g we performed the simulations for the systems with different chain lengths N . The dependence of T g on N for two different cooling rates is illustrated in Fig. 3.1(b). As we can see T g increases with the chain length, first strongly for small N , then weakly approaching an asymptotic value T ∞ g . The solid lines in Fig. 3.1(b) correspond to a fit to the equation:

T g (N ) = T ∞ g - C g N , (3.1)
where T ∞ g = 0.408 is the glass transition temperature for an infinitely long chain and C g = 0.0833 is a constant (exact values correspond to cooling rate Γ = 2 × 10 -5 ). This chain length dependence was first explained by Fox and Flory in [START_REF] Fox | [END_REF]76] based on freevolume arguments. Our data confirms earlier works [START_REF] Schnell | Simulated glassforming polymer melts: Glass transition temperature and elastic constants of the glassy state[END_REF][START_REF] Frey | Simulated glass-forming polymer melts: Dynamic scattering functions, chain length effects, and mode-coupling theory analysis[END_REF]. For further details see [START_REF] Baschnagel | Polymer Glasses, ch. Glass Transition and Relaxation Behavior of Supercooled Polymer Melts: An Introduction to Modeling Approaches by Molecular Dynamics Simulations[END_REF]. Equation (3.1) works for all cooling rates Γ studied (see Fig. this representation the slope is independent of cooling rate Γ. It means that most of the cooling rate dependence comes from T ∞ g . Thus we can slightly rewrite an Equation (3.1) in the form

T g (Γ, N ) = T ∞ g (Γ) 1 - K N , (3.2)
where now K is the constant independent on cooling rate.

Cooling rate dependence of T g

As were mentioned before T g depends on the rate with which the liquid is cooled into the glass (see [START_REF] Angell | Relaxation in glassforming liquids and amorphous solids[END_REF][START_REF] Baschnagel | Polymer Glasses, ch. Glass Transition and Relaxation Behavior of Supercooled Polymer Melts: An Introduction to Modeling Approaches by Molecular Dynamics Simulations[END_REF][START_REF] Cavagna | Supercooled liquids for pedestrians[END_REF]. This dependence can be rationalized as follows. Let us assume that, starting from some initial temperature T i ≫ T g , we continuously cool the liquid with the rate Γ,

T (t) = T i -Γt ⇐⇒ t(T ) = T i -T Γ (3.3)
As T decreases, the time scale t(T ) associated with the cooling protocol increases, but not as quickly as the molecular relaxation time τ α (T ) (the "α relaxation time") of the liquid. In the supercooled liquid τ α (T ) increases in a super-Arrhenius fashion [START_REF] Berthier | Theoretical perspective on the glass transition and amorphous materials[END_REF] which can be fitted by the Vogel-Fulcher-Tammann (VFT) equation:

τ α (T ) = τ ∞ exp B T -T 0 (3.4)
were τ ∞ is the asymptotic relaxation time at high temperature (formally, lim T →∞ τ α (T ) = τ ∞ ), B is a material-characteristic temperature scale, and T 0 is the "Vogel-Fulcher temperature" at which the relaxation time appears to diverge [START_REF] Donth | The Glass Transition[END_REF][START_REF] Berthier | Theoretical perspective on the glass transition and amorphous materials[END_REF][START_REF] Cavagna | Supercooled liquids for pedestrians[END_REF]. As T traverses the region near T g , τ α (T ) becomes comparable to t(T ). For T < T g , τ α (T ) quickly exceeds t(T ) by many orders of magnitude so that the liquid freezes in a (nonequilibrium) amorphous solid on the time scales accessible to experiment or simulation. Operationally, we can therefore define T g by the criterion: τ α (T g ) = t(T g ). With Eqs. 3.3 and 3.4 we then get (up to corrections of order ln[1 + 1/lnΓ]):

T g (Γ) = T 0 g - B ln(AΓ) (3.5)
where

A (= τ ∞ /[T i T 0 ]
) is a characteristic cooling rate and T 0 g (= T 0 ) is the glass transition temperature for infinitely slow cooling. With this line of reasoning we expect that T g has a weak logarithmic dependence on cooling rate, as a consequence of the super-Arrhenius increase of the α relaxation time. In applications of Equation (3.5) to experimental [START_REF] Brüning | Glass transition on long time scales[END_REF] or simulation data [79-81] T 0 g , A and B are treated as adjustable parameter. Excluding the N -dependence, Figure 3.2 shows the T ∞ g as a function of the logarithm of the cooling rate Γ. In accord with many experimental results [START_REF] Brüning | Glass transition on long time scales[END_REF] we see that Equation (3.5) can describe the simulation data. Nevertheless, in the range of cooling rates studied the change of T ∞ g is rather weak and may also be approximated by a linear extrapolation in terms of log Γ (this behavior is proposed by Narayanaswamy in [START_REF] Narayanaswamy | Thermoreological Simplicity in the Glass Transition[END_REF]):

T ∞ g (Γ) = T ∞ g (Γ 1 ) + B log(Γ/Γ 1 ), B = dT ∞ g d log Γ Γ 1 (3.6)
where Γ 1 (= 2 × 10 -6 ) is some reference rate from which the extrapolation is started, B is the constant (= 0.832 ). Using the conversion to SI units (see Sect. 2.2) gives the constant B ≈ 4 K, which is a reasonable value when compared to experiment [START_REF] Mckenna | Comprehensive Polymer Science[END_REF]. Of course the typical cooling rates employed in current simulations (10 -9 -10 -12 K/s) are orders of magnitude faster than in experiment (typical experimental values are 10 -3 -100 K/s).

Despite this huge difference, the discussion of this section shows that the variations of T g with chain length or cooling rate qualitatively agree between experiment and simulation.

This may appear paradoxical, but is (likely) related to the fact that both laboratory and computer studies probe the emergence of nonequilibrium phenomena when the monomer relaxation time becomes comparable to the "experimental" time scale set by the respective cooling process. Similarly, as T is lowered below T g , the laboratory and numerical systems lose the ability to relax and freeze in a glassy solid.

Film properties

Now we move to the characterization of the glass transition in our free standing model polymer films. Here we focus on the film thickness dependence of the glass transition temperature T g for one chain length N = 16.

Density profile ρ(z)

The monomer density along the z axis is defined as:

ρ(z) = 1 A nN i=1 δ(z -z i ) (3.7)
where N n is the total number of monomers in the system and A is the area of the simulation box.

Technically ρ(z) is calculated from the probability P (z) to find the particle with position z:

ρ(z) = N nP (z) dV = N (z -dz 2 , z + dz 2 ) Adz (3.8)
where N (z -dz/2, z + dz/2) is the number of monomers in slab [z -dz 2 , z + dz 2 ], and dz is the binning.

From the MD data for ρ(z) we determine the position of the left (z L ) and the right (z R ) interface, the "plateau" density in the center of the film ρ 0 and the apparent width of the interface ω by fitting to the data the approximation [START_REF] Peter | Structure et relaxation structurale des fondus de polymères vitrifiables en couches minces[END_REF]:

ρ(z) = ρ 0 2 1 -erf z -z R ω ρ(z) = ρ 0 2 1 + erf z -z L ω (3.9)
The film thickness h is then defined as h = z R -z L . Since the interface is antisymmetric around z 0 , the plateau density ρ 0 is equivalent to the average density of the film ρ =

1 h ∞ -∞ dzρ(z).

T g determination in films

Similar to the bulk system (Sect. 3.1.1), we determine T g by monitoring the plateau density as the film is continuously cooled with a cooling rate Γ = 2 × 10 5 . Figure 3 ρ 0 (T ) gradually changes slope upon crossing T g . Following the procedure applied to the bulk sample (see Sec. 3.1.1) we can determine T g as the intersection point of linear extrapolations from the liquid and glass sides. This method is rather simple because it has four explicit parameters for the two linear regions. However, there are additionally two implicit parameters in the choice made for the beginning and end of the transition region which are not included in the fit. As mentioned in [START_REF] Peter | Structure et relaxation structurale des fondus de polymères vitrifiables en couches minces[END_REF] the result of this procedure depends on this interval chosen for the fit. Thus, these implicit parameters create a problem to study systematically the shifts in T g when the shifts become relatively big (for small shifts these implicit parameters could be kept constant).

A less ambiguous definition of T g which does not hinge on hidden fit parameters, was proposed in [START_REF] Dalnoki-Veress | Molecular weight dependence of reductions in the glass transition temperature of thin, freely standing polymer films[END_REF]. This work used the empirical expression Eq. (3.10) for the film thickness h(T ), which is valid for any measurement that varies linearly with temperature in the melt and glass region. Thus, we can apply it to ρ 0 (T ):

ρ 0 (T ) = w M -G 2 ln cosh T -T g w + (T -T g ) M + G 2 + ρ g (3.10)
This expression has only five parameters that have a clear meaning: T g -the position of the transition, ρ g -the value of the film plateau density ρ 0 at T = T g , w -the width of the transition between melt and glass, and M and G are the dρ 0 /dT slope values for the melt and the glass, respectively. Equation (3.10) perfectly fits our data for ρ 0 (T ) (see red dashed line on Fig. 3.4) and gives T g values which are consistent with values obtained by the standard two-line method. However, since it is free of the choices made for the fit interval, we apply this procedure to determine T g (h).

Thickness dependence of T g : T g vs h

Solely from Figure 3.4 we can see that the glass transition temperature of the film is depressed with respect to the bulk value (T bulk g = 0.40189). Herminghaus et al. [START_REF] Herminghaus | The glass transition of thin polymer films: some questions, and a possible answer[END_REF] suggested a formula for the thickness dependence of T g :

T g = T bulk g 1 + h 0 h , (3.11)
where h = h g = h(T g ) and h 0 is a characteristic length scale and the only fitting parameter in this model, since the bulk glass transition T bulk g can be determined independently (as in our simulation). The underlying assumption of this model is that the relaxation in a film close to T g is due to a coupling of the viscoelastic bulk to capillary waves at the free surface.

Our data for T g (h) are well described by Equation (3.11) with h 0 = 1.69 (see black dashed line on the inset of the Fig. 3.5). This model determines the parameter h 0 as h 0 = γ/E, where γ is the surface tension at the free surface and E is Young's modulus of the film. For our films the surface tension γ ≈ 1.5 (see Figure 2.3), which suggests an elastic modulus of approximately E ≈ 1. This value is not unreasonable [START_REF] Peter | Structure et relaxation structurale des fondus de polymères vitrifiables en couches minces[END_REF]. Although Equation (3.11) fits our data well, we found that another widely employed empirical form (Eq. (3.12)), based on the assumption of the existence of a liquid-like layer at the free surface [START_REF] Mirigian | Theory of activated glassy relaxation, mobility gradients, surface diffusion, and vitrification in free standing thin films[END_REF][START_REF] Keddie | Size-dependent depression of the glass transition temperature in polymer films[END_REF], works slightly better for our data (see inset of Fig. 3.5)

T g = T bulk g 1 - h 0 h δ , (3.12)
with h 0 = 1.32, δ = 1.0. For the big film thickness h ≫ h 0 , Equation (3.12) coincides with Equation (3.11), the only difference appears for h ≈ h 0 where Equation (3.12) has much steeper decrease. For a deeper discussion about the differences of both models along with other existent results and models see [START_REF] Alcoutlabi | Effects of confinement on material behaviour at nanometer size scale[END_REF][START_REF] Baschnagel | Computer simulation of supercooled polymer melts in the bulk and in confined geometry[END_REF].

Position dependence of T g : T g vs z

One of the most common and logical explanations of the T g reduction in free-standing polymer films is related to the presence of the free surface, which possibly allows for an enhanced mobility of the monomers (see Inset of Fig. 4.20). T g (h) alone cannot confirm this explanation. The direct check of this assumption is to measure a layer resolved glass transition temperature T g (z).

Although it was stated that the dynamic glass transition appears simultaneously in the whole film [START_REF] Böddeker | Dynamics near free surfaces of molecular dynamics simulated Ni 0.5 Zr 0.5 metallic glass films[END_REF], a lot of studies show that it is possible to define a layer resolved "pseudothermodynamic" T g (z) [START_REF] Peter | Structure et relaxation structurale des fondus de polymères vitrifiables en couches minces[END_REF][START_REF] Baschnagel | Computer simulation of supercooled polymer melts in the bulk and in confined geometry[END_REF][START_REF] Mirigian | Theory of activated glassy relaxation, mobility gradients, surface diffusion, and vitrification in free standing thin films[END_REF]. Recording density profiles at different temperatures along the cooling run we use an equation similar to Equation (3.10) to define T g (z) (see Fig. as a fitting parameter. As we see the transition for the surface layer occurs at the lower temperature compared to the layer at the middle of the film.

To find an analytic form of T g (z) we make a few assumptions. First we assume, that T g (h) of the film can be written as the simple mean of T g (z):

T g (h) = 1 h z R z L dzT g (z) = 1 0 dz * T g (z * ) = T g (z) (3.13)
where z L and z R are the positions of the free interfaces. The second hypothesis is that the only effect of each surface is to decrease the glass transition by some value ∆T g (z)

from T bulk g . Since both interfaces are identical, the resulting T g (z) can be written as:

T g (z) = T bulk g -∆T g (z) -∆T g (h -z) (3.14)
Differentiating then Eq. (3.13) and using (3.11), (3.14) we obtain:

T g (z * ) = T bulk g 2    1 + 2h 0 hz * 1 + h 0 hz * 2 + 1 + 2h 0 h(1-z * ) 1 + h 0 h(1-z * ) 2    (3.15)
Equation (3.15) fits well our data on T g (z). The only concern is that the parameter h 0 from the fit, h 0 = 2.52, is much bigger than one from the fit of T g (h) (h 0 = 1.69).

This results in the significantly lower values of T g (h) calculated from Eq. (3.13) (see two solid lines on Figure 3.7). Nevertheless, some references suggest to use a weighting function in Eq. (3.13) instead of the linear arithmetic mean [START_REF] Mangalara | The relationship between dynamic and pseudo-thermodynamic measures of the glass transition temperature in nanostructured materials[END_REF]. However, we found that much simpler explanation works for us: by calculating the average T g of the film from the plateau density versus temperature curve we are averaging T g (z) not over the whole range z ∈ [0, h], but over a smaller region z ∈ [0 + ξ, h -ξ]. This is logical, since plateau density is a quantity that assembles information about the middle of the film, excluding some region (apparently of the size ξ) close to the surface.

Interestingly, the penetration depth of the surface effects for T g (z) is much bigger than the apparent width of the interface ω calculated from the density profiles (ω < 0.05h for all temperatures studied). This large scale of T g (z) reduction is consistent with findings of Binder et al about the growth of dynamic length scale when approaching T g [START_REF] Scheidler | Cooperative motion and growing length scales in supercooled confined liquids[END_REF]. Chapter 4

ρ vs T and h vs T for different film preparations

Mechanical properties of glass-forming polymer melts

This chapter is dedicated to the results on viscoelastic properties. As before, we start with the bulk system. Focusing on one chain length N = 4, we investigate the temperature dependence of the short (µ a ) and long time (G eq ) response functions. In this context we compare different methods and discuss the effect of the finite time of the simulation. We continue with the bulk system, but move to the chain length dependence of viscoelastic properties. In the end, we use the results on the bulk system to explain the dependence of viscoelastic properties on the thickness of the film. In MD simulation time averaging is often used instead of ensemble averaging. In equilibrium in the infinite time limit these two averages give the same result. But the simulation time is always finite. Therefore we must study the influence of the sampling time ∆t on the quantities of interest. Following [START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary conditions, temperature, and sampling time[END_REF] we use gliding average to define the sampling time dependent average of the quantity â:

a(∆t) = a(∆t, s) s = 1 n s (∆t) s ns(∆t)-1 s=0 a(∆t, s), (4.1)
a(∆t, s) = 1 ∆t t<s+∆t t=s â(t) (4.2) n s = t max -∆t + 1 (4.3)
where the average • s is taken over all shifts s of the time window ∆t within time series t max (see Fig. If â specifies a fluctuation of some property b (â = δ b2 ), then:

δ b2 (∆t, s) = 1 ∆t t<s+∆t t=s b(t) 2 - 1 ∆t t<s+∆t t=s b(t) 2 , (4.4)
for b = √ βV τ this defines µ f (∆t).

µ f (∆t) = 2 ∆t ∆t 0 dt 1 - t ∆t [C(0) -C(t)] , (4.5)
where C(t) = τ (t + s)τ (s) s is the stress autocorrelation function. For all temperatures µ f is slowly increasing with the sampling time ∆t from µ f = 0 for short times to their plateau value for long times. For the temperatures T much higher than T g (T = 0.5 on Figure 4.1) this plateau is equivalent to µ a , giving zero shear modulus (G eq (T > T g ) = 0) that corresponds to a liquid state. This is not the case for the temperatures below T g , where the plateau value is smaller than µ a , leading to the finite G eq and a solid-like behavior.

Some technical details concerning the choice of ∆t must be mentioned here. From Eq. (4.1) we can say that µ f (∆t) could have some dependence on ∆t tmax , since averaging is carried out over n s (∆t) = t max -∆t+1 different "independent" µ f (∆t, s). That is why our aspiration to take the biggest sampling time (∆t = t max ) must be taken with care since for this time the statistics is lost drastically (see Top curve on Fig. 4.1). This can also be seen from Figure 4.1 (bottom) where we present µ f (∆t) for two different t max using a fixed temperature T = 0.3 < T g . As we see for ∆t ≪ t max the plots match perfectly, showing the tendency of µ f to converge to some value in a long time limit. In the vicinity of t max we observe a change of trend that could be misinterpreted as a further increase of µ f towards µ a . But the increase of t max shows that this trend-change is nothing but an artifact due to the loss of statistics. The similar finite time effect is observed in another work [START_REF] Helfferich | Continuous-time random walk approach to supercooled liquids: I. different definitions of particle jumps and their consequences[END_REF] for the waiting-time distribution of monomer jumps in supercooled polymer melt. All this motivated us to fix the ratio tmax ∆t = 10 (with t max = 10 5 τ LJ ) for the definition of the stress fluctuation term µ f = µ f (∆t) both for the bulk and film simulations.

In the glassy state the time-averages sample only one fixed glass configuration, hence we additionally average over different realizations of the glassy state at fixed temperature We turn now to the characterization of the temperature dependence of the shear modulus G eq and its components. Figure 4.3 presents µ a (black circles and blue diamonds) and µ f (red squares) obtained for a system with N = 4 monomers per chain. Note that both definitions of µ a , from Eq. (2.30) and Eq. (2.32), perfectly match, which shows that our system is isotropic. As indicated by the green solid lines, the affine shear elasticity µ a (T ) decays linearly with T . As expected [START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary conditions, temperature, and sampling time[END_REF], µ a (T ) changes its slope near T = T g . Panel (b)

shows that in this temperature range µ a depends linearly on density ρ(T ). This means would mean that the system reacts on deformation by creating a stress in the direction of deformation and not oppose to it. We would rather expect µ a to vanish, when the density decreases to zero (which Eq. (2.32) also suggests). For checking this, we performed simulations for higher temperatures T = 0.9 -2. The similar scaling with density is also found for µ a of Lennard-Jones liquid along fluidsolid coexistence [START_REF] Heyes | Scaling of lennardjones liquid elastic moduli, viscoelasticity and other properties along fluidsolid coexistence[END_REF].

In the end, Figure 4.5 shows the temperature dependence of the equilibrium shear modulus G eq . For low temperatures G eq linearly decreases with temperature. A linear variation of the shear modulus with the temperature is in very good agreement with experiments and other simulations [START_REF] Klix | Discontinuous shear modulus determines the glass transition temperature[END_REF][START_REF] Nadal | Continuous model for the shear modulus as a function of pressure and temperature up to the melting point: Analysis and ultrasonic validation[END_REF]. Moreover we observe a clear step-like behavior, but still continuous, from the high-temperature liquid behavior with G eq = 0 to the lowtemperature solid behavior with G eq > 0, that gives the possibility to define T g solely from the elasticity. The exact shape of the transition region is a bit ambiguous because it is strongly dependent on the sampling time ∆t (see Sect. 4.1.4).

Comparison to other methods

In order to confirm our results on G eq (T ) we want to compare it to other methods. temperatures. Following Eq. (2.36) we then define G eq as G eq = 1.035ρk B T /g 0 (t s )). As we see the results from all methods are in good agreement with each other, which justifies the accuracy of the stress fluctuation method.

Finite time effects

Before moving to the chain length dependence of G eq (T ) we want to answer the question:

what would happen with G eq (T ) if we choose different sampling time ∆t in Sect. 4.1.1. with its standard deviation δG eq (T ) (Inset) for a broad range of sampling times ∆t. Here δG eq (T ) characterizes the fluctuations between 100 different configurations. We emphasize that G eq (T ) decreases continuously without a jump singularity for all ∆t sampled.

In contrast to this, the standard deviation δG eq (T ) is strongly non-monotonous with a remarkable peak near the glass transition temperataure T g . The liquid-solid crossover of G eq (T ) at T g becomes systematically sharper with increasing ∆t. Note that the data for the two largest sampling times ∆t are similar but not yet identical. Thus it is difficult The transition becomes more sudden with increasing ∆t. Inset: Corresponding standard deviation δG eq showing a peak at T ≈ T g which becomes sharper with increasing ∆t.

to speculate on the behaviour for much larger sampling times and a jump discontinuity for ∆t → ∞ cannot be excluded. While the liquid-solid transition characterized by the ensemble-averaged shear modulus G eq (T ) becomes increasingly better defined, the peak of δG eq (T ) gets more singular. Note that for ∆t = 10 5 the peak value, δG eq (T ≈ 0.35) ≈ 3.5

is about a third of the sudden drop of the modulus G eq (T ) between T = 0.34 and T = 0.38.

The liquid-solid transition is thus masked by very strong fluctuatuions between the configurations.

Distribution of G eq

The striking peak of δG eq below T g seen in Fig. 4.7 begs for a more detailed characterization of the distribution p(G eq ; T, ∆t) of the time-averaged shear modulus G eq . Focusing on our largest sampling time ∆t = 10 5 , the main panel of Fig. 4.8 presents normalized histograms obtained using 3 × m = 300 measurements. We emphasize that the histograms are unimodal for all T and ∆t. The T -dependence of G eq and ∆G eq below T g seen in Fig. 4.7 is thus not due to, e.g., the superposition of two configuration populations representing either solid states with finite G eq and liquid states with G eq ≈ 0. The maximum Inset: Difference G med eq -G eq of the median G med eq and the ensemble average G eq vs. T for two sampling times. The difference has a peak slightly below T g corresponding to very lopsided distributions.

G max eq of the (unimodal) distribution systematically shifts to higher values below T g , in agreement with its first moment G eq (Fig. 4.7), while the distributions become systematically broader and more lopsided, i.e. liquid-like configurations with small G eq remain relevant. The increase of ∆G eq with sampling time ∆t seen in the inset of Fig. 4.7 is due to the broadening of p(G eq ) caused by the growing weight of small-G eq configurations (not shown). For even smaller temperatures T ≪ T g , the distributions get again more focused around their maxima G max eq (as expected from the inset in Fig. 4.7) and less lopsided.

That the large standard deviations and the asymmetry of the distributions are related may be demonstrated by comparing the first moment G eq of the distribution, its median G med eq and its maximum G max eq . One confirms that 0 < G med eq -G eq < G max eq -G eq below T g for all ∆t. As seen from the inset of Fig. 4.8, G med eq -G eq has a peak similar to ∆G eq becoming also sharper with increasing ∆t.

Bulk: Chain length dependence

To study the chain length dependence we used only 10 independent cooling runs for every N ∈ {4, 8, 16, 32, 64}. This gives already reasonable errors in G eq (Fig. . Anyway this scaling not only shows the universality of G eq (T ) behavior, but also allows to compare the results of our simulation to the "real world experiments" (see Sect. 4.4).

Film: Thickness dependence

In the previous section we found a scaling that helps to merge all bulk data of G eq (T )

for different N on one master curve. Since this scaling eliminates the N dependence, we continue with fixed N = 16 for studying the influence of reduced dimension, or in other words the film thickness dependence.

µ a , µ f vs T

We start with presenting results on the affine elastic constant µ a as a function of temperature T for four different film thicknesses (see Fig. 4.13). As for the bulk, µ a changes slope around T g . In the liquid phase we see that µ a increases with film thickness h, approaching bulk value µ bulk a for big h. Panel (b) on the Figure 4.13 shows µ a vs 1/h for T = 0.5 (here the bulk system correspond to a point 1/h = 0). The behavior is clearly linear:

µ a (h, T > T g ) = µ bulk a (T ) 1 - h ω h (4.11)
with h ω being a fitting parameter, which is found to be slightly dependent on temperature.

For our example, T = 0.5, h ω ≈ 0.1607.

Surprisingly this h-dependence almost disappears in the glass state, where we observe µ a ≈ µ bulk a for all temperatures below T g .

The bulk result µ bulk a (ρ, T, N ) = ρk B T + µ B (ρ, N ) that we obtained in the previous section, suggests that µ a for a given N and T must be equal for systems with equal density and different for the systems with different densities. Panel (c) of Fig. 4.13 shows the plateau density ρ 0 of the films. As we see, we get a completely opposite behavior:

in the liquid state we have equal ρ 0 for all films but have a strong h-dependence of µ a , whereas in the glass ρ 0 varies with h but µ a does not. This could call into question the validity of the results of the previous section, when one applies them to the films.

However, the reason for these controversial result between film and bulk could be hidden in the inhomogeneity of the film density and the dependence of µ B density. Hence it is necessary to check whether the functional dependence of µ a (ρ), obtained for the bulk system, does also only hold for local properties in the film.

Following [START_REF] Van Workum | Local elastic constants in thin films of an fcc crystal[END_REF] we defined a local, or layer resolved, µ a (z):

µ B (z) = 1 A a<b ∂ 2 U ∂r 2 ab - 1 r ab ∂U ∂r ab (r α ab ) 2 (r β ab ) 2 r 2 ab 1 | z ab | Θ z -z a z ab × Θ z b -z z ab µ a (z) = µ B (z) -P T (z) + 2ρ(z)k B T (4.12)
where Θ is the Heaviside step function, P T (z) is the local tangential pressure defined by Eq. (2.21). Note that this expression does not depend on the volume of the system or the binning size. It relates the local stress at z to a homogeneous strain. dependence could be applied in the film, this gives µ film a (z) = µ bulk a (ρ(z)) (for all T and every N ). Indeed, Fig. 4.14 shows that µ bulk a (ρ(z)) (black solid lines) matches exactly the µ a (z) computed directly from Equation (4.12), which justifies our assumption and the results for a bulk system. In the inset of Fig. 4.14 we also show that the average of µ a (z) coincide with the results of the direct measurements of µ a for all temperatures T :

µ a = 1 h ∞ -∞ dzµ a (z) = 1 h ∞ -∞ dzµ bulk a (ρ(z)) (4.13)
Dividing the integral into "interface" (with a strong z-dependence of µ a ) and "plateau"

(where µ a (z) = const = µ bulk a (ρ 0 )) parts we get here 2ω stands for the width of the interface (see Fig. 3.3). If µ bulk a was linear in ρ, then 2ω dzµ a (z) = ωµ a (ρ 0 ), so there would be no h-dependence. Nevertheless,the dependence is not simply linear. By taking into account Eqs. (4.9 -4.10), which we rewrite for simplicity as µ bulk a (ρ) = aρ 3 + bρ, we get:

µ a = µ bulk a (ρ 0 ) 1 - 1 h 2ω -2 2ω
µ a (h) = µ bulk a (ρ 0 ) 1 - 2ω h aρ 3 0 (1 -2C 3 ) + bρ 0 (1 -2C 1 ) µ bulk a (ρ 0 ) C i = 1 2ω 2ω dz ρ(z) ρ 0 i , i = 1, 3 (4.15)
where C i depends on the shape of the film interface in the density profile. In our case the density is well approximated by Eq. (3.7), which gives C 3 ≈ 0.23625 and C 1 = 1 2 . Equation (4.15) explains the decrease of µ a with h found in the liquid state (Eq. (4.11)).

For T = 0.5 Equation (4.15) suggests h ω ≈ 0.16 ± 0.02 (ω ≈ 0.57, µ a (ρ 0 ) ≈ 89.32, a = μ∞ LJ = 22.237, ρ 0 = 1.00549), which is in good agreement with the value obtained from the fit of Eq. (4.11) (≈ 0.1607). Nevertheless the derivation suggests that we must have the h-dependence for all temperatures T , both in the liquid and the glass. The reason why we do not see it for T < T g , is that decrease of µ a with decreasing h is masked

The effect of the film preparation

The results for G eq for the film with h(T g ) = 8.6 are shown on the Figure 4.17. As wee see G eq for both preparations has the same behavior in the glass transition region, which is consistent with the results presented in Sect. 3.2.5 where we showed that T g is identical for the two preparations. Since the main difference of the two preparations is that they create different tangential pressures in the film (see Sect. 3.2.5), this motivated us to assume that tangential pressure is responsible for h-dependence of G eq .

Bulk: G eq vs P in the glass It is actually not a surprise that the equilibrium shear modulus depends on the pressure of the system (see [START_REF] Nadal | Continuous model for the shear modulus as a function of pressure and temperature up to the melting point: Analysis and ultrasonic validation[END_REF] for example). But in order to confirm our assumption we move back to the corresponding bulk system (N = 16). For the temperature T = 0.2 deep in the glass, we performed simulations with different pressures P ∈ [1, -1]. The same as before, we used the sampling time ∆t = 10 4 with the data acquisition time t max = 10 5 .

The number of independent cooling runs is equal to 10. The maximum strain in the system corresponds to P = -1 and equal to [L(P = -1) -L(P = 0)] /L(P = 0) = 0.0054. G eq (T ) for the film with h(T g ) = 8.50. Briefly, the first preparation fixes the area A of the film for all temperatures, while the second gives the film some time to relax the area (see Sect. 2.3.2 for more details). Black solid line correspond to the bulk system. where A ≈ 0.11 is constant.

G eq scaling with P T and T g (h)

Taking into account Equation (4.16) we tried to remove the influence of the pressure in our simulation. In order to do this, we plot the rescaled shear modulus G * eq = Geq 1+AP T as a function of unitless temperature T * = T /T g (see Fig. 4.19). Interestingly all four curves for G eq (T ) collapse onto one master curve. Moreover in the glass phase this master curve is equivalent to the bulk one. This means that the h-dependence of G eq (T ) in the glass is nothing else but a manifestation of the pressure dependence and is the consequence of the chosen film preparation in which P T ∝ 1/h. Nevertheless, the shape of the transition differs for the bulk and the films. The answer to this is that T g varies strongly with the position in the film, and is much suppressed close to the interface (see Sect. 3.2.4).

A glance at the layer resolved G eq

To get a feeling of why the shape of the transition differs, we look at the layer resolved G eq (T ) calculated with Eq. (2.36) (Sect. 2.4.4). Figure 4.20 depicts G eq (T ) for h(T g ) = 21.33 calculated at three positions in the film: at the interface, close to the interface and in the middle of the film. We used a modified mean square displacement g 0 (r) = 3 2 ∆x 2 + ∆y 2 , which was calculated in the slabs of thickness ∆z = 1. We have already showen that the decrease of G eq in the glassy stated is consequence of the negative tangential pressure in the film and Figure 4.20 one more time proves this. Since P T (z) is larger close to the interface (see Fig. 2.2 for a typical pressure profile), we also observe a huge drop in G eq (T ). In the middle of the film the pressure is still negative but way smaller, which gives a smaller decrease of G eq (T ). Interestingly, in the middle of the film the shape of the transition is equivalent to the bulk one. As for the interface, we see that the transition is shifted to lower temperatures. This shift in T g give rise to the less sudden transition, when G eq (z, T ) is averaged over the film. 4.12 ) and four films are compared to the literature data. Orange crosses show the recent experimental data on 2d colloidal glass [START_REF] Klix | Discontinuous shear modulus determines the glass transition temperature[END_REF]. Cyan pluses and blue solid line are simulation results on 3dKA and 2dpLJ models, respectively [START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary conditions, temperature, and sampling time[END_REF].

Chapter 5 Conclusion

Using a flexible bead-spring polymer model we performed molecular dynamics (MD) simulation to study the influence of confined geometry on the properties of glass-forming polymer melts. In this work we mainly focused on two important aspects of amorphous materials: the glass transition temperature and the mechanical properties. In both cases we start our analysis by looking at the bulk system and then move to the films.

Glass transition temperature

Starting with the bulk system, we measure the glass transition temperature T g via dilatometry, but using the density instead of the volume. We next confirmed the chain length dependence T g (N ) = T ∞ g 1-K N and showed that T ∞ g has a weak logarithmic dependence on cooling rate Γ [START_REF] Baschnagel | Polymer Glasses, ch. Glass Transition and Relaxation Behavior of Supercooled Polymer Melts: An Introduction to Modeling Approaches by Molecular Dynamics Simulations[END_REF]. By doing so, we obtain the well characterized "reference point" for our films.

Looking on the monomer density profile ρ(z) in the film we approximated it with erf function to obtain the film thickness h and the plateau density ρ 0 , which corresponds to the density in the middle of the film. Although T g in the films is usually determined from the temperature dependence of the film thickness h, we used the plateau density instead of h in order to be consistent with the bulk analysis. We showed that T g in the film is depressed with respect to the bulk, and the depression becomes more pronounced the power of the density-scaling of µ a .

• In extension of µ a analysis, it would be useful to split G eq into LJ and bond contributions, or, even more practical, to intra-and inter-chain contributions.

For the films:

• It would be interesting to check if the relationship µ a (z) = µ bulk a (ρ(z)) still holds for the simulated supported films and for the "real world" polymer films.

• The layer resolved G eq (z) must be calculated with stress fluctuation formalism. This will give a deeper understanding of the h-dependence. Moreover, it will allow to study more precisely the shape of G eq (T ) at the glass transition.

• Implement displacement field method for the films, since it is related to creep compliance. Implement then also biaxial deformation to have a direct connection with "nanobubble inflation" experiment [START_REF] O'connell | Novel nanobubble inflation method for determining the viscoelastic properties of ultrathin polymer films[END_REF]. The stress fluctuation formalism could be also applied to measure the bulk modulus of the system [START_REF] Schnell | Simulated glassforming polymer melts: Glass transition temperature and elastic constants of the glassy state[END_REF].

K eq = λ a + 2 3 G eq + P (5.1)

Similarly to Eq. (2.26) the λ is defined as

λ = λ a -λ f (5.2)
where again λ a specifies the instantaneous response, and λ f is characterizing the long time rearrangements of the particles in the system. Using the general equations from [START_REF] Schnell | Simulated glassforming polymer melts: Glass transition temperature and elastic constants of the glassy state[END_REF] give

λ a = µ B λ f = βV δ P 2 | γ (5.3)
where δ P 2 | γ is the spontaneous fluctuations of the stress in the NVγT ensemble. for the ∆t-dependence of the shear-stress fluctuations in agreement with the more direct demonstration given in Ref. [START_REF] Wittmer | Shear-stress relaxation and ensemble transformation of shear-stress autocorrelation functions[END_REF].

Figure 1 :

 1 Figure 1: Figure principal: Température de transition vitreuse T g en fonction de h.

Figure 2 :

 2 Figure 2: Partie (a): Module de cisaillement G eq fonction de la température T pour les systèmes de différentes longueurs de chaîne N . Partie (b): Module de cisaillement sansunités G eq /k B ρT g à l'équilibre en fonction de la température réduite T /T g , où T g est la température de la transition vitreuse calculée par dilatométrie.

Figure 3 :

 3 Figure 3: (a): Module de cisaillement à l'équilibre G eq en fonction de la température T pour le système isotrope et 4 films d'épaisseurs différentes (N = 16). (b): Module de cisaillement à l'équilibre renormalisé G * eq = G eq /(1 + AP T ), où P T est la pression tangentielle et A = 0.11, en fonction de la température adimensionnée T /T g . (c): Pression tangentielle P T = (P xx + P yy )/2 en fonction de la température pour tous les films.

Figure 1 . 1 :

 11 Figure 1.1: Idealized variations in the heat capacity C p , the volume V , the relaxation time τ and the equilibrium shear modulus G eq along the cooling through the glass transition T g . For the volume V , the difference between the glass transition and crystallization is stressed.

Fig. 1 . 1 )

 11 Fig. 1.1). This increase of the relaxation time is accompanied by the emergence of shear rigidity and viscoelastic behavior. In the other words, the equilibrium shear modulus G eq becomes finite in the glass (on the typical time scales), in contrast to a zero value in the liquid (see Panel (d) of Fig. 1.1).

Farago, 2 Model, Methods, Simulation Details 2 . 1

 221 and J. P. Wittmer In Connie B. Roth (Ed.), Polymer Glasses (pp. 55-105) by CRC Press.Most of the results in Section 4.1 will be published in the following article:• 2017 Shear modulus and shear-stress fluctuations in polymer glasses. Kriuchevskyi, I., Wittmer, J. P., Meyer H. and Baschnagel, J. In preparationChapter Molecular dynamics simulationMolecular dynamics (MD) simulations[START_REF] Frenkel | Understanding Molecular Simulation[END_REF][START_REF] Allen | Computer Simulation of Liquids[END_REF] consider a classical system of N particles with masses m i = m (i = 1, . . . , N ) that interact via effective forces, in our case determined by the coarse grained model potentials defined below (Sect. 2.2). At any time t the positions (r) and velocities (v) of all the particles fully specify the microscopic configuration of the system. Starting from initial configuration ({r i }, {v i }, i = 1, . . . , N ), MD obtains the time evolution of this state by solving Newton's equations of motion:

  26, 58-62]. Virnau et al. studied the phase separation kinetics of a mixture of hexadecane (C 16 H 34 ) and carbon dioxide (CO 2 ) [60, 61]. By matching the critical point of the liquidgas transition in hexadecane with that of beadspring chains (consisting of 5 monomers each) they obtain σ ≃ 4.5 × 10 -10 m and ǫ ≃ 5.8 × 10 -21 J (≃ 420 K). Assuming σ = 4 × 10 -10 m Paul and Smith converted τ LJ to seconds by comparing the late-time diffusive dynamics of chemically realistic models for nonentangled melts of polyethylene

  performed MD simulations of polymer systems consisting of n linear chains, each chain being composed of N monomers. More precisely, the studied values are N = 4 (n = 3072), N = 8 (n = 1536), N = 16 (n = 768), N = 32 (n = 384), and N = 64 (n = 192). The simulation box is cubic with periodic boundary conditions in all spatial directions.

Figure 2 . 1 :

 21 Figure 2.1: Panel (a): The scheme of the simulation procedure. The time ratios are not realistic and made just for the purpose of explanation (see Sec. 2.3 for exact values). The blue line shows the cooling in which the temperature decreases linearly with time according to Eq. (2.20). Blue squares are saved configurations. For each of these saved configurations there is a constant temperature run. It consists of: an NPT (P = 0) run, needed to relax the volume of the sample and measure the mean volume that corresponds to the pressure equal to zero; quick deformation to this mean volume; NVT run to relax residual stresses created by the deformation; NVT run used to sample the quantities of interest. Panel (b): The schematic behavior of the volume in the constant temperature run. After relaxation the mean volume is measured in an NPT run (red dashed line). The red circle shows an instantaneous volume of the sample which is generally not equal to the mean volume. From this value the volume is linearly changed to the mean value and fixed (green line). It is made by simultaneously scaling all three linear dimensions of the box. Panel (c): The schematic pressure behavior in the constant temperature run. First the pressure relaxes to the target value (red dashed line) and than oscillate around.

Figure 2 .

 2 Figure 2.2 shows the profiles P T (z) and P N (z) in the film. The normal pressure is zero both in the liquid and in the glass. As for the tangential pressure, in general it can be divided into two parts:

Figure 2 . 2 :

 22 Figure 2.2: Tangential and normal pressure profiles of the film.

Figure 2 . 3 :

 23 Figure 2.3: Integrals of the tangential and normal pressure profiles as a function of temperature T . Using Equation (2.23) gives the surface tension γ ≈ 1.5 close to the liquid-glass transition.

Figure 2 . 4 :

 24 Figure 2.4: Panel (a): The scheme of the simple shear experiment with γ being the shear strain, red arrow the affine displacement and green arrow additional non-affine displacement. Panel (b): Sketch of the Eq. (2.25) and Eq. (2.26) with G(t) being the shear relaxation modulus (black solid line), G eq the equilibrium shear modulus (blue solid horizontal line), µ a the affine shear elasticity (red dashed horizontal line), µ f the non-affine shear elasticity.

Figure 3 . 1 :

 31 Figure 3.1: Panel (a): Density ρ (= N n V ) versus temperature. The determination of the glass transition temperature T g as the intersection of linear extrapolations (solid lines) from the glass and liquid branches of ρ is indicated. The result was obtained with the cooling rate Γ = 2 × 10 -5 and the chain length N = 4 (circles) and N = 64 (squares). Panel (b): T g versus N . T g obtained from the continuous cooling run via dilatometry at pressure P = 0. The symbols correspond to different cooling rates Γ. The solid lines show a fit to Eq. (3.1).

  3.1(b) for two examples of Γ = 2 × 10 -5 and Γ = 2 × 10 -6 ). The coefficients T ∞ g and C g depend on the cooling rate.

Figure 3 .Figure 3 . 2 :

 332 Figure 3.2: Panel (a): Scaling plot of T g versus 1/N . All data follow the same linear master curve which is consistent with Eq. (3.1). Panel (b): T ∞ g versus logarithm of the cooling rate Γ. The symbols correspond to T ∞ g from the approximation of the T g vs N curves with Equation (3.1). Solid line is the linear fit (Eq. (3.6)) and the dashed line is a 1/ ln Γ fit based on the VFT equation (Eq. (3.5)).

Figure 3 . 3 :

 33 Figure 3.3: Density profile ρ(z) for a free-standing film at T = 0.50 (red circles). The film consists of n = 768 chains with chain length N = 16. The blue solid line is the approximation of the "right" interface with Eq. (3.9). The vertical green line shows the position of the interface z R . Two orange dashed lines (z R -ω and z R + ω) indicate the interface. The green dotted line correspond to the plateau density of the film ρ 0 . The black solid line is the bulk density ρ bulk at the same p and T , which is quite close to ρ 0 for this temperature.

Figure 3 .

 3 Figure 3.3 shows the equilibrated density profile of a film at T = 0.5 and the approximation of the one of the interfaces with Eq. (3.9). The vertical green line indicates the position of the interface. At this temperature the plateau density of the film ρ 0 is very close to the bulk density ρ bulk at the same p and T .

  .4 provides an example. It depicts the temperature dependence of ρ 0 for two films with a different starting thickness.

Figure 3 . 4 :

 34 Figure 3.4: The plateau density ρ 0 as a function of temperature T for two films with a different starting thicknesses. Red dashed lines are fits to the Eq. (3.10) with h g = 21.33, T g = 0.3745 and h g = 6.61, T g = 0.3258.

Figure 3 . 5 :

 35 Figure 3.5: Main panel: T g (h)/T bulk g

  3.6).

Figure 3 .

 3 Figure 3.7 shows T g (z * = z/h), indeed T g (z * ) is reduced close to the interface resulting in the smaller T g of the film comparing to the bulk.

Figure 3 . 6 :

 36 Figure 3.6: Main panel: Density profile for T = 0.5. Schematic view of the mapping from position z to rescaled positions z * = (z -z L (T ))/h(T ), where z L is the position of the "left" interface, z R -the position of the "right" interface. The mapping is made for every temperature. Here dz is the thickness of the layer in which the density is averaged. Inset: Layer resolved density as a function of temperature: blues circles correspond to the layer at the surface of the film, green squares -the layer in the middle of the film. Black solid lines correspond to a fit to equation similar to Equation (3.10), which gives T g (z * )

Figure 3 .

 3 Figure 3.8 shows the difference in plateau density ρ 0 and thickness h of the film as a function of temperature for two film preparations (see Sect. 2.3.2). The effect of the area relaxation can be summarized as follows: the area A decreases, film thickness h increases, and the plateau density ρ 0 stay the same in the liquid, but increases in the glass. As before, we determine T g from Eq.(3.10). The exact values are T g = 0.3745 and T A g = 0.3732 for the preparation with non-relaxed and relaxed area of the film, respectively. The difference between both values is negligible. Thus, we can conclude that our different preparations do not affect T g .

Figure 3 . 7 :Figure 3 . 8 :

 3738 Figure 3.7: Layer resolved glass transition temperature T g (z * ) for the film with h g = 21.33, where z * is the rescaled position in the film z * = z/h (see Figure 3.6; red diamonds: dz ≈ 0.08h; blue circles: dz ≈ 0.3h). Dashed line -bulk glass transition temperature T bulk g ,

4. 1 4 4. 1 . 1

 1411 Stress fluctuation formalism applied to bulk polymer melts: Case study N = Sampling time dependence: µ f , µ a , G eq

2 . 1 )

 21 , and n s (∆t) -the number of shifts for a fixed window size. Examples of such quantities are τ and μa (Eq. (2.31) and Eq. (2.29)).

Figure 4 .

 4 Figure 4.1 shows the sampling time dependence of the affine term µ a (∆t) and nonaffine term µ f (∆t). Being a simple mean the expectation value of the affine term µ a does not depend explicitly on the sampling time ∆t. In contrast µ f (∆t) is known to have strong dependence on ∆t (see Appendix C):

  ∆t)/n s (0) Eq.(1.1)

Figure 4 . 1 :

 41 Figure 4.1: Bottom: Affine µ a and non-affine µ f parts of the shear modulus G eq as a function of the sampling time ∆t. The chain length N = 4. The black circles and the red line corresponds to µ f (∆t) and µ a (∆t) respectively for T = 0.5 > T g . Green triangles and blue line corresponds to µ f (∆t) and µ a (∆t) for T = 0.3 < T g . For the above data t max = 10 5 . Magenta diamonds correspond to µ f (∆t) at T = 0.3 calculated over the smaller time series t max = 10 4 . Top: The number of shifts n s (Eq. (4.3)) as a function of sampling time (see Equation (4.1)). Dashed lines show the different choice for the sampling time ∆t. The effect of the choice on the temperature dependence of G eq is shown on Figure 4.7. Our choice of the sampling time ∆t = ∆t 3 = 10 4 with t max = 10 5 .

Figure 4 . 2 :

 42 Figure 4.2: The equilibrium shear modulus G eq as a function of sampling time ∆t for different T (T g (N = 4) ≈ 0.38). The vertical dashed lines correspond to different choices of the sampling time, which are the same as in Fig. 4.1 and Fig. 4.7.

Figure 4 . 3 :

 43 Figure 4.3: Panel (a): Affine µ a and non-affine µ f parts of the equilibrium shear modulus G eq as a function of temperature T (chain length N = 4). Black circles: µ B part of µ a calculated with Eq. (2.30); blue diamonds: µ B calculated from radial pair distribution function g(r) (Eq. (2.32)). Panel (b): µ a (T ) plotted versus density ρ(T ). In this range the dependence is clearly linear µ a = Aρ + B, with A = 117.85 and B = -40.458 (red solid line).

8 .Figure 4 . 5 :

 845 Figure 4.5: The temperature dependence of the equilibrium shear modulus G eq (T ) (chain length N = 4). Black dashed line correspond to the glass transition temperature T g defined from dilatometry (see Sect. 3.1.1). Solid green line indicates zero.

Figure 4 .Figure 4 . 6 :

 446 Figure 4.6 depicts the equilibrium shear modulus G eq from the previous section, calculated with the stress fluctuation formalism, along with the results from two other methods: the displacement correlations and the mean square displacement one (see Sect. 2.4 for a short review of these methods). The choice of t s for the later method is illustrated in Inset of Fig. 4.6, where we show the mean square displacement of the monomers (g 0 (t)) for different

Figure 4 .

 4 Figure 4.7 answers this question. It displays the value of G eq (T ) (Main panel) together

Figure 4 . 7 :

 47 Figure 4.7: Finite time effects. Orange dashed line specify T g calculated from dilatometry (see Sect. 3.1.1). Main Panel: Shear modulus G eq (T ) for different samping times ∆t.

5 Figure 4 . 8 :

 548 Figure 4.8: Main panel: Distribution p(G eq ) for ∆t = 10 5 for different temperatures as indicated. Inset: Difference G medeq -G eq of the median G med eq and the ensemble average G eq vs. T for two sampling times. The difference has a peak slightly below T g corresponding

Figure 4 . 10 :

 410 Figure 4.10: Panel (b): µ LJ B /ρ 3 as a function of temperature T . In this range of temperatures µ LJ B /ρ 3 with a good approximation is independent of T , which gives µ LJ B = Aρ 3 . Panel (a): Average in temperature µ LJ B /ρ 3 = A as a function of 1/N . Black solid line correspond to a linear fit: A(N ) = µ ∞ LJ 1 + 1 αN , with µ ∞ LJ = 22.237 and α = 2.02.

Figure 4 . 11 :

 411 Figure 4.11: Panel (a): µ f as a function of temperature T for different chain lengths N . Black solid lines correspond to data for µ a (T, N ). Panel (b): Average µ f of the glassy part versus 1/N .

Figure 4 . 12 :Fig. 4 .

 4124 Figure 4.12: Panel (a): Equilibrium shear modulus G eq as a function of temperature T for the systems with different chain length N . Panel (b): Unitless equilibrium shear modulus G eq /k B ρT g as a function of unitless temperature T /T g , where T g is the glass transition temperature calculated from dilatometry (see Sec. 3.1.1).

Figure 4 . 13 :

 413 Figure 4.13: Panel (a): µ a as a function of temperature T for films with four different thicknesses and the bulk system. Panel(b): µ a as function of inverse of the film thickness (1/h) for T = 0.5. Red solid line correspond to a linear fit µ bulk a (T ) 1 -hω h (Eq. (4.11)), with µ bulk a

Figure 4 .

 4 Figure 4.14 depicts µ a (z) calculated from Eq. (4.12) for T = 0.55 for two films with different thickness (h(T g ) = 8.50 and h(T g ) = 21.33). In the middle of the film µ a (z) is equivalent to the bulk µ a for both (and actually for all four) films. This agrees with our expectation since the plateau density, which is the density in the middle of the film, is equal in the liquid state for all films. Even more, if we now suppose that the µ bulk a (ρ, T, N )

Figure 4 . 14 :

 414 Figure 4.14: Main panel: The local µ a (z) calculated from Eq. (4.12) for two films with different thickness (green circles: h(T g ) = 8.50 and blue triangles: h(T g ) = 21.33) at temperature T = 0.55. Black solid lines correspond to µ a (z) = µ bulk a (ρ(z)). Red solid line corresponds to the bulk value of µ a at T = 0.55. Inset: µ a (T ) from the Figure 4.13 for the h(T g ) = 8.50 (green circles) is compared to the average of the local µ a (z) over the film (black solid line).

Figure 4 . 16 :

 416 Figure 4.16: The equilibrium shear modulus G eq as a function of temperature T for four films of different thickness (symbols) along with the bulk system (black solid line).

Figure 4 .

 4 Figure 4.18 shows that G eq increases linearly with pressure P :

Figure 4 . 17 :

 417 Figure 4.17: The influence of the film preparation on the equilibrium shear modulus

Figure 4 . 18 :Figure 4 . 19 :

 418419 Figure 4.18: The equilibrium shear modulus G eq of the bulk system as function of pressure P in the glass. The temperature T = 0.2, N = 16. Black solid line is the fit to Equation (4.16), with A ≈ 0.11.

Figure 4 . 20 :

 420 Figure 4.20: Main panel: Equilibrium shear modulus G eq as a function of temperature for three positions in the film (h(T g ) = 21.33): at the interface (black squares), close to the interface (green circles) and in the middle of the film (blue triangles). The red solid line corresponds to overall G eq (T ), same as in Fig. 4.16. Inset: Mean square displacement of the particles (g 0 (t)) for the three positions in the film, specified in the main panel, at T = 0.3.

  g )=21.33 h(T g )=8.50 h(T g )=6.61 h(T g )=4.82

Figure 4 . 21 :

 421 Figure 4.21: Unitless equilibrium shear modulus G eq /k B ρT g as a function of unitless temperature T /T g . The results for the bulk (black solid line corresponds to the data for all chain lengths N = 4, 8, 16, 32, 64, see Fig. 4.12 ) and four films are compared to the

  Our results on the glass transition temperature confirm essentially what is already known from previous simulation studies of flexible bead-spring model. The analysis of the mechanical properties in the whole temperature range is the new key result of this work.

Figure 5 . 2 :

 52 Figure 5.2: The temperature dependence of the equilibrium shear modulus G eq (T ) and the equilibrium bulk modulus K eq for the chain length N = 4. The average is performed over 100 independent coolings. Black dashed line corresponds to the glass transition temperature T g defined from dilatometry (see Sect. 3.1.1).

Figure 5 . 2 N - 1 s=0(

 521 Figure 5.2 shows the equilibrium shear modulus G eq (T ) and the equilibrium bulk modulus K eq for the chain length N = 4 averaged over 100 independent cooling runs.

) contributions:

l -1 r l ∂U LJ ∂r l

(r α l ) 2 (r β l ) 2 r 2 l (4. [START_REF] Angell | Relaxation in glassforming liquids and amorphous solids[END_REF] where the sums are taken over all bonded and non-bonded pairs of particles respectively.

As we see on Figure 4.4, µ bond B > µ LJ B for all densities(temperatures), meaning that the short time elastic response is dominated by bond interaction. The bond contribution is linear in density, and all non linearity of µ B comes from its LJ contribution µ LJ B , which can be fitted to Aρ 3 +Bρ. To summarize, we show here the set of equations that describes the temperature dependence of µ a : where A,B,C do not depend on the temperature (but could still have some chain length dependence, which will be discussed in Sect. 4.2.1) and pressure P = 0 in our simulation. both in the glass and liquid states (see Appendix A, where we show the direct impact of this reduction of the statistic).

4.2.1 ρ, µ a , µ f vs T As we saw in Figure 3.1 the density ρ depends not only on temperature T , but also increases with chain length N for a fixed temperature. This N -dependence is the consequence of the different distances between bonded and non-bonded monomers (see [START_REF] Frey | Propriétés viscoélastiques des fondus de polymères vitrifiables[END_REF]), which in our case are given by the equilibrium bond length (l 0 = 0.967) and the minimum of Lennard-Jones potential (r min = 2 1/6 σ LJ ) respectively. Hence the volume occupied by the end monomer V end is larger than the volume of the inner monomer V in , and the volume of the chain can be written as

Inserting this into the definition of density we get:

where we remind that M is the total number of monomers in the system and n = M/N is the number of chains.

In order to understand the chain length dependence of µ a = 2ρk ). These arguments suggest the form: temperature, but increasing with N . More precisely A is found to be linear in 1/N (see Panel (a) of Fig. 4.10 ), which gives:

where constants μ∞ LJ (= 22.237) correspond to µ LJ B /ρ 3 for infinite chains at density ρ = 1 and the physical origin of α(= 2.02) is not yet understood.

This N -dependence of the LJ term can be explained in the same way as for µ bond B , but taking into account that the contribution to the sum depends on the neighbor shells. For first neighbors there are

terms in the sum of Equation (4.6), where n 1 -is the number of first neighbors. This qualitatively explains the N -dependence of µ LJ B in Equation (4.10) and give the connection of the fitting parameter α to the physical quantity n 1 : α = n 1 2 -1. We used the word "qualitatively" because the value of α from the fit equal to 2.02 suggests n 1 ≈ 6, which is reasonable but still much less than the real value ≈ 13 (calculated from the integral of the radial distribution function g(r)).

Finally, Figure 4.11 depicts the µ f (T ) for different chain lengths N . For all N high temperature µ f is identical to µ a . The low-temperature glassy µ f is independent of T , but has a strong N dependence. The glassy plateau of the µ f is found to be linear in 1/N (see Fig. As for the bulk, the high temperature µ f is equivalent to µ a . Interestingly, in the glass phase µ f changes its h-dependence. For all temperatures below T g it decreases with increasing film thickness h. We will try to explain this h-dependence in the context of G eq below.

G eq vs T : h-dependence?

Now we move on to our key result on the h-dependence of the equilibrium shear modulus G eq . Using again Eq. (2.26) we obtain G eq for our four films. As wee see from Fig. 4.16 for all films, the liquid part G eq is equal to zero. In the glass phase G eq decreases with decreasing h strongly, which is the consequence of the h-dependence of µ f , since µ a does not depend on h for all temperatures below T g . On the other hand the slope of G eq (T ) is fully determined by µ a (T ).

In order to understand the origin of this h-dependence of G eq (T ) (and at the same time that of µ f (T )) in the glass, we compare the results for two film preparations (see Sect. 2.3.2).

Comparison to other models and experiment

This last section is dedicated to our aspiration to find a universality in the viscoelastic behavior of different systems. Here we compare our results on G eq for the bulk and the film systems to different simulated models and recent experiment from the literature. We consider two models: first is the three-dimensional Kob-Andersen model (3dKA) [START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary conditions, temperature, and sampling time[END_REF],

second is the two-dimensional polydisperse Lennard-Jones model (2dpLJ) [START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary conditions, temperature, and sampling time[END_REF] (we thank Joachim Wittmer and Hong Xu for providing the raw data). The experiment is on twodimensional colloidal binary mixture (we thank Peter Keim for sending us a raw data from his article [START_REF] Klix | Discontinuous shear modulus determines the glass transition temperature[END_REF]). G eq was calculated with the stress fluctuation formalism for both simulated models and with the displacement correlation method in the experiment.

We plot the unitless shear modulus Geq = Geq(T )/ρk B T T /Tg as a function of unitless temperature T = T /T g , as in Sect. 4.2.2, which gives the possibility to directly compare our results to the literature. Despite the differences in nature of these glass-forming liquids and in spatial dimension there is a remarkable similarity in the behavior. We want to stress a few points here. First, in the liquid state all models show Geq ≡ 0, except the experimental system. This nonzero value in experiment is due to two reasons [START_REF] Klix | Private communication[END_REF]: a) in the fluid, it is hard to stabilise the system; b) finite time effects occur dominantly on the fluid site, longer sampling times would reduce G eq (at least, shorter sampling times increase G eq ). Second, the Tdependence of G eq near T g is found to be similar for all systems. Of course, the exact form of the transition depends on the procedure and particularly on times used in the simulation/experiment. But if we would think in terms of G eq jump at T g , we observe that this jump is almost identical for all systems. Third, in the glass phase there is a clear separation between 2d and 3d systems. Deep in the glass (T → 0) our bulk system and 3dKA model converge to the same value of Geq . In the same time 2dpLJ model and 2d

experimental data are identical for all T < T g , and converge to a higher Geq for T → 0.

Interestingly, our films show the intermediate behavior, going from Geq closer to 3d value for big h to Geq equal to 2d value for the smallest h studied.

with decreasing h. The thickness dependence of T g can be well described by T g (h) = T bulk g /(1 + h 0 /h) [START_REF] Herminghaus | The glass transition of thin polymer films: some questions, and a possible answer[END_REF]. Recording the monomer density profiles ρ(z) along the cooling, we also obtain the layer resolved "pseudo-thermodynamic" T g (z). In agreement with recent experiments and simulations [START_REF] Peter | Structure et relaxation structurale des fondus de polymères vitrifiables en couches minces[END_REF][START_REF] Baschnagel | Computer simulation of supercooled polymer melts in the bulk and in confined geometry[END_REF][START_REF] Mirigian | Theory of activated glassy relaxation, mobility gradients, surface diffusion, and vitrification in free standing thin films[END_REF] we find that T g (z) decreases strongly with decreasing distance to the surface. This proves that the reduction of T g with decreasing h

in free-standing polymer films is related to the presence of the free surface.

Mechanical properties

Using the "stress fluctuation formalism" [START_REF] Squire | Isothermal elastic constants for argon: Theory and monte carlo calculations[END_REF] we focused on two different response functions:

the short-time (or infinite frequency) shear modulus µ a and the long time ("equilibrium") shear modulus G eq = µ a -µ f , where µ f measures the fluctuations of the stress.

Temperature dependence

We started our analysis from the case study of N = 4 in the bulk. We found that the sampling time effects are not important for the µ a but become crucial when the fluctuation term µ f is considered. We made a choice of the sampling time ∆t = t max /10, with the data acquisition time t max = 10 5 τ LJ , and fixed it for all simulations. The ensemble average was performed over m = 100 independent cooling runs.

Moving to the temperature dependence, we obtained that the short time response function, µ a (T ), can be well described in terms of the density as µ a (T ) = Aρ 3 (T )+Bρ(T ).

The LJ potential is found to be responsible for the ρ 3 part, whereas the part proportional to ρ is mostly comes from the bond potential. A similar scaling with density is also found for elastic moduli of Lennard-Jones liquid along fluid-solid coexistence [START_REF] Heyes | Scaling of lennardjones liquid elastic moduli, viscoelasticity and other properties along fluidsolid coexistence[END_REF].

As for the long time response, the equilibrium shear modulus G eq is found to be an unambiguous measure of the transition from the liquid. It shows a clear step-like behavior, from the liquid state with G eq = 0 to a solid state with G eq > 0. The exact shape of the transition depends on the sampling time. We stressed that the transition becomes sharper when the sampling time increases, but stays continuous. In addition, we showed that it is also possible to determine G eq from the monomer mean square displacement or the displacement correlations [START_REF] Klix | Elastic properties of 2d amorphous solids[END_REF][START_REF] Klix | Discontinuous shear modulus determines the glass transition temperature[END_REF][START_REF] Schnell | Simulated glassforming polymer melts: Glass transition temperature and elastic constants of the glassy state[END_REF][START_REF] Fuchs | Aspects of the dynamics of colloidal suspensions: Further results of the mode-coupling theory of structural relaxation[END_REF].

Chain length dependence

We performed simulations for five different chain lengths N = 4, 8, 16, 32, 64 in the bulk.

The same sampling time ∆t = 10 4 was used, but the ensemble average was taken over m = 10 independent cooling runs since it gives a reasonable results for a smaller computational cost (see Appendix A). We found that µ a increases with N . As before we split the total µ a into contributions from bond and LJ potentials. We found that the bond part is strongly increases with N as 1-1/N , whereas the LJ part is weakly decreases with N as 1+1/αN , with α ≈ 2.02. While the increase in the bond part, from N = 4 to N = 64, is ≈ 23% of the total µ a (N = 4), the drop in the LJ part barely reaches 3%. We managed to explain this N -dependence based on the simple argument of "counting" the terms in the sum of Eq. (2.30) (Sect. 4.2.1). In the end we presented G eq for different N . In the glass state it has a weak increase with N . We managed to collapse the data for all five chain lengths by plotting the unitless equilibrium shear modulus Geq = G eq (T )/ρk B T g as a function of unitless temperature T = T /T g .

Film thickness dependence

Moving to the films, we started our analysis by looking at µ a (T ) for four different film thickness h(T g ) = 21.33, 8.50, 6.61, 4.82 for N = 16. In agreement with [START_REF] Forrest | Brillouin light scattering studies of the mechanical properties of thin freely standing polystyrene films[END_REF] we found that this infinite frequency shear modulus does not depend on h in the glassy state.

Nevertheless, in the liquid state it increases with film thickness h. Measuring the local µ a (z), we showed that this increase is connected to a drop of µ a (z) close to the interface.

We showed that µ a (z) can be fully described in terms of the bulk results as µ a (z) = µ bulk a (ρ(z)). Using µ bulk a (ρ) ∼ ρ 3 + ρ we derived the h-dependence of µ a in the film. We stressed that ρ 3 term is responsible for decrease of the µ a with decreasing h. We also stated that in the glassy state this decrease of µ a is masked by the increase of the plateau density (which is in contradiction to [START_REF] Forrest | Brillouin light scattering studies of the mechanical properties of thin freely standing polystyrene films[END_REF], where the density is claimed to be constant in h).

We showed the h-dependence of G eq (T ). It is found to decrease with decreasing h, which is also observed both in amorphous polymer films [START_REF] Stafford | Elastic moduli of ultrathin amorphous polymer films[END_REF] and in fcc crystal films [START_REF] Van Workum | Local elastic constants in thin films of an fcc crystal[END_REF].

Here we showed that this h-dependence is a manifestation of the pressure dependence of G eq (T ). Thus it is the consequence of the particular film preparation, which in our case (where we fix the area of the film along the whole simulation) creates an h-dependent tangential stresses P T in the film. To prove this, we went back to the bulk system and performed simulations for different pressures P at T = 0.2 < T g . We obtained the pressure dependence of G eq (P ), which was found to be linear: G eq (P ) = G eq (0)(1 + AP ), with A ≈ 0.11. By plotting G eq (T )/(1 + AP T ) vs T /T g for our four films, it is possible to collapse all data onto one master curve.

Comparison to the literature

In the end, we compared our results on G eq for the bulk and the film systems to other sytems from the literature. We considered two model systems, the three-dimensional Kob-Andersen model (3dKA) [START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary conditions, temperature, and sampling time[END_REF] and the two-dimensional polydisperse Lennard-Jones model (2dpLJ) [START_REF] Wittmer | Shear modulus of simulated glass-forming model systems: Effects of boundary conditions, temperature, and sampling time[END_REF], and one experimental system, the two-dimensional colloidal binary mixture [START_REF] Klix | Discontinuous shear modulus determines the glass transition temperature[END_REF]. We showed, that by plotting the unitless shear modulus Geq = Geq(T )/ρk B T T /Tg as a function of unitless temperature T = T /T g it is possible to compare these largely different systems. We obtained a remarkable similarity in the behavior of Geq . If one would think in terms of G eq jump at T g , we observed that this jump is almost identical for all systems. Nevertheless, the behavior in the glass state is found to be dependent on the dimensionality of the system.

Outlook

Although our work gives some insight into the influence of the free surface on the properties of the glass-forming polymer films, there still remain a lot of open questions. Further simulation study could address the following questions.

For the bulk system:

• In [START_REF] Heyes | Scaling of lennardjones liquid elastic moduli, viscoelasticity and other properties along fluidsolid coexistence[END_REF] it is suggested that the power of the density scaling is obtained from the scaling of the distances r = rρ -1/3 . Reference [START_REF] Heyes | Scaling of lennardjones liquid elastic moduli, viscoelasticity and other properties along fluidsolid coexistence[END_REF] shows that the repulsive part of the potential gives rise to a ρ 5 dependence and the atttractive to a ρ 3 dependence of µ a (for the 12-6 LJ potential). Since the stress fluctuation formalism gives the possibility to directly split µ a into any possible contribution, it is possible to split the LJ part into the "attractive" and "repulsive" parts. This will show the origin of the observed ρ 3 scaling of µ a , and also answer the question if the ρ 5 + ρ 3 scaling suggested in [START_REF] Heyes | Scaling of lennardjones liquid elastic moduli, viscoelasticity and other properties along fluidsolid coexistence[END_REF] is also a possible answer in our case". To extend this, it would be interesting to directly see the influence of the power of the potential (U (r) When we switched from N = 4 case to the study of chain length dependence in Chapter 4 we decreased the number of independent cooling runs from 100 to 10. We commented that this gives "already reasonable errors in G eq ". Here we want to directly show the influence of this decrease of the statistics. Figure 5.1 depicts G eq averaged over 100 and 10 independent cooling runs. As wee see there is not much difference between these two cases. this result is then seen (assuming a sufficiently slow shear-stress barostat) to generalize to C τ τ (t) = G(t) -λG eq with G(t) being the shear-stress relaxation modulus.

Shear-stress fluctuations in self-assembled transient elastic networks is systematically scanned from the liquid limit (∆x ≫ 1) to the solid limit (∆x ≪ 1)

where the network topology is quenched and an ensemble average over m independent configurations is required. Generalizing previous work on permanent networks it is shown that the shear-stress relaxation modulus G(t) may be efficiently determined for all ∆x using the simple-average expression G(t) = µ a -h(t) with µ a = G(0) characterizing the canonical-affine shear transformation of the system at t = 0 and h(t) the (rescaled)

Appendix C

Relation between C(t) and µ f (∆t)