<figure><img src="image_1.png" /></figure>

HAL
open science

# Elastic Moduli of Crystalline Thin Films: How to get the Intrinsic Coefficients from Molecular Simulations?

Patrice Porion, Joël Puibasset

<figure><img src="image_2.png" /></figure>

## To cite this version:

Patrice Porion, Joël Puibasset. Elastic Moduli of Crystalline Thin Films: How to get the Intrinsic Coefficients from Molecular Simulations?. 2025. hal-05251815

HAL Id: hal-05251815

https://hal.science/hal-05251815v1

Preprint submitted on 12 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Elastic Moduli of Crystalline Thin Films: How to get the Intrinsic Coefficients from Molecular Simulations?

Patrice Porion and Joël Puibasset

Université d'Orléans, CNRS, ICMN, UMR 7374, Orléans, France.

(*Electronic mail: joel.puibasset@cnrs-orleans.fr)

(*Electronic mail: patrice porion@cnrs-orleans.fr)

(Dated: 4 September 2025)

Author to whom correspondence should be addressed: joel.puibasset@cnrs-orleans.fr

Molecular simulation is a powerful tool to model nanometric solids, in particular to predict the elastic properties of nanometric rods or films with realistic potentials. However, when implementing simulations, several questions arise. For instance, in a mechanical test, where should the forces be applied? On the outermost atoms, or on the virtual Gibbs surface? Which deformation should be considered? That deduced from the displacement of the outermost atoms, or from the displacement of the Gibbs surface? As a matter of fact, the responses differ depending on whether one is more interested in theoretical considerations or practical implementation in molecular simulations. This leads to two definitions of the elastic constants, with significant differences when considering small systems. This work proposes a simple formalism, allowing anyone is interested in elastic properties to switch between these two points of view. A simple Lennard-Jones model of crystalline thin film is used to illustrate this formalism, with an application to the study of the impact of the film thickness on its elastic constants at zero temperature.
# I. INTRODUCTION

The importance of recently developed nanometric devices (nanoelectromechanical systems, NEMS) cannot be underestimated.¹ This generates many experimental researches on the elastic properties of nanometric systems such as cantilevers, as well as the development of accurate models.²⁻⁸ In particular, many studies have focused on the influence of the size of the nanometric devices on their elastic properties. For such systems, presenting a large surface-to-volume ratio, a large effect is generally attributed to the surface stress induced by the free surfaces.⁹⁻²⁰ Theoretical approaches have been developed for a long time, in order to better understand experimental data.²¹⁻²⁸ From a fundamental point of view, the nanometric devices are generally seen as a piece of bulk matter with interfaces, leading to a simple decomposition of the elastic properties into bulk plus surface contributions. One consequence is that one expects a linear dependence of the elastic constants with the inverse of a characteristic distance measuring the system size (for instance its thickness in the case of a solid film). This point of view remains however controversial.¹³,¹⁷,²⁹⁻³² In particular, one may underestimate the fact that for a thin film, the two opposite solid-vacuum interfaces can actually influence each other if the film thickness is similar to or smaller than the range of the interactions.

Molecular simulation is a powerful tool to study these nanometric systems and predict their properties, thanks to the possibility to use accurate interatomic potentials, as well as many algorithms developed in recent years.²⁹,³³⁻⁴⁰ In molecular models, the measurement of elastic constants generally follows an experimental route that consists in applying forces and measuring the strain-stress curves for small deformations. However, several questions arise. In particular, one may wonder where the forces should be applied: directly on the outermost atoms constituting the free surfaces, or on the virtual surface (the so-called Gibbs dividing surface) that delimits the volume of the solid in the standard decomposition into bulk + surface contributions? Another interesting point is whether it is more appropriate to introduce the elastic constants deduced either (i) from the local deformation measured in the center of the system, or (ii) from the global deformation deduced from the displacement of the outermost atoms, or (iii) from the displacement of the Gibbs surface? The difficulty arises from the fact that the response differs depending on whether one is interested in theoretical considerations or practical implementation in molecular simulations, leading to two possible sets of elastic constants. These questions are also relevant in situations where the solid is in contact with another one, such as in composite or multilayers.⁹,¹⁰,¹³,³⁰,⁴¹⁻⁵¹ Such
considerations are not clearly stated in the literature, and deserve special attention, in particular if one is interested in comparing molecular simulation results with theory or experiments. This work proposes a simple formalism allowing anyone is interested in elastic properties to switch between these different points of view. It is in particular shown that this formalism is a powerful tool to convert the elastic constants as given by the simulations into the corresponding intrinsic quantities as defined from a theoretical point of view, enabling a simple analysis in terms of bulk + surface contributions. It will reveal particularly illuminating to treat the problem of the dependence of the elastic constants of a thin film with its thickness: the relations greatly simplify when presented in terms of the intrinsic elastic constant, while the relationships involving the constants from the simulations are more complex.

The paper is organized as follows: we first present the thin film model and introduce the important definitions relevant to the film mechanics: strain and stress, with particular emphasis on the theoretical and simulation viewpoints. The general formalism linking both is then presented. Finally, the simulation results are given, and, after transformation, are analyzed in terms of bulk + surface contributions. A simple model emerges, relying on a few numerical quantities, allowing one to predict the mechanical behavior of the film with an excellent accuracy.

## II. MODEL AND METHODS

### A. Notations

This section presents a list of the main quantities and notations used in this paper, with short definitions, so that the reader can easily refer to them later.

Geometry of the bulk FCC crystal relaxed at zero temperature:

$a_b$: bulk lattice parameter,

$\lambda_b = a_b/2$: distance between the bulk (001) planes.

Geometry of the thin film cut in the FCC crystal:

$L_x, L_y$: lateral dimensions of the strained film,

$L$: common value of $L_x$ and $L_y$ for the unstrained film,

$N$: number of atoms in the film,

$n$: number of crystallographic (001) planes,

$H'_e$: distance between the outermost (001) planes,
$H_{\varepsilon}^{G}$: distance between the two Gibbs dividing surfaces,

$\lambda_{\varepsilon}$: distance between the two innermost (001) planes,

$\Delta_{r}$: outermost layer displacement (surface relaxation).

Elasticity (bulk):

$\boldsymbol{\sigma}_{\alpha\beta} \in \{x,y,z\}$: stress tensor (general),

$\boldsymbol{\varepsilon}_{\alpha\beta} \in \{x,y,z\}$: symmetrized strain tensor (general),

$\boldsymbol{\sigma}_{1} = \boldsymbol{\sigma}_{xx}$, etc., $\boldsymbol{\varepsilon}_{1} = \boldsymbol{\varepsilon}_{xx}$, etc.: Voigt notation,

$C_{ii}^{b}$: bulk diagonal stiffness constant,

$C_{i \neq j}^{b}$: bulk off-diagonal stiffness constant,

$E^{b}$: bulk Young's modulus,

$\nu^{b}$: bulk Poisson's ratio.

Elasticity (thin film):

$\boldsymbol{\varepsilon}_{1}$, $\boldsymbol{\varepsilon}_{2}$: in-plane strains (Voigt notation),

$\boldsymbol{\varepsilon}_{3} = (\boldsymbol{\lambda}_{\varepsilon} - \boldsymbol{\lambda}_{b})/\boldsymbol{\lambda}_{b}$: out-of-plane strain in the center,

$\boldsymbol{\varepsilon}'_{3} = (H'_{\varepsilon} - H'_{0})/H'_{0}$: global out-of-plane strain,

$\boldsymbol{\sigma}_{i}$, $\boldsymbol{\sigma}'_{i}$: strains conjugate to $\boldsymbol{\varepsilon}_{i}$ and $\boldsymbol{\varepsilon}'_{i}$,

$C_{ij}$, $C'_{ij}$: stiffness coefficients relative to $\boldsymbol{\varepsilon}_{i}$ and $\boldsymbol{\varepsilon}'_{i}$,

$E_{i}$, $E'_{i}$: Young's moduli relative to $\boldsymbol{\varepsilon}_{i}$ and $\boldsymbol{\varepsilon}'_{i}$,

$\nu_{ij}$, $\nu'_{ij}$: Poisson's ratios relative to $\boldsymbol{\varepsilon}_{i}$ and $\boldsymbol{\varepsilon}'_{i}$,

$s_{ij}$: surface stress tensor.

Thermodynamics:

$T$: temperature,

$U_{f}$: internal energy of the film,

$u_{b,\varepsilon}$: energy per atom for the strained bulk FCC,

$A_{f}$, $A_{b}$: film and bulk Helmholtz free energy,

$\Omega_{f}$, $\Omega_{b}$: film and bulk Grand potential,

$N^{ex}$, $n^{ex}$, $\Omega_{f}^{ex}$, $A_{f}^{ex}$: excess quantities,

$\gamma$ or $\gamma_{\varepsilon}$: surface free energy of the strained film,

$\mu$: chemical potential of the film,

$F_{x}$, $F_{y}$: lateral forces along x and y.

Interatomic potentials:

$\boldsymbol{\varepsilon}_{LJ}$, $\boldsymbol{\sigma}_{LJ}$: Lennard-Jones potential parameters,
$r_c$: interaction cutoff.

## B. Thin film model

This study focuses on thin films with crystalline FCC structure, built from a piece of bulk crystal without defect. The lattice parameter of the relaxed (zero pressure) bulk crystal, at zero temperature, is denoted $a_b$. In the following, the subscript “b” will be used to specify any quantity belonging to the bulk solid phase. The thin film is obtained by cutting a slab along the dense crystallographic (001) planes. By convention, the axes are chosen parallel to the main crystallographic directions of the crystal. The x and y axes are those parallel to the film surfaces, while the z direction is normal to the free surfaces (see Fig.1). The film is infinite in the x and y directions, but finite along z, and characterized by its number of crystallographic (001) planes $n$. By convention, $n$ will be used as a convenient parameter to define the film thickness.

Since the film is homogeneous in the x and y directions parallel to the free surfaces, it is possible, for theoretical analysis as well as for numerical calculations, to consider any finite portion of the film of lateral dimensions denoted $L_x$ and $L_y$ (see Fig.1). This system then contains a finite number of atoms $N$, proportional to $n$, that has to remain constant upon film deformations. Therefore, upon an homogeneous deformation $\varepsilon_{xx}$ in the x direction, the system size becomes $L_x(\varepsilon_{xx}) = (1 + \varepsilon_{xx})L_x(0)$, where $L_x(0)$ is the initial system size at zero deformation. The same applies in the y direction. By convention, in this work, the deformations are measured with respect to the bulk relaxed crystal. Therefore, $L_x(0)$ or $L_y(0)$ correspond to the initial dimensions of the slab as cut from the bulk, prior to any spontaneous or imposed deformation. The $x \leftrightarrow y$ symmetry allows one to choose $L_x(0) = L_y(0) = L$ without any loss of generality. However, the deformations in the x and y directions are independent and may differ.

In the z direction, the system is highly inhomogeneous. In the central region of the film, the local geometry and the mechanical properties are expected to be close to those of the bulk crystal, while they most probably significantly depart from the bulk while approaching the free surfaces. In particular, one observes that the crystallographic planes close to the free surfaces undergo relaxation of their position, materialized by a small displacement. This relaxation results from the imbalance in forces applying to the crystallographic planes close to the free surfaces, due to the removal of atoms when creating the film. We will see that the amplitude of this relaxation depends on the pressure that may apply on the film, or, equivalently, on the deformation in the core of the
film. On the other hand, for a film thick enough, the crystallographic planes in the central region
are equally spaced as in the bulk phase. The interplanar spacing in the center is denoted $\lambda_{\epsilon}$, where
the subscript “$\epsilon$” is a generic notation to point out that the system is deformed. $\lambda_{\epsilon}$ is expected
to have the same value as in the bulk phase submitted to the same strain. In particular, for the
reference film with $\epsilon_{xx} = \epsilon_{yy} = 0$ and free to relax along the z axis (implying $\epsilon_{zz} = 0$ in the center
of the film), $\lambda_{\epsilon=0} = \lambda_{b} = a_{b}/2$.

### C. Thermodynamics of surfaces

The presence of two solid-vacuum interfaces for the thin film affects its thermodynamic and
mechanical properties. According to the thin wall approximation, any extensive quantity calcu-
lated for the film can be decomposed into (i) a volume contribution from a piece of matter with
properties similar to the bulk ones, and (ii) a surface contribution corresponding to the so-called ex-
cess quantity. This decomposition relies on a conventional definition of the volume of the system,
in particular the position of the frontiers in the z direction where local quantities exhibit variations.
One possibility is to define the film thickness as the distance between the two outermost layers,
denoted $H'_{\epsilon}$ (see Fig. 2), where the subscript “$\epsilon$” is again used to remind that the film is possibly
deformed. This quantity can easily be measured or imposed during the course of the simulation
and can be seen as a practical definition of the film thickness for simulation purposes. Note that
$H'_{\epsilon}/\lambda_{\epsilon}$, where $\lambda_{\epsilon}$ is the interplanar spacing in the center of the film, would not give the number $n$
of crystallographic planes in the film. The difference is, by definition, the excess number of crys-
tallographic planes $n^{ex} = n - H'_{\epsilon}/\lambda_{\epsilon}$. From a theoretical point of view, a more natural definition of
the film thickness follows the usual prescription of two symmetric Gibbs equimolar dividing sur-
faces positioned so that the excess quantity of matter vanishes (i.e. $n^{ex} = 0$). The corresponding
film thickness, defined as the distance between the two symmetric Gibbs surfaces, is denoted $H_{\epsilon}^{G}$.
From $n^{ex} = 0$, one gets $H_{\epsilon}^{G} = n\lambda_{\epsilon}$ (see Fig. 2). An important point has to be emphasized. $H'_{\epsilon}$ does
depend on the amplitude of the surface relaxation, while $H_{\epsilon}^{G}$ and the corresponding position of the
Gibbs dividing surfaces are insensitive to it. However, both are affected if the film is deformed in
the z direction. It is also important to keep in mind that, even without relaxation, $H_{\epsilon}^{G}$ is slightly
larger than $H'_{\epsilon}$ by approximately one interplanar spacing.

Another important quantity is the surface tension $\gamma_{\epsilon}$, defined as the surface excess per unit area
of the grand potential of the film, $\gamma_{\epsilon} = \Omega_{f}^{ex}/(2L_{x}L_{y})$, where a factor two is introduced to take into
account the two faces of the film. The total grand potential of the film is $\Omega_f = A_f - \mu N$, where $A_f$ is its Helmholtz free energy, $\mu$ the chemical potential and $N$ the number of atoms in the film. By construction, $\gamma_\epsilon$ is independent of the position of the Gibbs dividing surface. At zero temperature, $A_f = U_f$, the internal energy of the film. Since $\mu$ is the variation of the energy when introducing one atom at constant surface area, this atom has to be introduced in the core of the film. This is equivalent to say that the chemical potential is that in the central bulk-like region of the film, i.e. $\mu = u_{b,\epsilon}$, the bulk crystal energy per atom (at the same strain as in the core of the film). Therefore,

$$ \gamma_{\epsilon} = \frac{U_{f}^{ex} - N^{ex}u_{b,\epsilon}}{2L_{x}L_{y}} = \frac{U_{f} - Nu_{b,\epsilon}}{2L_{x}L_{y}}, \qquad (1) $$

where the excess quantities $U_f^{ex}$ and $N^{ex}$ in the numerator can be replaced by $U_f$ and $N$ since the volume contributions cancel.

In liquids, an immediate manifestation of the surface tension is the shrinkage of the free surfaces towards a minimum compatible with the total volume of the liquid. On the other hand, in solids, the surface reduction at constant number of atoms induces a strain in the core, with a rapid increase of the corresponding elastic energy. In this case, more appropriate quantities can be introduced, the so-called surface stresses $s_{xx}$ and $s_{yy}$, obtained by analyzing the variations of the excess grand potential $\Omega_f^{ex}$ with respect to the film dimensions in the $x$ and $y$ directions. For instance along the $x$ direction, the variation of the total grand potential $\partial\Omega_f/\partial L_x$ is by definition the $x$ component $F_{xx}$ of the total lateral force exerted by the film in the $x$ direction. Note that this component is proportional to $L_y$: $F_{xx} = f_{xx}L_y$. This force decomposes into bulk plus surface contributions. The bulk contribution can be deduced from the bulk elastic properties: $\partial\Omega_b/\partial L_x = \sigma_{xx}H_\epsilon^G L_y$, where the stress tensor coefficient $\sigma_{xx}$ gives the $x$ component of the force per unit area due to a deformation along $x$. Therefore, $\partial\Omega_f^{ex}/\partial L_x = \partial(\Omega_f - \Omega_b)/\partial L_x = (f_{xx} - \sigma_{xx}H_\epsilon^G)L_y$. The $x$ component of the surface stress tensor along the $x$ direction is then defined as

$$ s_{xx} = \frac{1}{2} \frac{1}{L_y} \frac{\partial \Omega_f^{ex}}{\partial L_x} \bigg|_{T,\mu} = \frac{1}{2} (f_{xx} - \sigma_{xx} H_\epsilon^G), \qquad (2) $$

where a factor 1/2 is introduced to take into account the two faces of the film. Note that it is always $H_\epsilon^G = n\lambda_\epsilon$ that appears in Eq. (2), regardless of the choice of the position of the Gibbs dividing surface that can be different from the equimolar surface: $s_{xx}$ is an intrinsic quantity. One can similarly define $s_{yy}$ from the force per unit length along the $y$ direction $f_{yy} = F_{yy}/L_x$ and the bulk contribution $\sigma_{yy}H_\epsilon^G$.
Using $\Omega_f^{ex} = 2\gamma L_x L_y = 2L^2\gamma(1 + \epsilon_{xx})(1 + \epsilon_{yy})$, one gets the Shuttleworth's relation$^{52}$

$$s_{xx} = \gamma + \partial\gamma/\partial\epsilon_{xx}. \qquad (3)$$

The same relation holds for $s_{yy}$ in the y direction.

## D. Elasticity theory

Hooke's law states that, in the limit of small deformations, the stress tensor is linear with the strain tensor. The elastic stiffness $C$ is defined as the fourth-rank tensor that relates stress to strain. Using the Voigt notation to simplify the equations, the generalized Hooke's law writes$^{53}$

$$\sigma_i = \sum_{j=1}^{6} C_{ij} \epsilon_j \qquad (4)$$

with $\sigma_1 = \sigma_{xx}$, $\sigma_2 = \sigma_{yy}$, $\sigma_3 = \sigma_{zz}$, $\sigma_4 = \sigma_{yz}$, $\sigma_5 = \sigma_{xz}$, $\sigma_6 = \sigma_{xy}$, $\epsilon_1 = \epsilon_{xx}$, $\epsilon_2 = \epsilon_{yy}$, $\epsilon_3 = \epsilon_{zz}$, $\epsilon_4 = 2\epsilon_{yz}$, $\epsilon_5 = 2\epsilon_{xz}$, and $\epsilon_6 = 2\epsilon_{xy}$, where $\sigma_{\alpha\beta}$ and $\epsilon_{\alpha\beta}$ are the stress and strain tensors, respectively, for which $\alpha$ and $\beta$ equal $x$, $y$, or $z$. The $C_{ij}$'s are the coefficients of the elastic stiffness matrix in the Voigt notation. Since in the following we will consider only mechanical tests without shear, the indices $i$ and $j$ in Eq. (4) take values $\le 3$, giving

$$\begin{pmatrix} \sigma_1 \\ \sigma_2 \\ \sigma_3 \end{pmatrix} = \begin{pmatrix} C_{11} & C_{12} & C_{13} \\ C_{21} & C_{22} & C_{23} \\ C_{31} & C_{32} & C_{33} \end{pmatrix} \begin{pmatrix} \epsilon_1 \\ \epsilon_2 \\ \epsilon_3 \end{pmatrix}. \qquad (5)$$

Equation (5) shows that the stiffness coefficient $C_{ij}$ may be obtained directly by imposing a non-zero deformation $\epsilon_j$ only along the $j$ axis and measuring the average $i$ component of the stress tensor $\sigma_i$ (finite-strain method). Conversely, it might be more convenient to use the finite-stress method, where a constant stress is applied to the system and the average deformation is measured in different directions. The obtained tensor is the compliance matrix, that can be written in terms of Young's moduli ($E_i$) and Poisson's ratios ($\nu_{ij}$) as

$$\begin{pmatrix} \epsilon_1 \\ \epsilon_2 \\ \epsilon_3 \end{pmatrix} = \begin{pmatrix} 1/E_1 & -\nu_{21}/E_2 & -\nu_{31}/E_3 \\ -\nu_{12}/E_1 & 1/E_2 & -\nu_{32}/E_3 \\ -\nu_{13}/E_1 & -\nu_{23}/E_2 & 1/E_3 \end{pmatrix} \begin{pmatrix} \sigma_1 \\ \sigma_2 \\ \sigma_3 \end{pmatrix}. \qquad (6)$$
The compliance matrix is the inverse of the stiffness one in the thermodynamic limit where the iso-
choric and isostress statistical ensembles, used in the finite-strain and finite stress methods respec-
tively, are equivalent.³³ In our case, since this study is done at zero temperature, the equivalence
holds even for a finite system.

Not all elastic coefficients are independent. The matrices are symmetric, and, in the case of a
bulk crystal with cubic symmetry, only two elastic coefficients are independent if shear is ignored.
For the stiffness matrix, these independent coefficients are $C_{11} = C_{22} = C_{33}$ (this quantity will be
denoted $C_{ii}^b$ in the following), and $C_{12} = C_{13} = C_{23}$ (denoted $C_{i\neq j}^b$). For the compliance matrix,
the two quantities are the elastic modulus $E_1 = E_2 = E_3$ (denoted $E^b$) and the Poisson's ratio
$\nu_{12} = \nu_{13} = \nu_{23}$ (denoted $\nu^b$). For a film cut from a bulk crystal, with cubic symmetry, along
crystallographic planes, the tetragonal symmetry imposes that $C_{ij}$, for $i$ and $j \le 3$, has only four
independent coefficients, $C_{11} = C_{22}$, $C_{12}$, $C_{13} = C_{23}$, and $C_{33}$. For the compliance matrix, one has
the four quantities $E_1 = E_2$, $E_3$, $\nu_{12}$, and $\nu_{13} = \nu_{23}$. The remaining Poisson's ratios are given by
$\nu_{21} = \nu_{12}$, $\nu_{31} = \nu_{13}E_3/E_1$ and $\nu_{32} = \nu_{23}E_3/E_2$, with $\nu_{31} = \nu_{32}$ by tetragonal symmetry.

E. Strain and stress definitions for the film

The strain in the film is defined as follows. In the x and y directions parallel to the free surfaces, the film being homogeneous, the strain is equal to the relative deformation of the corresponding film dimension, $\epsilon_1 = (L_x - L)/L$ and $\epsilon_2 = (L_y - L)/L$, where $L_x$ and $L_y$ are the dimensions of the strained system, while $L$ corresponds to the unstrained (reference) bulk system.

In the z direction, two definitions of the strain are possible:

(i) The first one follows from the position of the two symmetric Gibbs equimolar surfaces, $\epsilon_3 = \delta H_\epsilon^G/H_0^G$, where $\delta H_\epsilon^G$ measures the Gibbs surface displacement due to film deformation. $\delta H_\epsilon^G$ is calculated by requiring that the excess quantity of matter $N^{ex}$ (or equivalently $n^{ex}$) remains constant during the deformation, and equals to zero for the Gibbs equimolar dividing surface. The position of the Gibbs surface then moves in inverse proportion to the density in the center, i.e. $\delta H_\epsilon^G/H_0^G = \delta\lambda_\epsilon/\lambda_b$, the relative variation of the interplanar distance in the center of the film upon deformation. This leads to a practical definition in terms of the distance $\lambda_\epsilon$ between the two innermost crystallographic planes in the deformed film, $\epsilon_3 = (\lambda_\epsilon - \lambda_b)/\lambda_b$, where the reference state is the relaxed bulk.

(ii) The second definition follows from the fact that, in simulations, it is not possible to control the
deformation in the core, whereas the position of the outermost layers can be monitored. Therefore, another definition of the film deformation along the z direction is introduced, $\epsilon'_3 = (H'_\epsilon - H'_0)/H'_0$, where $H'_0$ is the distance between the two outermost crystallographic planes for the reference film ($\epsilon_1 = \epsilon_2 = 0$ and free to relax along z), and $H'_\epsilon$ is the same distance when the film is deformed.

From a theoretical point of view, $\epsilon_3$ is the good quantity to measure the deformations. However, $\epsilon'_3$ is required for practical calculations using atomistic models. The objective of this work is to analyze the consequences of these two different definitions on the values of the elastic constants, and give the theoretical framework to convert the simulation results into quantities more suitable for the theory, and conversely.

Let us now consider the definition of the stress in the film. The stress $\sigma_1$ along the x direction derives from the x component $F_{xx}$ of the force that applies to any section (of area $S$) of the film normal to the x axis, $\sigma_1 = F_{xx}/S$. The area $S$ is defined by the Gibbs dividing surface, $S = L_y H_\epsilon^G$. A similar equation holds for $\sigma_2 = F_{yy}/(L_x H_\epsilon^G)$ with $F_{yy}$, the y component of the force along the y direction. Note that $\sigma_1$ and $\sigma_2$ explicitly depend on the position of the Gibbs surface: a choice different from the Gibbs equimolar dividing surface would lead to different values of the stress. However, the associated extensive quantities, such as the elastic energy or the total forces on the film, are unchanged. In the same way, the stress tensor $\sigma_3$ derives from the z component $F_{zz}$ of the force that applies to the free surface, $\sigma_3 = F_{zz}/(L_x L_y)$. As for $\sigma_1$ and $\sigma_2$, this definition implicitly depends on the position of the Gibbs surface, since it supposes that the force is applied exactly at the position of the Gibbs surface.

An important point to mention is that the measured forces, and therefore stresses, do depend on the variables used to describe the degrees of freedom of the system (generalized coordinates), in particular regarding the choice between $\epsilon_3$ and $\epsilon'_3$. In practice, the simulations provide the forces corresponding to the deformations $\epsilon'_3$. More details are given at the end of section II H.

## F. Relations between quantities derived from $\epsilon_3$ and $\epsilon'_3$

With the two definitions (i) and (ii) of the film deformations $\epsilon_3$ and $\epsilon'_3$ given in section II E, one can in principle calculate the two corresponding stress tensors, and finally the two sets of elastic constants of the film from Eq. (5). Obviously, depending on the adopted definition, one expects different values for the two sets. We now explicit the practical consequences of choosing definition (i) or (ii), by introducing a simple formalism allowing one to switch from one definition to the
other. For clarity, the prime symbol is used for the quantities measured in the atomistic simulations, where $\epsilon'_3$ is controlled through the monitoring of the outermost crystallographic planes, while non-prime quantities are those that would be obtained by imposing $\epsilon_3$. It is possible to relate prime and non-prime quantities, since $H'_e$ and $H'_e$ used to define $\epsilon_3$ and $\epsilon'_3$ are not independent. Let us introduce $\Delta_r$, the displacement of the outermost crystallographic plane (due to surface relaxation) with respect to the extrapolated position from the center of the film with the interplanar distance $\lambda_e$ (see Fig. 2). This surface relaxation is expected to depend on the internal deformation in the film, since, upon deformation, the interplanar distances in the core vary, with consequences on the forces applying to the outermost layers. By definition,

$$ \Delta_r(\epsilon) = 0.5[H'_e - (n-1)\lambda_e], \quad (7) $$

where the factor 0.5 takes into account the two free faces, and $(n-1)$ counts the number of interplanar spacings of the film. For the reference film, without deformation, $\Delta_r(0) = 0.5[H'_0 - (n-1)\lambda_b]$. Using the definitions (i) and (ii) above, the change in relaxation amplitude upon core deformation $d\Delta_r = \Delta_r(\epsilon) - \Delta_r(0)$ is given by $d\Delta_r = 0.5[H'_0\epsilon'_3 - (n-1)\lambda_b\epsilon_3]$. Introducing the partial derivatives with respect to $\epsilon$, $d\Delta_r = \sum_{i=1}^3 \frac{\partial\Delta_r}{\partial\epsilon_i}\epsilon_i$, one gets

$$ \epsilon'_3 = \frac{1}{H'_0} \left\{ 2 \frac{\partial \Delta_r}{\partial \epsilon_1} \epsilon_1 + 2 \frac{\partial \Delta_r}{\partial \epsilon_2} \epsilon_2 + \left[ (n-1)\lambda_b + 2 \frac{\partial \Delta_r}{\partial \epsilon_3} \right] \epsilon_3 \right\}. \quad (8) $$

This is a linear relation with respect to the $\epsilon_i$'s, with coefficients involving the partial derivatives calculated at zero deformation. Conversely, introducing the partial derivatives with respect to $\epsilon'$ at zero deformation, $d\Delta_r = \sum_{i=1}^3 \frac{\partial\Delta_r}{\partial\epsilon'_i}$, one has

$$ \epsilon_3 = \frac{1}{(n-1)\lambda_b} \left\{ -2 \frac{\partial \Delta_r}{\partial \epsilon'_1} \epsilon'_1 - 2 \frac{\partial \Delta_r}{\partial \epsilon'_2} \epsilon'_2 + \left[ H'_0 - 2 \frac{\partial \Delta_r}{\partial \epsilon'_3} \right] \epsilon'_3 \right\}. \quad (9) $$

It is useful to introduce the matrix formalism to switch easily between prime and non-prime deformations. With the additional definitions $\epsilon'_1 = \epsilon_1$ and $\epsilon'_2 = \epsilon_2$, Eq. (8) gives

$$ \begin{pmatrix} \epsilon'_1 \\ \epsilon'_2 \\ \epsilon'_3 \end{pmatrix} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ \frac{2}{H'_0} \frac{\partial \Delta_r}{\partial \epsilon_1} & \frac{2}{H'_0} \frac{\partial \Delta_r}{\partial \epsilon_2} & \frac{H''_0}{H'_0} \end{pmatrix} \begin{pmatrix} \epsilon_1 \\ \epsilon_2 \\ \epsilon_3 \end{pmatrix} \quad (10) $$

with $H''_0 = (n-1)\lambda_b + 2\frac{\partial\Delta_r}{\partial\epsilon_3}$. From Eq. (9),

$$ \begin{pmatrix} \epsilon_1 \\ \epsilon_2 \\ \epsilon_3 \end{pmatrix} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ \frac{-2}{H_0''} \frac{\partial \Delta_r}{\partial \epsilon_1'} & \frac{-2}{H_0''} \frac{\partial \Delta_r}{\partial \epsilon_2'} & \frac{H_0''}{H_0''} - \frac{2}{H_0''} \frac{\partial \Delta_r}{\partial \epsilon_3'} \end{pmatrix} \begin{pmatrix} \epsilon'_1 \\ \epsilon'_2 \\ \epsilon'_3 \end{pmatrix} \quad (11) $$
with $H_0''' = (n-1)\lambda_b$.

For small deformations, these relations are exact, with matrix coefficients that are independent of $\epsilon$ or $\epsilon'$ (the derivatives $\frac{\partial \Delta_r}{\partial \epsilon_i}$ and $\frac{\partial \Delta_r}{\partial \epsilon'_i}$ are calculated at $\epsilon_1 = \epsilon_2 = \epsilon_3 = 0$ and $\epsilon'_1 = \epsilon'_2 = \epsilon'_3 = 0$). We introduce the notation $T_i^j$ for the coefficients of the matrix $T$ that transforms the $\epsilon_i$ into the $\epsilon'_i$, and $T'^j_i$ for the coefficients of the matrix $T'$ that transforms the $\epsilon'_i$ into $\epsilon_i$,

$$ \epsilon'_i = \sum_{k=1}^{3} T_i^k \epsilon_k, \quad \epsilon_i = \sum_{k=1}^{3} T'^k_i \epsilon'_k. \qquad (12) $$

Obviously, $T' = T^{-1}$. These coefficients also correspond to the following partial derivatives,

$$ T_i^j = \frac{\partial \epsilon'_i}{\partial \epsilon_j}, \quad \text{and} \quad T'^j_i = \frac{\partial \epsilon_i}{\partial \epsilon'_j}. \qquad (13) $$

From Eq. (10), one has

$$ \begin{aligned} T_1^1 &= T_2^2 = 1 \\ T_1^2 &= T_1^3 = T_2^1 = T_2^3 = 0 \\ T_3^1 &= \frac{2}{H_0'} \frac{\partial \Delta_r}{\partial \epsilon_1}, \quad T_3^2 = \frac{2}{H_0'} \frac{\partial \Delta_r}{\partial \epsilon_2}, \quad T_3^3 = \frac{H_0''}{H_0'}. \end{aligned} \qquad (14) $$

Using Eq. (11), similar relationships are obtained for $T'^j_i$, involving the partial derivatives with respect to $\epsilon'_i$:

$$ \begin{aligned} T_1'^1 &= T_2'^2 = 1 \\ T_1'^2 &= T_1'^3 = T_2'^1 = T_2'^3 = 0 \\ T_3'^1 &= -\frac{2}{H_0'''} \frac{\partial \Delta_r}{\partial \epsilon_1'}, \quad T_3'^2 = -\frac{2}{H_0'''} \frac{\partial \Delta_r}{\partial \epsilon_2'}, \quad T_3'^3 = \frac{H_0'}{H_0'''} - \frac{2}{H_0'''} \frac{\partial \Delta_r}{\partial \epsilon_3'}. \end{aligned} \qquad (15) $$

This formalism will prove extremely useful to rationalize the molecular simulation results. For instance, transforming the measured $\sigma'_i$ into the $\sigma_i$ is straightforward thanks to the following property. Since the stresses are first derivatives of the free energy $A$ with respect to deformations, $\sigma_i = (1/V)\partial A/\partial \epsilon_i$ and $\sigma'_i = (1/V)\partial A/\partial \epsilon'_i$, where $V$ is the volume, they transform like the derivative operator $\partial/\partial \epsilon_i$. Differential summation rule gives

$$ \frac{\partial}{\partial \epsilon_i} = \sum_{k=1}^{3} \frac{\partial \epsilon_k'}{\partial \epsilon_i} \frac{\partial}{\partial \epsilon_k'} = \sum_{k=1}^{3} T_k^i \frac{\partial}{\partial \epsilon_k'}. \qquad (16) $$

Therefore, $\sigma_i = \sum_{k=1}^3 T_k^i \sigma'_k$.

The matrices $T$ and $T'$ involve the partial derivatives of $\Delta_r$ with respect to $\epsilon_i$ and $\epsilon'_i$. The derivatives relative to $\epsilon'_i$ can be straightforwardly determined from our simulation results using the
finite difference method. However, getting the derivatives relative to $\epsilon_i$ would be more involved. Therefore, we chose a general procedure that applies to any quantity, for partial derivatives of first and second order, and with respect to either $\epsilon'_i$ or $\epsilon_i$. The procedure consists in calculating the quantity of interest, denoted $Q$ (in this work, $Q = \Delta_r$ or $\gamma$), for a given set of imposed small deformations $\epsilon'_i$, and fit the data with a multivariate polynomial of second order,

$$Q(\epsilon'_1, \epsilon'_2, \epsilon'_3) = Q(0) + \sum_{i=1}^{3} Q^{(i)} \epsilon'_i + \frac{1}{2} \sum_{i,j=1}^{3} Q^{(i)(j)} \epsilon'_i \epsilon'_j \quad (17)$$

where $Q(0)$ means $Q(0,0,0)$, $Q^{(i)}$ and $Q^{(i)(j)}$ are constants corresponding to the first and second derivatives at origin, $Q^{(i)} = \partial Q / \partial \epsilon'_i$, and $Q^{(i)(j)} = \partial^2 Q / \partial \epsilon'_i \partial \epsilon'_j$. In practice, to map the $(\epsilon'_1, \epsilon'_2, \epsilon'_3)$ space, five values are chosen for each of the three deformations $\epsilon'_i \in \{0, \pm 5 \times 10^{-5}, \pm 10^{-4}\}$, resulting in a set of 125 points to fit the $Q$ data. The tetragonal symmetry of the film implies that $Q^{(1)} = Q^{(2)}$, $Q^{(1)(1)} = Q^{(2)(2)}$, $Q^{(1)(2)} = Q^{(2)(1)}$, and $Q^{(1)(3)} = Q^{(3)(1)} = Q^{(2)(3)} = Q^{(3)(2)}$.

In order to obtain the partial derivatives relative to $\epsilon_i$, the same procedure can be applied. The first operation consists in calculating the set of deformations $\epsilon_i$ corresponding to the $\epsilon'_i$. It is straightforward for $\epsilon_1$ and $\epsilon_2$. For $\epsilon_3$, one can either measure directly the deformation in the center of the film (with the two innermost layers), or use Eqs. (9), (11) or (12) if the $T'_i$ have been determined. Once the corresponding set of deformations $\epsilon_i$ is obtained, the $Q$ values are fitted with

$$Q(\epsilon_1, \epsilon_2, \epsilon_3) = Q(0) + \sum_{i=1}^{3} Q^{(i)} \epsilon_i + \frac{1}{2} \sum_{i,j=1}^{3} Q^{(i)(j)} \epsilon_i \epsilon_j \quad (18)$$

where $Q(0) = Q(0,0,0)$, $Q^{(i)} = \partial Q / \partial \epsilon_i$, and $Q^{(i)(j)} = \partial^2 Q / \partial \epsilon_i \partial \epsilon_j$. The same symmetry rules apply, $Q^{(1)} = Q^{(2)}$, $Q^{(1)(1)} = Q^{(2)(2)}$, $Q^{(1)(2)} = Q^{(2)(1)}$, and $Q^{(1)(3)} = Q^{(3)(1)} = Q^{(2)(3)} = Q^{(3)(2)}$. Note that the relations between the first derivatives $Q^{(i)}$ and $Q^{(i)}$ are given by Eq. (16),

$$Q^{(i)} = \sum_{k=1}^{3} T_k^i Q^{(k)}. \quad (19)$$

The second derivatives transform according to

$$\frac{\partial^2}{\partial \epsilon_i \partial \epsilon_j} = \sum_{k,l=1}^{3} \frac{\partial \epsilon'_k}{\partial \epsilon_i} \frac{\partial \epsilon'_l}{\partial \epsilon_j} \frac{\partial^2}{\partial \epsilon'_k \partial \epsilon'_l} = \sum_{k,l=1}^{3} T_k^i T_l^j \frac{\partial^2}{\partial \epsilon'_k \partial \epsilon'_l} \quad (20)$$

Therefore, one expects

$$Q^{(i)(j)} = \sum_{k,l=1}^{3} T_k^i T_l^j Q^{(k)(l)}. \quad (21)$$

Note that the quantity $Q$ and its derivatives may also depend on the film thickness $n$. In this case, we will write $Q_n$, $Q_n^{(i)}$, $Q_n^{(i)(j)}$, $Q'_n^{(i)}$, and $Q'_n^{(i)(j)}$.
G. Imposed strain method

For the bulk crystal, the two independent elastic parameters $E^b$ and $v^b$ are determined in a single mechanical test, consisting in imposing a finite deformation $\epsilon_1$ along the x direction while the box length in the y and z directions are kept unchanged (zero deformation with respect to the initial relaxed solid). The stresses $\sigma_i$ are measured in the three directions. By symmetry, one expects $\sigma_2 = \sigma_3$, and the ratio $\sigma_2/\epsilon_1 = \sigma_3/\epsilon_1$ gives the off-diagonal stiffness coefficient $C_{i\neq j}^b$, while $\sigma_1/\epsilon_1$ gives the diagonal coefficient $C_{ii}^b$. The bulk elastic modulus is deduced from $E^b = (C_{ii}^b + 2C_{i\neq j}^b)(C_{ii}^b - C_{i\neq j}^b)/(C_{ii}^b + C_{i\neq j}^b)$ and the bulk Poisson's ratio from $v^b = C_{i\neq j}^b/(C_{ii}^b + C_{i\neq j}^b)$.

For the thin film, we have seen that we have two possible definitions of the strain ($\epsilon_1, \epsilon_2, \epsilon_3$) and ($\epsilon'_1, \epsilon'_2, \epsilon'_3$) and the corresponding stresses ($\sigma_1, \sigma_2, \sigma_3$) and ($\sigma'_1, \sigma'_2, \sigma'_3$). We therefore have two sets of stiffness coefficients, $C_{ij}$ and $C'_{ij}$, defined by analogy with Eq. (5),

$$ \sigma_i = \sum_{j=1}^{3} C_{ij} \epsilon_j \quad \text{and} \quad \sigma'_i = \sum_{j=1}^{3} C'_{ij} \epsilon'_j. \qquad (22) $$

The tetragonal symmetry imposes that, in each case, only four coefficients are independent. Therefore, two mechanical tests are required.

The first mechanical test consists in imposing a finite deformation $\epsilon'_1 = \epsilon_1$ along the x direction, while the box length in the y direction is unchanged ($\epsilon'_2 = \epsilon_2 = 0$) and the outermost crystallographic planes are frozen ($\epsilon'_3 = 0$) (see Fig. 3). Measuring the three different stresses $\sigma'_1, \sigma'_2$, and $\sigma'_3$ in the x, y, and z directions, as given by the simulations, allows one to calculate the three independent coefficients $C'_{11}, C'_{21}$, and $C'_{31}$, deduced from $\sigma'_i = C'_{i1} \epsilon'_1$ for $i \in \{1, 2, 3\}$ [Eq. (22)]. It is emphasized that, even if it were possible to measure the $\sigma_i$, the corresponding $C_{i1}$ could not be deduced from $\sigma_i/\epsilon_1$ since $\epsilon_3 \neq 0$.

The second test consists in imposing a finite deformation $\epsilon'_3$, i.e. an imposed displacement of the two outermost atomic planes, while the box dimensions in the x and y directions are unchanged, $\epsilon'_1 = \epsilon_1 = 0$ and $\epsilon'_2 = \epsilon_2 = 0$ (see Fig. 3). The stress $\sigma'_3$, given by the simulations, allows one to calculate $C'_{33} = \sigma'_3/\epsilon'_3$. The two identical stresses $\sigma'_1$ and $\sigma'_2$ give $C'_{13} = \sigma'_1/\epsilon'_3 = C'_{23} = \sigma'_2/\epsilon'_3$. Although one can easily measure the displacement of the two innermost layers to deduce $\epsilon_3$, it is not possible to obtain $\sigma_1, \sigma_2$, or $\sigma_3$ directly from the simulations. Therefore, as previously, the coefficients $C_{i3} = \sigma_i/\epsilon_3$ are not directly accessible to simulations.

The symmetry of the $C'_{ij}$ imposes that $C'_{31}$, obtained with the first mechanical test, should be equal to the $C'_{13}$ obtained with the second test. This is used to check the consistency between the
two mechanical tests and to validate our simulation results.

## H. Numerical Details

The questions raised in this work are generic to nanometric systems, and can, at least in a first approach, be studied with simple pairwise additive interactions such as the Lennard-Jones (LJ) potential.⁵⁴⁻⁵⁹ Considering this potential is also useful since it is frequently integrated into widely used force fields.⁵⁵,⁵⁹ Of course, more realistic potentials exist, such as EAM, TBSMA or Tersoff.⁶⁰⁻⁶² Such potentials are able, for instance, to predict the correct sign of the displacements associated to surface relaxation, with quantitative agreement. They will be considered in further studies. For this work, the simple LJ model is sufficient to qualitatively reproduce the experimental observation that the free surfaces relax due to the imbalance of forces caused by the missing atoms above the surface. It shows that surface relaxation occurs over a few atomic layers, in relation with the interaction range, which is expected to have a large effect in ultrathin films and in multilayers. This model is used to study the elastic constants of a thin solid film as a function of its thickness. It is shown that the formalism introduced in section II F is a powerful tool to convert the elastic constants as given by the simulation into the corresponding intrinsic quantities as defined from a theoretical point of view, enabling a simple analysis in terms of bulk + surface contributions.

In this work, interatomic forces derive from the (12-6) LJ potential

$$V_{LJ} = 4\epsilon_{LJ} \left[ \left( \frac{\sigma_{LJ}}{r} \right)^{12} - \left( \frac{\sigma_{LJ}}{r} \right)^{6} \right] \quad (23)$$

where $\epsilon_{LJ}$ is the energy well and $\sigma_{LJ}$ is the van der Waals or kinetic diameter. The potential is cut at $r_c = 4\sigma_{LJ}$, and linearly shifted so that the forces are continuous. Similar potential and cutoff have been used in the literature, allowing comparisons.²⁷,³⁵ In this work, all calculations are performed at zero temperature, by energy minimization (conjugate gradient method) using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) molecular dynamics software.⁶³ At this temperature, the stable LJ solid phase is the cubic FCC crystal.

The LJ parameters are used as units: $\sigma_{LJ}$ for distances, $\epsilon_{LJ}$ for energies, $\epsilon_{LJ}/\sigma_{LJ}$ for forces, $\epsilon_{LJ}/\sigma_{LJ}^2$ for surface tension and surface stress, $\epsilon_{LJ}/\sigma_{LJ}^3$ for stresses and stiffness coefficients, and the inverse for moduli. Depending on the context, the quantities are either given with their units, in expressions such as $r_c = 4\sigma_{LJ}$, or given in dimensionless form (in figures for instance).

We first consider the bulk solid. Using periodic boundary conditions minimizes finite size effects. At zero temperature, using a system at least twice as large as the interactions cutoff makes
<table><caption>TABLE I. Bulk quantities for the LJ crystal at zero temperature.</caption><thead><tr><th>property</th><th>simulation result</th></tr></thead><tbody><tr><td>a<sub>b</sub></td><td>1.54554 &sigma;<sub>LJ</sub></td></tr><tr><td>&lambda;<sub>b</sub> = a<sub>b</sub>/2</td><td>0.77277 &sigma;<sub>LJ</sub></td></tr><tr><td>u<sub>b</sub></td><td>-8.11508 &epsilon;<sub>LJ</sub></td></tr><tr><td>C<sup>b</sup><sub>ii</sub></td><td>100.57 &epsilon;<sub>LJ</sub>/&sigma;<sup>3</sup><sub>LJ</sub></td></tr><tr><td>C<sup>b</sup><sub>i&ne;j</sub></td><td>57.492 &epsilon;<sub>LJ</sub>/&sigma;<sup>3</sup><sub>LJ</sub></td></tr><tr><td>E<sup>b</sup></td><td>58.751 &epsilon;<sub>LJ</sub>/&sigma;<sup>3</sup><sub>LJ</sub></td></tr><tr><td>&nu;<sup>b</sup></td><td>0.36372</td></tr></tbody></table>

it is possible to reach the exact values for the mechanical properties, without finite size effect. We use seven unit cells in the three space directions (the axes are parallel to the crystallographic directions), giving a total number of atoms $N = 1372$ and resulting in a box size $\ge 10\sigma_{LJ}$, larger than $2r_c$. The solid is relaxed at zero temperature and pressure. The measured lattice parameter is $a_b = 1.54554\sigma_{LJ}$, the interplanar spacing is $\lambda_b = 0.5a_b = 0.77277\sigma_{LJ}$, and the box edges measure $L = 7a_b = 10.8188\sigma_{LJ}$. The potential energy per atom is $u_b = -8.11508\epsilon_{LJ}$. This cubic crystal is used to calculate the bulk elastic constants. The diagonal and off-diagonal stiffness coefficients are $C_{ii}^b = 100.57\epsilon_{LJ}/\sigma_{LJ}^3$ and $C_{i\ne j}^b = 57.492\epsilon_{LJ}/\sigma_{LJ}^3$. The Young modulus is $E_b = 58.751\sigma_{LJ}^3/\epsilon_{LJ}$ and the Poisson's ratio is $\nu_b = 0.36372$. These values are in agreement with previous independent works.$^{27,35}$ They are reported in Table I.

This relaxed crystal is then used to build the films by extracting a slab parallel to the crystallographic (001) planes, and characterized by its number of crystallographic planes $n$. The thinnest film considered in this study contains four crystallographic planes. Since it was desirable to evidence how the bulk properties are approached, we also consider very large films, thicker than the initial bulk solid comprising only seven unit cells in the $z$ direction (14 crystallographic planes). These thick films are built thanks to the periodic images along $z$ of the initial crystal. The thickest film reaches two hundred atomic layers. To simulate an isolated film, the periodic boundary conditions in the $z$ direction perpendicular to the surfaces are broken by introducing an empty space much larger than the interaction range. On the other hand, the simulation box is not modified in the $x$ and $y$ directions parallel to the surfaces, and the periodic boundary conditions are maintained to minimize finite size effects. The free surfaces of the film are allowed to relax along $z$, while the $x$
and y dimensions are kept unchanged. This results in a negative pressure in the x and y directions, consequence of the surface stress. This partially relaxed system is taken as the reference state to measure deformations in the film. In this reference state, for a film large enough, the atomic lattice in the central region of the film is identical to that in the initial bulk crystal: therefore, the deformations $\epsilon_i$ are all equal to zero. This condition would not be met if the film was fully relaxed (i.e. including relaxation of the lateral dimensions $L_x$ and $L_y$).

The stiffness coefficients are determined following an experimental route, i.e. as the initial slope (linear regime) of the stress-strain curves.³⁵ Typically, finite relative deformations around $2 \times 10^{-5}$ are applied. Such deformations are small enough to be in the linear regime, but still large enough to avoid errors arising when calculating differences that are too small. The calculations have been limited to simple deformations without shear, at zero temperature, by energy minimization (conjugate gradient method) using the LAMMPS molecular dynamics software.⁶³ In particular, the stress tensor is obtained with the default compute-pressure command. One has to be aware of the fact that the pressure tensor given by the LAMMPS software involves the total volume of the simulation box. It is thus necessary to rescale the result to the volume of the film as defined by the Gibbs surfaces ($\sigma'_i = -P_i^{\text{lammps}} L_z^{\text{box}} / H_0^G$, where $P_i^{\text{lammps}}$ is the pressure given by LAMMPS and $L_z^{\text{box}}$ is the size of the simulation box along z).

It is emphasized that, since the mechanical tests are performed by monitoring the deformations $\epsilon'_i$, the obtained stresses correspond to $\sigma'_i$. This has been numerically checked with two different procedures. In the first one, a constant outward force $f_{zz}$, typically $= 2 \times 10^{-3} \epsilon_{LJ} / \sigma_{LJ}$, was applied to each atom of the two opposite outermost layers of the film. This results in the stress $\sigma'_3 = 2f_{zz}/a_b^2$ (two atoms per lattice area). The numerical value is compared to the stress directly given by the simulations, after normalization by the volume of the film. Note that, in this case, the volume of the film is the one comprised exactly between the two outermost layers, since the forces are applied on the atoms and not on the Gibbs surface. One gets numerically exactly the same values as $\sigma'_3 = 2f_{zz}/a_b^2$ within the machine precision. The second procedure consists in comparing the stress directly given by the simulations when imposing a small relative deformation $\epsilon'_3 = 2 \times 10^{-5}$ to the stress $\sigma'_3$ deduced from the evolution of the internal energy of the film $U_f$ by numerical derivation at the same deformation: $\sigma'_3 = (1/V) \partial U_f / \partial \epsilon'_3$. As can be seen in Fig. 4, the stresses given by the simulations correspond exactly to the $\sigma'_3$ obtained by derivation. For comparison, the corresponding $\sigma_3$ are also reported, showing a noticeable difference.
### III. RESULTS

#### A. Surface relaxation

Upon film creation by cutting a piece of bulk crystal with the lateral dimensions $L_x$ and $L_y$ kept identical to the initial bulk value $L$, one observes a small outward displacement of the surface crystallographic planes along $z$ (surface relaxation). Quantitatively, the displacements are measured with respect to the ideal position of the planes in the bulk phase, i.e. calculated with the interplanar distance $\lambda_e = \lambda_b$ in this case since the strain in the center of the film is zero. Eq. (7) gives the formula for the outermost plane, but it is possible to generalize to any crystallographic plane. These generalized surface relaxation displacements are denoted $\Delta(m)$, where $m$ labels the crystallographic plane, numbered from the surface towards the center of the film ($m = 1$ for the outermost layer). The results are given in Fig. 5 for three film thicknesses $n = 10, 20$ and $30$, and for the first six layers (limited to five for the thinnest film with $n = 10$ because of the symmetry). The first point to note, is that all data superimpose, in agreement with the expected behavior for a pure surface effect independent of the film thickness. Furthermore, as can be seen, the amplitude of the displacements decreases rapidly when $m$ increases. It is negligible beyond the third lattice plane ($\le 10^{-4}$). Plotting the data in semi-log scale (see inset in Fig. 5) shows that the amplitude of the displacements decays exponentially with the layer label $m$. We introduce the following equation,

$$ \Delta(m) = \Delta_r e^{-(m-1)/m_r} \qquad (24) $$

where $(m-1)$ is introduced in the exponent so that $\Delta_r$ corresponds to the displacement of the outermost layer $m = 1$, in coherence with the definition given in Eq. (7). The parameter $m_r$ measures the typical exponential decay distance in units of interplanar distance. The solid lines in Fig. 5 are best fits with Eq. (24) performed with data from the two largest systems $n = 20$ and $30$, giving $\Delta_r = 2.06137 \times 10^{-2} \sigma_{LJ}$ and $m_r = 0.530$. For the thinnest film $n = 10$, the values are slightly different, but the differences are less than $10^{-4}$: the finite size effects are negligible. The first point to be mentioned is that the displacement of the outermost plane associated to the surface relaxation is quite small relative to the interplanar spacing: $\Delta_r/\lambda_b = 2.66 \times 10^{-2}$. Furthermore, $m_r = 0.53$ means that the amplitude of the displacement decreases by a factor 6.6 between layers $n$ and $n+1$. Thus, the surface relaxation is significant only for the very first layers.

When the film is deformed, the decay coefficient $m_r$ is only marginally varying, while the
amplitude of the relaxation $\Delta_r$ is significantly affected. We thus consider $m_r$ constant and focus on the variations of $\Delta_r$, as defined by Eq. (7), with the deformations. This quantity is very important for the following, in particular to calculate the coefficients $T_i^j$ and $T'_i^j$ which involve the first order derivatives of $\Delta_r(\varepsilon)$ [see Eqs. (10), (11), (14) or (15)]. These partial derivatives, relative to $\varepsilon$ and $\varepsilon'$, are calculated using Eqs. (17) and (18), and for different film thicknesses, $n = 14, 20, 30$ and $50$. They are denoted $\Delta_{r,n}^{(i)} = \partial\Delta_{r,n}/\partial\varepsilon_i$ and $\Delta'_{r,n}^{(i)} = \partial\Delta_{r,n}/\partial\varepsilon'_i$ to be consistent with the notations introduced in Eqs. (17) and (18). The results are reported as a function of the inverse number of layers $1/n$ in Fig. 6.

The first observation is that the derivatives relative to directions $x$ or $y$ are perfectly undistinguishable, as expected from the film lattice symmetry. One also observes that, although $\Delta_r$ depends weakly on $n$ (relative variations are $\le 10^{-5}$, see Fig. 5 for $m = 1$), the derivatives $\Delta_r^{(i)}$ and $\Delta'_r^{(i)}$ may exhibit noticeable film thickness dependence, in particular for the derivatives relative to the deformations $\varepsilon'_i$. At first sight, all variations are essentially proportional to $1/n$. Furthermore, extrapolations of the derivatives to infinite thickness ($n \to \infty$) seem to converge towards identical values for prime and non-prime quantities. This is perfectly in line with the idea that, for very thick films ($n \to \infty$), $\varepsilon'_3$ tend to become identical to $\varepsilon_3$. Another important observation is that the first derivatives relative to $\varepsilon_i$ are insensitive to the film thickness, $\Delta_{r,n}^{(i)} = \Delta'_r^{(i)}$. Therefore, it is much more convenient to use the coefficients $T_i^j$ given by Eqs. (14), where the involved derivatives $\Delta_r^{(i)}$ are constants, instead of the $T'_i^j$ that involve the $n$-dependent $\Delta_{r,n}^{(i)}$ [see Eqs. (15)]. However, in some cases, the $T'_i^j$ may be required. In this case, instead of using Eqs. (15), it is more efficient to perform a direct matrix inversion of the $T_i^j$ from Eq. (10),

$$ \begin{pmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3 \end{pmatrix} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ -\frac{2}{H_0''} \frac{\partial \Delta_r}{\partial \varepsilon_1} & -\frac{2}{H_0''} \frac{\partial \Delta_r}{\partial \varepsilon_2} & \frac{H_0'}{H_0''} \end{pmatrix} \begin{pmatrix} \varepsilon'_1 \\ \varepsilon'_2 \\ \varepsilon'_3 \end{pmatrix} \quad (25) $$

where now the dependence with the film thickness is entirely in $H_0' = (n-1)\lambda_b + 2\Delta_r$ and $H_0'' = (n-1)\lambda_b + 2\frac{\partial \Delta_r}{\partial \varepsilon_3}$. The Eqs. (15) are then replaced by

$$ \begin{align} T'_1 &= T'_2 = 1 \nonumber \\ T'_3 &= T'_1 = T'_2 = T'_3 = 0 \tag{26} \\ T'_3 &= -\frac{2}{H_0''} \frac{\partial \Delta_r}{\partial \varepsilon_1}, \quad T'_3 = -\frac{2}{H_0''} \frac{\partial \Delta_r}{\partial \varepsilon_2}, \quad T'_3 = \frac{H_0'}{H_0''}. \nonumber \end{align} $$

Figure 7 shows the evolution of the coefficients $T_i^j$ and $T'_i^j$ with the inverse of the film thickness
$1/n$. Only the non trivial coefficients ($\neq 0$ or 1) are shown, i.e. $T_3^i$, given by Eqs. (14), and $T_3^{i'}$ given either by Eqs. (15) or by Eqs. (26) corresponding to the $(T^{-1})_3^i$. As can be seen, the $x \leftrightarrow y$ symmetry is respected ($T_3^1 = T_3^2$ and $T_3^{1'} = T_3^{2'}$). Furthermore, $(T^{-1})_3^i$ are perfectly superimposed on $T_3^{i'}$, confirming that Eqs. (26) are equivalent to Eqs. (15). One can also see that $T_3^3$ and $T_3^{3'} \to 1$ for $n \to \infty$, while the other coefficients vanish. Therefore, for infinitely large films, the two matrices $T$ and $T'$ converge to the unit matrix.

Before giving the numerical values for the first derivatives $\Delta_r^{(i)}$, it is important to realize that the $\Delta_r^{(i)}$ and $\Delta_{r,n}^{(i)}$ are linked by the general relations given by Eq. (19), involving the $\Delta_r^{(i)}$ themselves. A self consistent fit can therefore be performed on the data displayed in Fig. 6, where the $\Delta_r^{(i)}$ are supposed to be constants independent of $n$, and are taken as the free parameters of the fit. The $n$-dependence of the $\Delta_{r,n}^{(i)}$ is then entirely determined by inversion of Eq. (19). The corresponding best fit appears as solid lines in Fig. 6, showing an excellent global agreement with the data. The obtained values of the free parameters $\Delta_r^{(i)}$ are given in Table II. Note that they also correspond to the values of $\Delta_{r,n}^{(i)}$ for $n \to \infty$, as can be seen in Fig. 6, and from the fact that $T'$ converges to the unit matrix for $n \to \infty$. This last property implies that for any $n$-dependent quantity $Q_n$, the first and second derivatives relative to $\epsilon$ or $\epsilon'$ are identical in the $n \to \infty$ limit. We therefore adopt the unique notation $Q_\infty^{(i)}$ and $Q_\infty^{(i)(j)}$, even for the derivatives relative to $\epsilon'$. It is important to note that the knowledge of these two constants ($\Delta_{r,\infty}^{(1)} = \Delta_{r,\infty}^{(2)} = 0.104773\sigma_{LJ}$ and $\Delta_{r,\infty}^{(3)} = 0.075100\sigma_{LJ}$, see Table II) determines entirely the matrices $T_i^j$ and $T_i^{j'}$ thanks to Eqs. (14) and (26), and for any film thickness $n$.

## B. Thin film elastic constants as obtained from simulations

The two different mechanical tests (imposed strain method, see Fig. 3) are applied to the film to determine the elastic coefficients. It is reminded that, since the monitored deformations in simulations are the $\epsilon'_i$, the measured stresses are the $\sigma'_i$, and the measured stiffness coefficients are the $C'_{ij}$. The results of our simulations are given in Fig. 8 as a function of the inverse film thickness $1/n$. The data for the infinite system ($1/n = 0$) are those obtained previously for the bulk crystal, $C_{ii}^{tb} = C_{ii}^b = 100.57\epsilon_{LJ}/\sigma_{LJ}^3$ and $C_{i\neq j}^{tb} = C_{i\neq j}^b = 57.492\epsilon_{LJ}/\sigma_{LJ}^3$, see Table I.

All data exhibit smooth variations, consistent with the bulk values in the $n \to \infty$ limit. The anisotropy of the film (tetragonal symmetry) appears as a splitting of the bulk coefficients, $C_{ii}^b$ into $C'_{11} = C'_{22} \neq C'_{33}$, and $C_{i\neq j}^b$ into $C'_{12} \neq C'_{13} = C'_{23}$. The gap increases when the film is getting
<table><caption>TABLE II. First and second derivatives of the two dimensionless quantities $Q = Δ_r/σ_{LJ}$ and $Q = γσ_{LJ}^2/ε_{LJ}$, extrapolated to the infinitely thick Lennard-Jones film $n → ∞$ from fits performed for finite films with Eqs. (17) and (18).</caption><thead><tr><th>Q</th><th>Δ<sub>r</sub>/σ<sub>LJ</sub></th><th>γσ<sup>2</sup><sub>LJ</sub>/ε<sub>LJ</sub></th></tr></thead><tbody><tr><td>Q<sub>∞</sub>(0)</td><td>0.0206137</td><td>2.71731</td></tr><tr><td>Q<sub>∞</sub><sup>(1)</sup> = Q<sub>∞</sub><sup>(2)</sup></td><td>0.104773</td><td>-2.22208</td></tr><tr><td>Q<sub>∞</sub><sup>(3)</sup></td><td>0.075100</td><td>-0.000104</td></tr><tr><td>Q<sub>∞</sub><sup>(1)(1)</sup> = Q<sub>∞</sub><sup>(2)(2)</sup></td><td>0.13373</td><td>-14.897</td></tr><tr><td>Q<sub>∞</sub><sup>(3)(3)</sup></td><td>0.76064</td><td>-31.298</td></tr><tr><td>Q<sub>∞</sub><sup>(1)(2)</sup> = Q<sub>∞</sub><sup>(2)(1)</sup></td><td>1.29895</td><td>15.614</td></tr><tr><td>Q<sub>∞</sub><sup>(1)(3)</sup> = Q<sub>∞</sub><sup>(3)(1)</sup></td><td>0.49283</td><td>-17.889</td></tr></tbody></table>

thinner. Quite remarkably, the variations of the coefficients with $1/n$ is close to linear over the whole thickness range, i.e. down to four atomic layers only.

The stiffness coefficients can be converted into the Young's moduli $E'_i$ and Poisson's ratios $v'_{ij}$ by matrix inversion (stiffness to compliance matrix inversion). The results are displayed in Fig. 9. One first observes that the two elastic moduli $E'_1$ and $E'_3$ decrease when the film thickness decreases. The variation of the out-of-plane modulus $E'_3$ is limited to 10% for the thinnest film $n = 4$, while the in-plane modulus $E'_1$ decreases by 40%. One also observes that, unlike for the $C_{ij}$ coefficients, the variations of the moduli are not linear with the inverse film thickness. On the other hand, the variations of the Poisson's ratios are remarkably linear with $1/n$, with strongest variations for the in-plane coefficient $v'_{12}$ (+75% for the thinnest film $n = 4$). Another remarkable point is that the mean $(v'_{12} + v'_{31} + v'_{13})/3$ is constant over the whole thickness range, from $n = 4$ to $∞$. We already noticed this property in a previous paper,⁶⁴ that seems to be also compatible with data obtained in continuum approaches.⁶⁵
C. Intrinsic elastic constants of the thin film

The elastic coefficients $C'_{ij}$ given by the simulations are certainly not the most appropriate to characterize the mechanical properties of the film, in particular when considering their dependence with the thickness. We have indeed seen for instance that the derivatives of $\Delta_r$ exhibit no dependence with $n$ only when they are calculated with respect to the $\epsilon_i$ deformations. This is why it is desirable to calculate the $C_{ij}$.

The $C_{ij}$ can easily be deduced from the $C'_{ij}$. Since the stress tensor is the first derivative of the Helmholtz free energy $A$ with respect to the strain tensor, the stiffness matrix is a second order derivative of $A$ with respect to the deformations,

$$C'_{ij} = \frac{1}{V} \frac{\partial^2 A}{\partial \epsilon'_i \partial \epsilon'_j} \quad \text{and} \quad C_{ij} = \frac{1}{V} \frac{\partial^2 A}{\partial \epsilon_i \partial \epsilon_j}, \qquad (27)$$

where $V$ is the volume of the sample. Therefore, the transformation rules for the stiffness coefficients have to follow Eq. (21). Therefore, the $C_{ij}$ can be deduced from the $C'_{ij}$ according to

$$ \begin{aligned} C_{11} &= C'_{11} + 2T_3^1 C'_{31} + T_3^1 T_3^1 C'_{33} \\ C_{12} &= C'_{12} + T_3^1 C'_{32} + T_3^2 C'_{13} + T_3^1 T_3^2 C'_{33} \\ C_{13} &= T_3^3 C'_{13} + T_3^1 T_3^3 C'_{33} \\ C_{33} &= T_3^3 T_3^3 C'_{33}. \end{aligned} \qquad (28) $$

The results are given in Fig. 10. As previously, the data for a finite thickness exhibit a splitting due to the lower symmetry of the film, with a vanishing amplitude when approaching the $n \to \infty$ limit. The bulk values for the infinite crystal are obviously unaffected by the transformation since it equals unity for $1/n = 0$. It is emphasized that the conversion of $C'_{ij}$ into $C_{ij}$ has not removed the dependence with the film thickness, unlike what was observed for the $\Delta_r^{(i)}$.

As previously, the $C_{ij}$ allow one to calculate the Young's moduli and the Poisson's ratios by matrix inversion. They are displayed in Fig. 11. The evolution of the data with the film thickness follows essentially the same behavior as for the $E'_i$ and $v'_{ij}$. In particular, for the Poisson's ratios, one observes again that the mean $(v_{12} + v_{31} + v_{13})/3$ is constant over the whole thickness range, although one may note a small departure for the thinnest films.
IV. DISCUSSION

The main idea in this discussion is to decompose the elastic constants of the film into a sum of a volume bulk-like term and a surface contribution. Starting from the general relation giving the elastic constants as second order strain derivatives [Eq. (27)], and decomposing the free energy A into the bulk elastic term $A_{el} = 0.5H_0^G L^2 \sum_{kl} C_{kl}^b \epsilon_k \epsilon_l$ plus the surface free energy term $2\gamma L_x L_y$, one gets

$$C_{ij} = C_{ij}^{b} + \frac{2}{H_{0}^{G}} \frac{\partial^{2}}{\partial \epsilon_{i} \partial \epsilon_{j}} [\gamma(1 + \epsilon_{1})(1 + \epsilon_{2})] \quad (29)$$

and

$$C'_{ij} = \sum_{kl} T'_{k}{}^{i} T'_{l}{}^{j} C_{kl}^{b} + \frac{2}{H_{0}^{G}} \frac{\partial^{2}}{\partial \epsilon'_{i} \partial \epsilon'_{j}} [\gamma(1 + \epsilon'_{1})(1 + \epsilon'_{2})] \quad (30)$$

Both sets of equations require the knowledge of $\gamma$ and its first and second derivatives with respect to $\epsilon_i$ or $\epsilon'_i$. We calculate $\gamma$ with Eq. (1) for a complete set of deformations $\epsilon'_i \in \{0, \pm 5 \times 10^{-5}, \pm 10^{-4}\}$, giving 125 points for the fits with Eqs. (17) and (18), and for four film thicknesses $n = 14, 20, 30$ and $50$. The obtained values of $\gamma$ and its first derivatives at zero deformation ($\epsilon = 0$) are given in Fig. 12, while the second derivatives at $\epsilon = 0$ are displayed in Fig. 13. As can be seen in Fig. 12, at zero deformation, $\gamma$ is insensitive to the film thickness. One gets $\gamma(0) = 2.71731\epsilon_{LJ}/\sigma_{LJ}^2$ for all system sizes, with differences $\le 10^{-8}$. One also observes that the $x \leftrightarrow y$ symmetry is respected for the first derivatives, and that the values are essentially independent of the film thickness $n$ (pay attention to the scale). Another point that needs to be emphasized is that the derivatives $\gamma^{(3)}$ and $\gamma'^{(3)}$ are very close to zero. Furthermore, the first derivatives relative to $\epsilon_i$ and $\epsilon'_i$ are observed to be identical, $\gamma^{(i)} = \gamma'^{(i)}$, $i \in \{1, 2, 3\}$. This is a consequence of Eq. (19) applied to $\gamma$, giving $\gamma_n^{(3)} = T_3^3 \gamma_n'^{(3)} \simeq \gamma_n'^{(3)} \simeq 0$, and $\gamma_n^{(1)} = \gamma_n'^{(1)} + T_3^1 \gamma_n'^{(3)} \simeq \gamma_n'^{(1)}$ and $\gamma_n^{(2)} = \gamma_n'^{(2)} + T_3^2 \gamma_n'^{(3)} \simeq \gamma_n'^{(2)}$. The average values for the three largest thicknesses are taken as good estimates of the corresponding extrapolations to infinite thickness. They are shown as horizontal lines in Fig. 12 and given in Table II. As for $\Delta_{r,\infty}^{(i)}$, we adopt the unique notation $\gamma_{\infty}^{(i)}$ to denote the common value of $\gamma_n^{(i)}$ and $\gamma_n'^{(i)}$ for $n \to \infty$. Before moving on to the second derivatives, it is interesting to see that these results can be used to calculate the surface stress $s_{xx}$ from the Shuttleworth relation, Eq. (3): $s_{xx} = \gamma_{\infty}(0) + \gamma_{\infty}^{(1)} = 0.4952\epsilon_{LJ}/\sigma_{LJ}^2$. For comparison, the quantity given by Eq. (2) has also been calculated, at zero deformation, from the forces given by the simulations. Note that, since we are working at zero strain in the film, the bulk contribution vanishes. The result is independent of the film thickness to within $10^{-4}$, in agreement with the definition of a surface quantity. The
obtained value, $s_{xx} = 0.4953\epsilon_{LJ}/\sigma_{LJ}^2$, is in perfect agreement with that obtained previously from the Shuttleworth relation.

The second order derivatives, shown in Fig. 13, exhibit a behavior reminiscent of what was observed for the $\Delta_r^{(i)}$ and $\Delta_{r,n}^{(i)}$. In particular, the variations are almost linear with $1/n$ for the $\gamma_n^{(i)(j)}$, while the corresponding non prime $\gamma_n^{(i)(j)}$ are essentially constant, both converging towards identical values at infinite thickness $1/n = 0$. As for the derivatives of $\Delta_{r,n}$, it is possible to perform a global fit of the data, with the $\gamma_n^{(i)(j)}$ taken as free parameters independent of $n$, and predicting the evolution of the $\gamma_n^{(i)(j)}$ with the help of

$$Q^{(i)(j)} = \sum_{k,l=1}^{3} T'_{k}{}^{i} T'_{l}{}^{j} Q^{(k)(l)} \quad (31)$$

[obtained by inversion of Eq. (21)] applied to $Q = \gamma$. As for $\Delta_r$, these constant values (free parameters) are denoted $\gamma_{\infty}^{(i)(j)}$ since they also correspond to the limiting values of $\gamma_n^{(i)(j)}$ and $\gamma_n^{(i)(j)}$ for $n \to \infty$. The relations are slightly more complex than for the $\Delta_r$,

$$ \begin{aligned} \gamma_n^{(1)(1)} &= \gamma_{\infty}^{(1)(1)} + 2T_3^{\prime 1}\gamma_{\infty}^{(3)(1)} + T_3^{\prime 1}T_3^{\prime 1}\gamma_{\infty}^{(3)(3)} \\ \gamma_n^{(1)(2)} &= \gamma_{\infty}^{(1)(2)} + T_3^{\prime 1}\gamma_{\infty}^{(3)(2)} + T_3^{\prime 2}\gamma_{\infty}^{(1)(3)} + T_3^{\prime 1}T_3^{\prime 2}\gamma_{\infty}^{(3)(3)} \\ \gamma_n^{(1)(3)} &= T_3^{\prime 3}\gamma_{\infty}^{(1)(3)} + T_3^{\prime 1}T_3^{\prime 3}\gamma_{\infty}^{(3)(3)} \\ \gamma_n^{(3)(3)} &= T_3^{\prime 3}T_3^{\prime 3}\gamma_{\infty}^{(3)(3)} \end{aligned} \quad (32) $$

The corresponding results are shown as lines in Fig. 13. As can be seen, the agreement is quite good. The obtained free parameters, which also correspond to the common values of the derivatives relative to $\epsilon_i$ and $\epsilon_i'$ for $n \to \infty$ are given in Table II.

We are now in a position to interpret our numerical results for the $C_{ij}$. From Eq. (29), one has

$$ \begin{aligned} C_{11} &= C_{ii}^b + \frac{2}{H_0^G} \left( \gamma_{\infty}^{(1)(1)} + 2\gamma_{\infty}^{(1)} \right) \\ C_{12} &= C_{i \neq j}^b + \frac{2}{H_0^G} \left( \gamma_{\infty}^{(1)(2)} + \gamma_{\infty}^{(1)} + \gamma_{\infty}^{(2)} + \gamma_{\infty} \right) \\ C_{13} &= C_{i \neq j}^b + \frac{2}{H_0^G} \left( \gamma_{\infty}^{(1)(3)} + \gamma_{\infty}^{(3)} \right) \\ C_{33} &= C_{ii}^b + \frac{2}{H_0^G} \gamma_{\infty}^{(3)(3)} \end{aligned} \quad (33) $$

with $H_0^G = n\lambda_b$. The corresponding relations for the $C_{ij}'$ are deduced from Eq. (30) where the $T_i^{\prime j}$ are given by Eqs. (26) and the $\gamma_n^{(i)(j)}$ are converted into $\gamma_{\infty}^{(i)(j)}$ thanks to Eqs. (32). It is also
possible to deduce the $C'_{ij}$ from the $C_{ij}$ using Eq. (31) applied to the $C_{ij}$. One gets

$$ 
\begin{align} 
C'_{11} &= C_{11} + 2T'^{1}_{3}C_{31} + T'^{1}_{3}T'^{1}_{3}C_{33} \nonumber \\
C'_{12} &= C_{12} + T'^{1}_{3}C_{32} + T'^{2}_{3}C_{13} + T'^{1}_{3}T'^{2}_{3}C_{33} \tag{34} \\
C'_{13} &= T'^{3}_{3}C_{13} + T'^{3}_{3}T'^{1}_{3}C_{33} \nonumber \\
C'_{33} &= T'^{3}_{3}T'^{3}_{3}C_{33}, \nonumber
\end{align} 
$$

where the $T'^{j}_{i}$ are taken from Eqs. (26) and the $C^{j}_{i}$ from Eqs. (33). In both cases, one gets

$$ 
\begin{align} 
C'_{11} &= C_{ii}^{b} + \frac{2}{H_{0}^{G}} (\gamma_{\infty}^{(1)(1)} + 2\gamma_{\infty}^{(1)}) - \frac{4\Delta_{r,\infty}^{(1)}}{H_{0}''} C_{i\neq j}^{b} - \frac{8\Delta_{r,\infty}^{(1)}}{H_{0}^{G}H_{0}''} (\gamma_{\infty}^{(1)(3)} + \gamma_{\infty}^{(3)}) \nonumber \\
&\quad + 4 \left(\frac{\Delta_{r,\infty}^{(1)}}{H_{0}''}\right)^{2} C_{ii}^{b} + \frac{8}{H_{0}^{G}} \left(\frac{\Delta_{r,\infty}^{(1)}}{H_{0}''}\right)^{2} \gamma_{\infty}^{(3,3)} \nonumber \\
C'_{12} &= C_{i\neq j}^{b} + \frac{2}{H_{0}^{G}} (\gamma_{\infty}^{(1)(2)} + \gamma_{\infty}^{(1)} + \gamma_{\infty}^{(2)} + \gamma_{\infty}) - \frac{4\Delta_{r,\infty}^{(1)}}{H_{0}''} C_{i\neq j}^{b} - \frac{8\Delta_{r,\infty}^{(1)}H_{0}'}{H_{0}^{G}H_{0}''^{2}} (\gamma_{\infty}^{(1)(3)} + \gamma_{\infty}^{(3)}) \tag{35} \\
&\quad + 4 \left(\frac{\Delta_{r,\infty}^{(1)}}{H_{0}''}\right)^{2} C_{ii}^{b} + \frac{8}{H_{0}^{G}} \left(\frac{\Delta_{r,\infty}^{(1)}}{H_{0}''}\right)^{2} \gamma_{\infty}^{(3,3)} \nonumber \\
C'_{13} &= \left\{ C_{i\neq j}^{b} + \frac{2}{H_{0}^{G}} (\gamma_{\infty}^{(1)(3)} + \gamma_{\infty}^{(3)}) - \frac{2\Delta_{r,\infty}^{(1)}}{H_{0}''} C_{ii}^{b} - \frac{4\Delta_{r,\infty}^{(1)}}{H_{0}^{G}H_{0}''} \gamma_{\infty}^{(3,3)} \right\} \frac{H_{0}'}{H_{0}''} \nonumber \\
C'_{33} &= \left\{ C_{ii}^{b} + \frac{2}{H_{0}^{G}} \gamma_{\infty}^{(3,3)} \right\} \left(\frac{H_{0}'}{H_{0}''}\right)^{2}, \nonumber
\end{align} 
$$

where $H_{0}^{G} = n\lambda_{b}$, $H_{0}' = (n-1)\lambda_{b} + 2\Delta_{r}$, and $H_{0}'' = (n-1)\lambda_{b} + 2\Delta_{r,\infty}^{(3)}$. Note that, even for the $C'_{ij}$, the relations are given as a function of the constants $\Delta_{r,\infty}^{(i)}$, $\gamma_{\infty}^{(i)}$ and $\gamma_{\infty}^{(i)(j)}$ given in Table II and not $\Delta_{r,n}^{(i)}$, $\gamma_{n}^{(i)}$ or $\gamma_{n}^{(i)(j)}$. The results for the predicted $C'_{ij}$ and $C_{ij}$ are shown as dashed and solid lines in Fig. 8 and Fig. 10 respectively. As can be seen, the agreement is quite good, even for extremely thin films (down to four layers).

The dependence of the elastic constants of the film with its thickness can now be discussed in the light of Eqs. (33) and 35. It is clear from Eqs. (33) that the intrinsic stiffness coefficients $C_{ij}$ are expected to exhibit a linear dependence with the inverse Gibbs thickness of the film $H_{0}^{G}$, or, equivalently, its number of layers, since $H_{0}^{G} = n\lambda_{b}$. This results from the property that $\gamma_{\infty}$, $\gamma_{\infty}^{(i)}$, and $\gamma_{\infty}^{(i)(j)}$ are independent of the film thickness. This extremely simple behavior is well confirmed by the molecular simulations, see Fig. 8. On the other hand, the evolution of the coefficients $C'_{ij}$ involve three thicknesses, $H_{0}^{G}$, $H_{0}'$, and $H_{0}''$. A fairly good approximation would be to replace $H_{0}'$ and $H_{0}''$ by $H_{0}^{G}$. In that case, Eqs. (35) exhibit a simpler form, with terms linear, quadratic and cubic
with $1/n$. Although the quadratic and cubic terms are relatively small, neglecting them would alter
the quality of the fitting with the data.

In order to appreciate the differences between the $C_{ij}$'s and the $C'_{ij}$'s, these two quantities have
been reported on the same graph. For clarity, the data have been omitted, and only the theoretical calculations, which perfectly fit the data without loss of information, have been reported (see
Fig. 14). As previously noticed, the $C_{ij}$ coefficients exhibit a linear dependence with the inverse film thickness $1/n$, while the $C'_{ij}$ are quasi linear. The diagonal coefficients $C_{ii}$ and $C'_{ii}$ all decrease for thinner films. On the other hand, the off-diagonal coefficients have different behaviors. The out-of-plane off-diagonal coefficients $C_{13}$ and $C'_{13}$ decrease, while the in-plane off-diagonal coefficients $C_{12}$ and $C'_{12}$ increase. In all cases, transforming the $C'_{ij}$ into the $C_{ij}$ results in slightly larger values, the difference being essentially the same for all coefficients, with increasing amplitude when the film is thinner (reaching 7 to 8 $\epsilon_{LJ}/\sigma_{LJ}^2$ for the thinnest film $n=4$). The same comparison can be done for the elastic moduli and the Poisson's ratios. As can be seen in Fig. 15, the moduli $E_1$ and $E'_1$ on one hand, and the Poisson's ratios $v_{12}$ and $v'_{12}$ on the other hand, are perfectly superimposed: the in-plane elastic constants are unaffected by the definition $\epsilon_3$ or $\epsilon'_3$ of the deformation along $z$. This is because these constants are obtained at zero pressure along $z$ (free surfaces), and does not require the knowledge of the deformation along $z$. Interestingly, the out-of-plane moduli $E_3$ and $E'_3$ are also almost identical, the difference being always smaller than 1%. On the other hand, the out-of-plane Poisson's ratios $v_{13}$ and $v'_{31}$ significantly differ from $v'_{13}$ and $v'_{31}$. Therefore, using the imposed stress method (not presented in this work, it would consist in imposing stresses on the film and measure deformations) to get directly, by simulations, the elastic moduli $E_1$ and $E_3$ of the film, would give identical results regardless of the choice of the definition of the deformation along $z$. The in-plane Poisson's ratio $v_{12}$ would also be correctly calculated. On the other hand, the remaining out-of-plane ratios $v_{13}$ and $v'_{31}$ require to pay attention to the definition of the deformation along $z$.

V. CONCLUSION

This work focuses on a rigorous definition and calculation of the elastic constants of a thin film
in the three space directions. Two issues are considered. The first one concerns the definition of
the deformation of the film in the out-of-plane direction, leading to two possible definitions $\epsilon_3$ and
$\epsilon'_3$, see section II E. The second concerns the interpretation of the stress tensor as given by the
simulations, which depends in particular on the protocol used to control the deformation along z.
In practice, simulations can easily monitor the position of the outermost crystallographic plane,
leading to $\sigma'_i$ values. A general formalism is developed, allowing anyone to switch between, on one
hand, the $\varepsilon'_i$ and $\sigma'_i$ as given by the simulations, and, on the other hand, the $\varepsilon_i$ and $\sigma_i$ as generally
considered in theoretical developments. This allows one to calculate the two corresponding sets
of elastic constants $C_{ij}$ and $C'_{ij}$.

In this framework, we have shown that the evolution, with the film thickness, of the intrinsic
constants $C_{ij}$, can be easily rationalized in terms of bulk + surface contributions, the latter being
derived from the surface free energy $\gamma$. It is shown that the simulation results can be reproduced
with an excellent accuracy for ultrathin films, down to only four crystallographic planes. It is also
shown that the corresponding Young's moduli $E_i$ and $E'_i$ are actually essentially identical, meaning
that they are insensitive to the definition $\varepsilon_3$ or $\varepsilon'_3$ of the film. The in-plane Poisson's ratio $\nu_{12}$ is
also identical to $\nu'_{12}$. Only the out-of-plane Poisson's ratio $\nu_{13}$ and $\nu_{31}$ deserve special attention
since their values depend on the strain definition $\varepsilon_3$ or $\varepsilon'_3$.

These results were obtained for a simple Lennard-Jones model of thin film. However, most
results are thought to be quite general, in particular the formalism used to switch between the two
systems ($\varepsilon_i$, $\sigma_i$) and ($\varepsilon'_i$, $\sigma'_i$), and the simple interpretation in terms of bulk + surface contributions
for the $C_{ij}$, leading to a linear dependence of the $C_{ij}$'s with the inverse film thickness. It is in
particular emphasized that, despite the fact that we have made the calculations for an initially
undeformed film, the formalism can be transposed to any initial state, at the expense of replacing
the absolute strains or stresses by their variations. This work has been done at zero temperature.
Of course, the temperature is expected to have quantitative effects, but, qualitatively, one should
observe similar behavior, in particular because the theoretical developments are based on the free
energy $\gamma$, and, therefore, valid for non-zero temperature. It is however mentioned that, at non-zero
temperature, the position of the crystallographic planes should be defined by the average position
of their atoms. Finally, shear has not been considered in this work, as well as a general treatment
of the crystallographic orientation of the free surfaces. This will be considered in future work, as
well as more realistic potentials to describe metallic, silicon or carbon films.
## ACKNOWLEDGMENTS

Authors are grateful to the CNRS interdisciplinary "Défi Needs" through its "MiPor" program (Project DARIUS) and the financial support of Agence Nationale de la Recherche through the project ANR-23-CE30-0028-03 (NanoCav).

## AUTHOR DECLARATIONS

### Conflict of interest

The authors have no conflicts to disclose.

### Author contributions

Patrice Porion: Conceptualization ; Data curation ; Writing & editing. Joël Puibasset: Conceptualization ; Data curation ; Writing & editing.

## DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

¹K. L. Ekinci and M. L. Roukes, Review of Scientific Instruments **76**, 061101 (2005).

²M. J. Gordon, T. Baron, F. Dhalluin, P. Gentile, and P. Ferret, Nano Letters **9**, 525 (2009).

³X. Li, T. Ono, Y. Wang, and M. Esashi, Applied Physics Letters **83**, 3081 (2003).

⁴S. H. Park, J. S. Kim, J. H. Park, J. S. Lee, Y. K. Choi, and O. M. Kwon, Thin Solid Films **492**, 285 (2005).

⁵H. Sadeghian, H. Goosen, A. Bossche, B. Thijsse, and F. van Keulen, Journal of Physics D: Applied Physics **44**, 072001 (2011).
⁶R. Janknecht, R. Hahn, N. Koutná, J. Todt, M. Meindlhumer, A. Davydok, H. Riedl, J. Keckes, and P. H. Mayrhofer, *Materials & Design* **252**, 113720 (2025).

⁷H. Sun, L. Chen, S. Sun, and T.-Y. Zhang, *Science China Technological Sciences* **61**, 687 (2018).

⁸J.-G. Guo and Y.-P. Zhao, *Journal of Applied Physics* **98**, 074306 (2005).

⁹R. C. Cammarata, *Progress in Surface Science* **46**, 1 (1994).

¹⁰A. I. Fedorchenko, A.-B. Wang, and H. H. Cheng, *Applied Physics Letters* **94**, 152111 (2009).

¹¹P. Gumbsch and M. S. Daw, *Physical Review B* **44**, 3934 (1991).

¹²W. Haiss, *Reports on Progress in Physics* **64**, 591 (2001).

¹³R. E. Miller and V. B. Shenoy, *Nanotechnology* **11**, 139 (2000).

¹⁴P. Müller and A. Saúl, *Surface Science Reports* **54**, 157 (2004).

¹⁵P. Müller, A. Saùl, and F. Leroy, *Advances in Natural Sciences: Nanoscience and Nanotechnology* **5**, 013002 (2014).

¹⁶K. Van Workum and J. J. de Pablo, *Physical Review E* **67**, 031601 (2003).

¹⁷D. Wolf, *Applied Physics Letters* **58**, 2081 (1991).

¹⁸L. G. Zhou and H. Huang, *Applied Physics Letters* **84**, 1940 (2004).

¹⁹T.-Y. Zhang, Z.-J. Wang, and W.-K. Chan, *Physical Review B* **81**, 195427 (2010).

²⁰T.-Y. Zhang, H. Ren, Z.-J. Wang, and S. Sun, *Acta Materialia* **59**, 4437 (2011).

²¹F. W. De Wette, L. H. Fowler, and B. R. A. Nijboer, *Physica* **54**, 292 (1971).

²²C. Domb and L. Salter, *Philosophical Magazine* **43**, 1083 (1952).

²³G. K. Horton, *American Journal of Physics* **36**, 93 (1968).

²⁴M. L. Klein and W. G. Hoover, *Physical Review B* **4**, 539 (1971).

²⁵P. Mohazzabi and F. Behroozi, *Physical Review B* **36**, 9820 (1987).

²⁶D. D. Phuong, N. T. Hoa, V. Van Hung, D. Q. Khoa, and H. K. Hieu, *The European Physical Journal B* **89**, 84 (2016).

²⁷K. V. Tretiakov and K. W. Wojciechowski, *Computer Physics Communications* **189**, 77 (2015).

²⁸K. Westera and E. R. Cowley, *Physical Review B* **11**, 4008 (1975).

²⁹V. B. Shenoy, *Physical Review B* **71**, 094104 (2005).

³⁰F. H. Streitz, R. C. Cammarata, and K. Sieradzki, *Physical Review B* **49**, 10699 (1994).

³¹C. T. Sun and H. Zhang, *Journal of Applied Physics* **93**, 1212 (2003).

³²D. Wolf, *Surface Science* **225**, 117 (1990).

³³G. Clavier, N. Desbiens, E. Bourasseau, V. Lachet, N. Brusselle-Dupend, and B. Rousseau, *Molecular Simulation* **43**, 1413 (2017).
<sup>34</sup>N. Miyazaki and Y. Shiozaki, JSME International Journal Series A Solid Mechanics and Material Engineering **39**, 606 (1996).

<sup>35</sup>D. J. Quesnel, D. S. Rimai, and L. P. DeMejo, Physical Review B **48**, 6795 (1993).

<sup>36</sup>J. R. Ray, Journal of Applied Physics **53**, 6441 (1982).

<sup>37</sup>J. R. Ray, Computer Physics Reports **8**, 109 (1988).

<sup>38</sup>M. Sprik, R. W. Impey, and M. L. Klein, Physical Review B **29**, 4368 (1984).

<sup>39</sup>D. R. Squire, A. C. Holt, and W. G. Hoover, Physica **42**, 388 (1969).

<sup>40</sup>K. L. Van Workum, *Mechanical properties of solids and nanoscopic structures via molecular simulations*, Ph. d. thesis (2002).

<sup>41</sup>J. Q. Broughton and G. H. Gilmer, The Journal of Chemical Physics **79**, 5105 (1983).

<sup>42</sup>R. C. Cammarata and K. Sieradzki, Physical Review Letters **62**, 2005 (1989).

<sup>43</sup>I. L. Chang, S.-H. Chang, and J.-C. Huang, International Journal of Solids and Structures **44**, 5818 (2007).

<sup>44</sup>J. C. Hamilton and W. G. Wolfer, Surface Science **603**, 1284 (2009).

<sup>45</sup>P. Lu, L. H. He, H. P. Lee, and C. Lu, International Journal of Solids and Structures **43**, 4631 (2006).

<sup>46</sup>F. H. Streitz, R. C. Cammarata, and K. Sieradzki, Physical Review B **49**, 10707 (1994).

<sup>47</sup>F. H. Streitz, K. Sieradzki, and R. C. Cammarata, Physical Review B **41**, 12285 (1990).

<sup>48</sup>G. W. Tęcza, Physical Review B **46**, 15447 (1992).

<sup>49</sup>F. Wang, Y. Fu, S. Xiao, B. H. Skallerud, S. Liu, S. Luo, H. Kristiansen, J. He, and Z. Zhang, ACS Applied Nano Materials **6**, 11491 (2023).

<sup>50</sup>D. Wolf and J. F. Lutsko, Physical Review Letters **60**, 1170 (1988).

<sup>51</sup>X. Zhou, H. Ren, B. Huang, and T. Zhang, Science in China Series E: Technological Sciences **57**, 680 (2014).

<sup>52</sup>R. Shuttleworth, Proceedings of the Physical Society A **63**, 444 (1950).

<sup>53</sup>L. D. Landau and E. M. Lifshitz, *Theory of elasticity* (Pergamon Press, London, 1959).

<sup>54</sup>E. R. Cowley, Physical Review B **28**, 3160 (1983).

<sup>55</sup>H. Heinz, R. A. Vaia, B. L. Farmer, and R. R. Naik, The Journal of Physical Chemistry C **112**, 17281 (2008).

<sup>56</sup>T. Iwaki, JSME International Journal Series A Solid Mechanics and Material Engineering **39**, 346 (1996).

<sup>57</sup>M. Li and W. L. Johnson, Physical Review B **46**, 5237 (1992).
⁵⁸D. J. Quesnel, D. S. Rimai, and L. P. DeMejo, Solid State Communications **85**, 171 (1993).

⁵⁹K. Kanhaiya, S. Kim, W. Im, and H. Heinz, npj Computational Materials **7**, 17 (2021).

⁶⁰M. S. Daw and M. I. Baskes, Physical Review B **29**, 6443 (1984).

⁶¹V. Rosato, M. Guillope, and B. Legrand, Philosophical Magazine A **59**, 321 (1989).

⁶²J. Tersoff, Physical Review B **39**, 5566 (1989).

⁶³A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in 't Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, Computer Physics Communications **271**, 108171 (2022).

⁶⁴J. Puibasset, The Journal of Physical Chemistry C **123**, 15027 (2019).

⁶⁵J. Li, B. Narsu, G. Yun, A. Wang, and Z. Gao, Applied Physics A **124**, 813 (2018).

⁶⁶A. Stukowski, Modelling and Simulation in Materials Science and Engineering **18**, 015012 (2010).
<figure><img src="image_4.png" /><figcaption>FIG. 1. Schematic representation of the thin film, using the OVITO$^{66}$ software, and showing the convention used for the axes parallel to the crystallographic directions. The spheres represent the atoms of the FCC structure. The color code is used to highlight the two innermost (blue) and the two outermost (green) crystallographic planes used to define $\epsilon_3$ and $\epsilon'_3$ respectively (see text). The remaining atoms are in red.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIG. 2. Right: side view of the thin film shown in Fig. 1. The circles represent important atomic layers with the same color code as in Fig. 1: the two innermost layers (blue), the two outermost layers (green) and two sub-surface layers (red).</figcaption></figure>

The atomic layers in-between are shown as dots. The horizontal dashed lines passing through the centers of the atoms materialize the crystallographic (001) planes parallel to the free surfaces. The film thickness $H'_\epsilon$ is defined as the distance between the two crystallographic planes defined by the outermost atoms (green circles). $\lambda_\epsilon$ is the distance between the two innermost atomic planes (blue circles). Left: ideal bulk crystal (infinite along $z$), built with a constant interplanar distance equal to $\lambda_\epsilon$. This crystal is used as the reference to measure the displacements of the atomic planes due to surface relaxation (in particular, $\Delta_r$ is the displacement of the outermost plane). $H_\epsilon^G$ is the distance between the two opposite Gibbs equimolar dividing surfaces of the film (dotted lines). These surfaces are positioned exactly midway between two adjacent atomic planes in the ideal crystal (left), and delimit exactly the same number of atoms as in the film (right).
<figure><img src="image_6.png" /><figcaption>FIG. 3. Schematic representation of the two mechanical tests used to determine the elastic coefficients of the film. In each case, a uniaxial strain ($\epsilon'_1$,0,0) or (0,0,$\epsilon'_3$) is imposed (red arrows), while the stresses $\sigma'_i$ are measured in the three directions. This set of measurements allows one to calculate the four independent stiffness coefficients $C'_{ij}$.</figcaption></figure>
<figure><img src="image_7.png" /></figure>

FIG. 4. Triangles: out-of-plane dimensionless stresses (i.e. divided by the unit stress $\varepsilon_{LJ}/\sigma_{LJ}^3$), for three film thicknesses $n = 14, 30$ and $50$ layers, deduced from the variations of the internal energy of the film with respect to $\varepsilon_3$ (giving $\sigma_3$, up-triangles) and with respect to $\varepsilon'_3$ (giving $\sigma'_3$, down-triangles); the lines are guides to the eye through the data. Circles: out-of-plane stress as given by the stress tensor calculated in simulations.

<figure><img src="image_8.png" /></figure>

FIG. 5. Dimensionless (i.e. divided by the unit length $\sigma_{LJ}$) surface relaxation displacements (counted positively outwards), of the crystallographic (001) planes labeled from the surface ($m = 1$ for the outermost plane). Three film thicknesses are shown, $n = 10$ (down triangles), 20 (circles) and 30 (up triangles). For $n = 10$, only five planes are shown by symmetry. Lines are exponential fits with Eq. (24). Insert: same data plotted in semi-log scale.
<figure><img src="image_9.png" /><figcaption>FIG. 6. Dimensionless (i.e. divided by the unit length $σ_{LJ}$) first derivatives of the amplitude $Δ_{r,n}$ of displacement of the outermost atomic layer (surface relaxation), relative to $ε'_i$ and $ε_i$, as obtained from Eqs. (17) and (18) respectively. The solid lines are single parameter best fits (see text).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIG. 7. Symbols: evolution, with the inverse film thickness $1/n$, of the coefficients $T_3^i$, given by Eqs. (14), and $T_3^{i'}$, given either by Eqs. (15) or by Eqs. (26) corresponding to the $(T^{-1})_3^i$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIG. 8. Symbols: dimensionless stiffness coefficients $C'_{ij}$ (i.e. $C'_{ij} \sigma_{LJ}^3 / \epsilon_{LJ}$) of the film, as obtained from the mechanical tests performed on the atomic LJ films, as a function of the inverse of its thickness. The symbols at $1/n = 0$ are those of the bulk crystal. Dashed lines: predictions from Eqs. (35).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>FIG. 9. Symbols: dimensionless Young's moduli $E'_i$ (i.e. $E'_i\sigma_{LJ}^3/\epsilon_{LJ}$) and Poisson's ratios $v'_{ij}$ deduced from the stiffness coefficients $C'_{ij}$ shown in Fig. 8. Dashed lines: predictions from Eqs. (35) after numerical transformation into Young's moduli and Poisson's ratios. The mean $(v'_{12} + v'_{31} + v'_{13})/3$ is also reported in the figure.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>FIG. 10. Symbols: Evolution with the inverse film thickness 1/n of the dimensionless intrinsic stiffness coefficients $C_{ij}$ (i.e. $C_{ij} \sigma_{LJ}^3 / \epsilon_{LJ}$) as obtained after transformation of the $C'_{ij}$, given by the simulations, with Eqs. (28). The symbols at 1/n = 0 are those of the bulk crystal. Lines: predictions from Eqs. (33).</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>FIG. 11. Symbols: dimensionless Young's moduli $E_i$ (i.e. $E_i\sigma_{LJ}^3/\epsilon_{LJ}$) and Poisson's ratios $v_{ij}$ deduced from the stiffness coefficients $C_{ij}$ shown in Fig. 10. Lines: predictions from Eqs. (33) after numerical transformation into Young's moduli and Poisson's ratios. The mean $(v_{12} + v_{31} + v_{13})/3$ is also reported in the figure.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>FIG. 12. Dimensionless surface tension $\gamma$ and its first derivatives $\gamma^{(i)}$ and $\gamma^{(i)}$ with respect to $\varepsilon_i'$ and $\varepsilon_i$, at zero deformation, and as a function of the inverse film thickness $1/n$. All quantities are displayed in reduced units, i.e. they are divided by the unit of surface tension $\varepsilon_{LJ}/\sigma_{LJ}^2$. In each case, all data superimpose. Horizontal lines are best fits with constant values.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>FIG. 13. Second derivatives of $\gamma$ at zero deformation, with respect to $\varepsilon'_i$ and $\varepsilon_i$, as a function of the inverse film thickness $1/n$. All quantities are displayed in reduced units, i.e. they are divided by the unit of surface tension $\varepsilon_{LJ}/\sigma_{LJ}^2$. Lines are best fits [see text and Eq. (32)].</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>FIG. 14. Solid and dashed lines: dimensionless $C_{ij}$ and $C'_{ij}$, as a function of the inverse film thickness $1/n$, predicted by Eqs. (33) and 35 and reproduced from Figs. 8 and 10 with the same color and line code.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>FIG. 15. Solid and dashed lines: dimensionless $E_i$ and $E'_i$, and Poisson's ratios $v_{ij}$ and $v'_{ij}$, as a function of the inverse film thickness $1/n$, predicted by Eqs. (33) and 35 after numerical matrix inversion, and reproduced from Figs. 9 and 11 with the same color and line codes.</figcaption></figure>