Finite Element Modeling of the Dielectric Response of Metal/Metal Oxide Nanocomposites: Coarse-Graining the Quantum Response

Nanocomposite materials with metallic or ceramic inclusions show great promise as highly-tunable functional materials, particularly for applications where high dielectric permittivities are desirable, such as charge-storage or energy-storage materials. These applications present a challenge for computational approaches, as the field response of the nanoscale inclusion is quantum in nature, yet any representative sample of the material must encompass hundreds if not thousands of atoms. As currently implemented, finite element methods offer some predictive power for macroscale matrix-inclusion composites. However, their applicability cannot necessarily be extended to few-nanometer, molecular scale inclusions, where quantum and interfacial effects gain importance in the overall response to the applied field. Here, we develop an adjustable finite element method approach to calculate the low frequency dielectric constant of composites consisting of a metal-oxide matrix with molecular-scale silver inclusions, by introducing an interfacial layer in the general model. A process for coarse-graining atomistic ab initio results to generate best fit finite element models is also laid out in this work. We show that a continuum model informed by ab initio results can capture many of the relevant polarization effects in a metal/metal oxide nanocomposite, at a fraction of the computational cost.


I. Introduction
Materials with high dielectric permittivity (high-κ dielectrics) are important research and development targets, as gate components in metal-oxide semiconductor transistors [1,2] and energy storage materials in high energy density capacitors [3,4]. To pursue industry demands of ever-increasing performance in these devices, capacitance has generally been increased by decreasing dielectric layer thickness. However, the thickness of traditional SiO 2 gate dielectric layers can not be decreased below roughly 1 nm, [2] since increasing leakage currents negate the improved capacitance. Using a thicker dielectric with a higher permittivity allows device footprints to continue shrinking while maintaining low leakage currents.
One method of producing high-permittivity dielectrics is to add highly polarizable inclusions to an existing lowerκ material [3,5]. As with other composites, the properties of the mixed material could be tunable. By adjusting the composition, size, and loading of inclusions, materials with desired properties such as static dielectric permittivity, breakdown voltage, and dielectric loss can be fabricated for a given application, while maintaining the favorable properties of the matrix, such as processability or interfacial fit with other device components.
Several useful theoretical methods have emerged to model and explain the properties of dielectric composites. Mean field and effective medium theories produce equations such as the Maxwell and Bruggeman models, based on the assumption of dielectric spheres in a homogeneous medium [15][16][17]. Similar analytical models have also been developed to model general ellipsoidal inclusions [18]. However, experimental results indicate that interfacial effects between the inclusions and matrix can cause significant deviations from model predictions [19][20][21], which researchers have attempted to address by including parameters for capturing the structure and polarizability of this interface [17,22].
In addition to to analytical models, the finite element and finite difference methods have provided significant insights into composite dielectric behavior for elongated and irregularly-shaped inclusions [23][24][25], threephase mixtures [24], and random distributions of in-clusions [26,27]. These models often include an explicit interfacial layer to account for the effects mentioned above [24,25]. Results from finite element methods have prompted modifications to effective medium models and helped elucidate the conditions under which effective medium theories lose their validity.
The applicability of effective medium theories is challenged in nanostructured materials, where the quantum nature of electronic polarization effects becomes important due to the small feature sizes. Finite element analysis is more broadly applicable: Although not directly addressing these quantum effects, the method accounts for interactions between inclusions and shape effects that become important for non-trivial loading. However, both methods rely heavily on parametrization from experimental results, and the unknown transferability of these parameters restricts their application for materials discovery in the nanoscale regime.
Ab initio electronic structure methods provide a way to explore the quantum behavior of nanocomposites (NCs) without a priori knowledge from experiments. However, these methods scale poorly with the number of electrons in the system, meaning they can only be applied to composites where the inclusions are sufficiently small to allow a unit cell on the order of a few hundred atoms. In a previous publication from our group (Hally et al. [28]), Car-Parrinello Molecular Dynamics (CPMD) simulations were used to study NC dielectrics by optimizing the electronic and nuclear coordinates under an applied electric field, with reasonable efficiency for an ab initio method. The method was applied to modelling the dielectric response of NCs with a metal oxide matrix and small noble metal NP inclusions [28]. The addition of the molecular silver NPs to the magnesium oxide matrix increased the theoretical permittivity by up to 30% for atomic loadings of 3.7% silver. This increase was attributed to both the high polarizability of these inclusions and an increased polarizability of matrix ions near the inclusion due to greater electronic delocalization -a quantum effect.
While elegant, this method is computationally expensive, requiring extended time access to supercomputing resources for larger than single-atom inclusions. In contrast, the finite element method (FEM) is a low-cost alternative that can be effectively utilized on a modern laptop. In this work, we use the method devised by Hally et al. to generate quantum mechanical data for silver inclusions in a magnesium oxide matrix and use it to parametrize a continuum model of these inclusions whose effective permittivity can be calculated via FEM. We investigate the model parameters that must be introduced to capture some of the quantum size and shape as well as interfacial effects of the molecular scale inclusions, and discuss the transferability of such a model to inclusion types in a similar chemical subspace. The model is then used to investigate the effects of the shape and loading of inclusions on the effective dielectric permittivity of Ag/MgO NCs.

II. Methods and Systems
A. Nanocomposite systems

Ag 8 /MgO Models
Throughout this work, we examine NCs composed of small, molecular silver inclusions within a magnesium oxide matrix. Using Car-Parrinello molecular dynamics and the modern theory of polarization, Hally et al. calculated the permittivity of bulk MgO and an Ag 8 /MgO NC to be 9.26 and 11.94, respectively, at an inclusion atomic loading of 3.7% [28]. Here, we apply a similar procedure to cubic MgO supercells of 160 and 512 atoms, with the Ag 8 inclusions, in order to simulate the % loading dependence of the permittivity. Detailed Information about these three systens can be found in Section S1.

Anisotropic inclusions
Two additional composites are considered to investigate shape effects: Ag 12 /MgO and Ag 18 /MgO, with rod-like and disk-like inclusions, respectively. In a real NC, metal clusters created by sputtering or synthetic pathways will have a certain degree of polydispersity, and individual inclusions may not be well approximated by spheres. Moreover, by way of experimental techniques that can finely control NP shape, rod and disk-like inclusions can be created to introduce anisotropic responses for certain applications. Here, we test the applicability of the model derived for spherical inclusions for the two anisotropic inclusions shown in Figure 1(a), together with the "spherical" Ag 8 inclusion.
The 12-atom rod was obtained by replacing 2 ion by 2 ion MgO sites in 3 adjacent layers of MgO, and optimizing the Ag-MgO supercell (atomic and lattice degrees of freedom) using CPMD. The 18-atom disk was created by replacing 3 ion by 3 ion MgO sites in two adjacent layers of MgO, then optimizing. In both cases, the supercells were generated to have 4-5 layers of MgO between NP images in any direction, giving a 288 atom supercell for the Ag 12 rod and a 384 atom supercell for Ag 18 disk. In (c), the atomic geometry is reduced to a dielectric ellipsoid within a dielectric matrix. An interfacial layer of thickness T int is added around the ellipsoid to capture interaction effects with the matrix ions. In the continuum model, the components are described by bulk permittivities, with the matrix, inclusion, and interface relative permittivities being ε mat , ε inc , and ε int , respectively.

B. Isolated Nanoparticles
We also applied the continuum model to calculate the static polarizabilities of gas phase silver NPs. Relaxed geometries of 4-to 55-atom Ag NPs were generated by annealing in SIESTA [29], with a PBE [30]/DZP approach and Troullier-Martins pseudopotentials obtained from the Siesta GGA Pseudopotential Database. Cluster polarizabilities α ii along Cartesian direction i were then calculated in SIESTA by applying a finite electric field E i of 0.0001 a.u. and calculating the change in NP dipole moment µ i relative to the zero-field state: Unless otherwise noted, the reported polarizabilities are averaged over the diagonal elements of the polarizability tensor: The response was assumed to be in the linear regime due to the small applied field [31], and interactions with image NPs were minimized by using a cubic cell with 100 Å side length. The electron delocalization volume for each NP was calculated as the volume of the 0.001 a.u. isosurface of the optimized electronic charge density, in line with previous work [32].

C. Car-Parrinello Molecular Dynamics
Car-Parrinello molecular dynamics were performed on NC structures using a plane wave basis set in the Quantum Espresso code to calculate their effective permittivities. Starting geometries were generated by placing inclusion atoms in matrix lattice sites. After an electronic minimization to reach the ground state, both the ionic coordinates and cell size were relaxed using damped dynamics to a force convergence threshold of 1.0 × e −4 a.u. Next, the vanishing field polarization was calculated according to the modern theory of polarization before applying a homogeneous electric field of 0.001 a.u. in the positive zdirection [33]. Finally, the ionic and electronic degrees of freedom were relaxed in the presence of the external field using damped dynamics until the polarization had converged. The difference between the relaxed ion dipole moment p relaxed-ion and vanishing field cell dipole moment p zero-field yields the static dielectric response. The change in polarization density ∆P of our NCs can then be written as: where Ω is the simulation cell volume. The static dielectric constant ε r of the NC is calculated from where E is the magnitude of the applied electric field.

D. Finite Element Analysis
When coarse-graining the NC systems, the system is considered to be a union of matrix and inclusion continuum materials, combined in the form of a representative volume element (RVE). A tetragonal RVE is generated [as shown Figure 1(c)], with side lengths equal to the side lengths of the optimized cell from the Car-Parrinello simulations. Then, the geometry of the inclusion is approximated as an ellipsoid with the same aspect ratio and volume of the atomistic inclusion, which are determined using Bader analysis as explained in Section II G. To take advantage of the cell symmetry, most results were obtained by simulating just one-eighth of the full RVE.
In order to compute the effective permittivity of the composite, the first step is to solve the electrostatic Laplace equation, which describes the electrostatic potential V under the assumption that the material is a neutral dielectric with volume charge density of 0.
Here, ε is the total dielectric permittivity, the product of the dielectric constant ε r and the permittivity of free space ε 0 . To obtain a particular solution for V , we consider the problem domain to be our RVE, and apply appropriate conditions on its boundaries, ∂ Ω. We apply Dirichlet boundary conditions on the top and bottom faces of the RVE, setting the value of the electrostatic potential for all points in space r on each of these faces.
Setting this condition models the composite as the dielectric layer in a parallel plate capacitor. In addition, we enforce the periodicity of the solution in the x and y directions with Neumann boundary conditions on the 4 sides of the RVE:n · ∇V (r) = 0 ∀r ∈ ∂ Ω side (8) wheren is the unit normal vector pointing outward from the boundary. Using the definition of the electric field E = −∇V , this boundary condition amounts to specifying that the electric field at the sides of the RVE is strictly parallel to those sides. The full set of boundary conditions model our NC as the dielectric of an infinite-area capacitor, with a potential difference of V 1 −V 0 applied across its plates.
Once the potential V is obtained from the finite element method, the energy U stored in the dielectric can be calculated as Here, ε(r) is the electric permittivity at each point in space r, which in our model will have different values for the elements comprising the matrix, interface, and inclusion. We can alternatively write the energy stored in a parallel-plate capacitor in terms of its capacitance C and its overall permittivity ε.
where C is the capacitance, A is the plate area, and d is the distance between the two plates. Solving for ε eff and then substituting the energy integral, we have an expression for the overall permittivity of a composite dielectric that can be calculated via numerical integration from the finite element solution.
All RVEs were created and meshed in Gmsh version 4.7.0 [34]. Meshes are first order order tetrahedral meshes optimized using the Netgen extension within Gmsh. Furthermore, meshes were refined locally on the surface of the inclusion based on the local curvature until suitable convergence was obtained, as described in Section S2. Finite element calculations were all performed using the finite element solver GetDP, version 3.3.1 [35].

E. Validation and Benchmarking
To validate our finite element approach, we computed the effective dielectric constant for a series and parallel composite geometry which have known analytical solutions. The serial case consists of a two-phase laminar composite with the electric field applied normal to the layers and has an effective dielectric constant given by [36] 1 where subscript m denotes the matrix, subscript i denotes the inclusion, and f is the volume fraction.
The parallel case consists of a two-phase composite where all interfaces between matrix and inclusion are strictly orthogonal to the applied field. This system has an effective dielectric constant given by [36] ε eff = f m ε m + f i ε i (14) For the series and parallel geometries, cubic meshes were generated with volume loadings of high-κ material between 5% and 50%. The matrix was assigned a dielectric constant of 1 and the inclusion a dielectric constant of 10, 20, and 30. In all cases, there was no mesh dis-cretization error, and the FEM results matched analytical predictions to within floating point error [See Figure S2].
For the parallel configuration, structures with cylindrical inclusions were also generated. The cylindrical case introduces error in the form of the discretization of the curved interface. The error in the resulting effective permittivity of the composite was found to be directly proportional to the volume discrepancy between an ideal cylinder and the discretized one. In addition, the discretization error was found to be well suppressed by specifying a sufficient number of nodes per radian of curvature in the faces of our geometry. Figure S3 shows the construction of these cylindrical geometries and an analysis of the error due to discretization of the curved surface.

F. Partitioning the Dielectric Response
Calculation of the effective permittivity of the composite in Equation (12) requires knowledge of the electric permittivity ε(r) at all points in the RVE. The bulk matrix response can be calculated computationally via CPMD or taken from experimental data. However, the response of the molecular nanoparticle inclusion and that of the interface layer need to be determined: the former is expected to differ significantly from either bulk metal of nanoparticle-in-vacuum values, whereas the latter is unknown. We parametrize these values in a coarse-graining procedure discussed in Sections III A and III C, using CPMD-calculated dielectric permittivities and the overall polarization, partitioned into inclusion and matrix contributions.
The response of the CPMD model can be partitioned using the polarization of Wannier centers within the composite as a heuristic mapping of the polarization strength of the inclusion, as distinct from the interface and matrix. Wannier centers calculated by Quantum Espresso were all found to lie very close to either oxygen or silver atom centers (the Mg were fully ionized and had no remaining valence electrons when using the pseudopotential approach with 10 core electrons). Thus, the centers could be unambiguously assigned to either the inclusion or to specific matrix ions. The total change in the dipole moment of the simulation cell, ∆p can be calculated as where ∆w i are the Wannier center displacements, 2e is the charge of the two electrons occupying each Wannier function, ∆R i are ionic displacements with the corresponding charges Z i . In all cases, the ∆ refers to the difference between the finite and vanishing electric field states.
For the inclusion specifically, the change in dipole moment can then be calculated as the sum of the ionic and electronic dipole moments over the silver ions alone and their associated Wannier Functions: [28] To obtain the corresponding dipole moments for the continuum model, the polarization density P as defined by classical electrostatics can be integrated over the desired model component. The projection of this polarization density onto the applied field direction i takes the following form at each point r in space: where E i is the total electric field along direction i.
Integrating Equation (17) over the entire simulation cell yields the total dipole moment, analogous to the sum over ions and Wannier centers in the atomistic picture. Therefore, to obtain a quantity comparable to the inclusion dipole moment defined in Equation (16), we integrate the polarization density over the inclusion.
Note that we consider only the component of p inc along the direction of applied field, and we add the symbol to distinguish this from the CPMD-calculated value in Equation 16. Unlike Equation 16, the ∆ is dropped here because, unlike in the modern theory of polarization, there is no dipole moment in the absence of an applied potential (no permanent dipole). This is because the continuum model used includes only neutral dielectrics without free or localized charges. Figure 9 depicts the partitioning of dipole moments among different components of the composite in both the continuum and atomistic models.
To directly compare the dipole moments from the atomistic and the continuum pictures, the applied field must be the same in both cases. However, in the continuum model, the linearity of the response in weak fields ensures that regardless of the applied field, the dipole moment of the inclusion always constitutes the same proportion of the total dipole moment. Therefore, to simplify the implementation, the applied potential in FEM was held constant, and both the CPMD and FEM inclusion dipole moments were normalized by dividing by the total cell dipole moment calculated by the respective method. Thus, ω, the error in the normalized dipole moment for the FEM inclusion is given These normalized dipole moments were then compared to help differentiate continuum models with corresponding effective permittivities. The model with both a low error in effective permittivity and the best agreement inclusion dipole moments was considered the "best-fit" continuum model.

G. Volume Loading Definition
To reasonably draw parallels between quantum calculations of dielectric response and the FEM -calculated permittivity, it is important to rigorously define the inclusion volume in the two methods. The % loading in ab initio simulations is defined in an atomistic fashion, and has to be directly compared to volume loadings in continuum FEM models. As inclusion atoms roughly replace MgO matrix atoms in a one-to-one ratio, one option is to consider the percent atomic loading as a stand-in for volume loading [28,37]. However, this approach tends to underestimate the true electron delocalization volume that can be assigned to the inclusion, as shown in Figure 2(a). Alternatively, we define the inclusion volume as the summed Bader volumes of the silver atoms [38], calculated using the Henkelman Group code [39]. The Bader volume of a given atom is defined by the zero-flux isosurface of that atom within the material. For the Ag NPs considered here, the resulting volume broadly envelops the full inclusion [see Figure 2(c)]. This treatment significantly improves the fit of a spherical inclusion model to the CPMD calculated results at low loadings, as shown in Figure 2(b), where different volume loadings of dielectric spheres with permittivities 5, 10, and 1000 times larger than the matrix permittivity (ε r = 9.26) were simulated using FEM. Notably, using a model of a metallic sphere does not significantly improve the fit to the CPMD data points relative to the most polarizable dielectric sphere shown at the volume loadings considered, and its curve would essentially be overlaid over the ε inc /ε mat = 1000 curve. Higher loadings in CPMD are not well treated with a single spherical inclusion model and require careful treatment of the local polarization of matrix near the inclusion, as discussed in Section III C. In panels (a) and (b), CPMD-calculated effective permittivities are plotted in black "x" marks against the atomic and volume loading percent of silver, respectively. FEM-calculated effective permittivities are shown for inclusion to matrix ratios ε i /ε m = 5, 20, and 1000 in blue, orange and green lines, respectively. The Bader volume of the Ag 8 inclusion is plotted in panel (c) for the 3.7% atomic loading / 6.8% volume loading case.

A. Spherical NPs and the Metal Sphere model
Several DFT-based studies have shown the static polarizability of silver and gold NPs to be directly proportional to the electron delocalization volume [31,32,40]. These analyses considered small and medium sized nanoparticles with up to several hundred atoms and is thus applicable 1-17 | 6 even to clusters with discrete energy levels. To confirm that the trend is also reflected by the DFT-estimated polarization response in the small molecular inclusions considered in our studies, we examined the polarizabilities of thirteen distinct silver NPs with 4 to 55 atoms, using the process described in Section II B. For small clusters up to N = 13, we found close agreement to the structures and polarizabilites calculated by Pereiro and Baldomir [31], who searched extensively for global minumum structures. Larger clusters determined through annealing are expected to be local minima, but the polarizability trends remain, showing the expected volume dependence. The resulting relationship between volume and polarizability of the NP in vacuum is plotted in Figure 3.
The proportionality between the NP polarizability and volume is similar to the behaviour of a metallic sphere with radius r, which has a polarizability of r 3 .
[41] For a dielectric sphere with dielectric constant ε r , the polarizability α is also directly proportional to its volume: [41][42][43] Note that in SI units, this expression would contain a factor of 4πε 0 , which we omit, reporting the polarizabilities in units of Bohr 3 . While NP polarizability increases linearly with volume, the polarizability per atom is found to, in general, decrease as NP size increases. This is due in part to the larger fraction of surface atoms in smaller clusters, and the distinct contribution from surface and core atoms to properties such as polarization [32,40,44]. The resulting dependence of polarizability on the number of atoms is plotted in Figure 3(a). In addition to the downward trend, there are large fluctuations in the polarizability per atom from cluster to cluster, with more elongated particles displaying larger polarizabilities than their higher-symmetry neighbors. We note that the elongated particles also tend to have a larger percentage of surface atoms. Thus, instead of considering only the number of atoms, or electrons, in a cluster, we consider the electronic delocalization volume, as discussed in Section II G. This volume depends on the electronic structure of the cluster [see Figure 3(b)], and in effect accounts for the distinction between surface (larger volume) and core (smaller volume) atoms.
We used a least squares procedure to fit the data in Figure 3(b) to Equation 20, assuming spherical volumes. The NPs were found to have a fitted dielectric constant of 12.92. While dependent on the precise definition of particle volume, this value seems to be a consistent characteristic of all particles surveyed when using our volume definition, yielding R 2 = 0.997. Notably, the vacuum polarizability of the series of molecular silver NPs falls short of that of a metallic sphere, represented by the gray line in Figure 3(b). In most cases, the electron delocalization volume of the optimized NPs was roughly spherical in shape, becoming more spherical as the number of atoms increased [see This analysis considers only the electronic polarizability, neglecting the ionic contribution that is present when these NPs are dispersed in a matrix. However, the results provide a good justification for a dielectric volume model for describing the average particle polarizability. The following sections evaluate the suitability of such a model for treating anisotropic gas-phase NPs and NPs embedded in a host matrix to determine the transferability of fitted dielectric constants and the importance of particle shape and volume in the polarization behavior of each.

B. Particle anisotropy
Not all NPs considered here are well described by the spherical model. In the dielectric sphere analysis of Figure 3, we smoothed over any NP anisotropy by considering the mean polarizability over the three principal axes of the unit cell. However, to capture the influence of nonsphericity on NP response (for rod-like or disk-like inclusions, for example), an anisotropic model can be derived by, for example, fitting an ellipsoidal model onto the volume of the atomistic NP. As discussed in Section II G, the proper volume of a cluster of atoms is not simply defined. Moreover, a Bader volume cannot be calculated for the NP in vacuum. In this case, a first approximation is to consider an oriented bounding box that encloses the electronic delocalization volume of the inclusion. The principal axes of the model ellipsoid are parallel and proportional to the oriented bounding box. This approach incorporates information about the aspect ratio and orientation of the particles, and accounts for most of the difference between the individual polarizability components calculated with DFT and the predictions of a spherical model (Figure 4).
Oriented bounding boxes calculated to enclose the 0.001 au electronic density isosurface were considered in Figure 4 for the set of thirteen NPs with 4 to 55 atoms discussed above. The bounding box side orientations were then taken to be the principal axes of an ellipsoid and scaled so that the volume of the ellipsoid was equal to that of the 0.001 au electronic density isosurface. Appendix A describes how the classical theory for a dielectric sphere in Equation 20 can then be modified to treat ellipsoids oriented arbitrarily relative to an applied field.
The dielectric constant of the ellipsoidal model was determined by fitting the Cartesian projections of the NP polarizability to the DFT data, for all NPs in Figure 3. A least squares fit of Equation (24) to the polarizability data plot- ted in Figure 4 was used to obtain ε XX = 13.13, ε YY = 13.08, and ε ZZ = 11.78. These values were averaged to obtain a dielectric constant of 12.72 for the ellipsoids, which is very close to the value calculated earlier for the spherical model (12.92). However, the ellipsoidal model produces a significantly better fit to the DFT data as the shape of the NP is better described by the more general model [see Figure  4 Previous work on nanoparticle polarizability in vacuum used a cylindrical jellium model to account for shape effects, with similar outcomes in reproducing DFT results in a classical model [45]. However, the ellipsoidal dielectric model where anisotropy is captured in the electron delocalization volume is readily incorporated into the FEM model for composite systems described in Section II D. The validity of this model was also tested on gold nanoparticles with clusters ranging from 4 to 55 atoms. As shown in Figure S10, the linear dependence on electronic delocalization volume holds well, and the ellipsoidal approximation is important for particles with larger aspect ratios.
To develop an FEM model for anisotropic inclusions, we first examined the validity of the dielectric ellipsoid model for NPs in vacuum for the Ag 12 and Ag 18 NPs [ Figure 1(a)]. The particles were removed from the MgO lattice and reoptimized in SIESTA using a conjugate gradient minimization, which resulted in increased bond lengths, but little overall reconfiguration. The polarizabilities of the particles along each Cartesian axis were then calculated as in Section III A for both the re-optimized particles and the particles as extracted from the MgO without re-optimization.
Fitted ellipsoidal approximations for each particle were initially obtained using the oriented bounding box procedure, and using the side lengths of the bounding box to set the ellipsoid aspect ratios. However, for the more elongated rod particles, this procedure produced ellipsoids that were too wide. Using a least squares fitting procedure to find the best-fit ellipsoid [46], followed by scaling the principal axes to attain an ellipsoid volume equal to the 0.001 au electron density isosurface volume of the NP [see Figure 4(b)] provided a better fit. The dielectric constant of the ellipsoidal dielectric inclusions was set to 12.72, the optimal value found for the ellipsoidal model, and the method outlined in Appendix A was used to calculate particle polarizabilities. Figure 5 shows the diagonal elements of the NP polarizability tensor in vacuum calculated using CPMD, and the two FEM inclusion models considered in this work: the spherical and the ellipsoidal particle approximations. On the whole, the ellipsoidal particle model reproduces anisotropic polarizability reasonably well, with better matching of the continuum and quantum models can be achieved when particle geometries were relaxed in vacuum. Because the unrelaxed particles had been compressed in the nanocomposite, the relaxation also resulted in increased particle volume, leading to a corresponding increase in polarizability. Figure 5 Anisotropic polarization of anisotropic NPs in vacuum. The diagonal elements of the polarizability tensor are shown, calculated with the electric field applied parallel to the corresponding axis, in the CPMD (columns) and FEM ellipsoidal (stars) and spherical (black dashed lines) dielectric models.

(d)
Results for NP structures identical to those of the inclusion in the optimized NC ("unrelaxed", UR) or followed by structural relaxation in vacuum ("relaxed", R) are shown. In the continuum models, the ellipsoidal particles had the fitted dielectric constant of 12.72, and the spherical particle model had a dielectric constant 12.92.

C. Modelling the Interfacial Layer
As shown in Figure 2, a single spherical inclusion model in FEM fails to accurately treat NC polarization at high volume loadings, regardless of the dielectric permittivity assigned to the inclusion. The polarization of the silver inclusion most significantly impacts the electron density on the nearby oxygen ions. This results in a stronger polarization within a shell of matrix ions around the inclusion. [28,37] There are two contributions to this effect: The local field enhancement due to the polarization of the inclusion, and a higher polarizability of the interfacial matrix layer, due to changes in the matrix electronic structure through its interaction with the polarized NP. In the Ag 8 /MgO system of Ref. [28], for a matrix shell of about 3 Å near the inclusion, coupling between the electrons in the mNP and those of the O anions leads to spread of the O-sp 3 orbitals and thus increased matrix polarizability. For larger anisotropic inclusions, the polarization enhancement was found to be mostly localized within the two layers of MgO ions closest to the inclusion [37].
To capture the localized coupling of the matrix and the inclusion polarization response within the continuum 1-17 | 9 model, we considered a distinct interfacial layer to describe this shell of higher polarization. The layer was assigned a permittivity value by fitting the continuum model to the CPMD results.

Fitting Procedure in Coarse-Graining the Quantum Response
Fitting of the FEM-derived NC permittivity to the CPMD-calculated one was done by a sweep over three parameters in the FEM model: The permittivity of the inclusion, that of the interfacial shell, and the interface thickness. Throughout the simulations, the inclusion volume was held constant at the Bader volume of the inclusion, as described in Section II G, and the permittivity of the matrix was held fixed at 9.26, the value obtained during CPMD simulations of the matrix alone. Figure S8 depicts the parameter search space along with the average errors in effective permittivity relative to CPMD values for the three systems shown in Figure 2. This error, τ, is calculated as τ = |ε CPMD −ε FEM | ε CPMD . From this initial analysis, an interface thickness of 1 Å was identified as producing almost universally lower errors than all other thicknesses [ Figure S8. Figure 6 depicts τ for the 1 Å interface thickness case, over different values of the inclusion and interface susceptibility, normalized to the atomic density n k (χ k = χ k n k , where k indicates either the interface or the inclusion). A region of this plot from the upper left to lower right shows relatively constant values of τ below 8%, which is within the error margins of the CPMD method with an LDA functional. This indicates that combinations with either (i) low inclusion susceptibility and high interface susceptibility or (ii) high inclusion susceptibility and low interface susceptibility parametrizations reproduce the CPMD effective permittivity equally well.
Different choices of parameters along the low permittivity error line lead to different behaviours in the FEM model, particularly when simulating anisotropic inclusions. An additional feature can be chosen in order to define a parameter appropriate for a particular model. In the case of spherical and close-to-spherical inclusions, the induced dipole moment of the inclusion can be used as a distinguishing property, in order to more closely reproduce the polarization behavior of the atomistic model. We calculate the optimal fit of the inclusion dipole by partitioning the polarization of the CPMD model using the method described in Section II F, and determining the relative error of the FEM inclusion dipole from Equation 19. This procedure was performed for the 216-atom supercell using polarization data from Hally [28] to obtain a single dipole moment error ω for each tile in the heatmap of Figure 6. A pictorial representation, as well as details of the partitioning scheme in the CPMD and FEM models are provided in the SI and Figure S9. Along the constant-τ curve [the region of black shading in Figure 6(a)], the point presenting the smallest discrepancy between the CPMD-and FEM-calculated inclusion dipoles (ω ∼ 1.0%) corresponds to very similar, high susceptibilities for both the inclusion and interface relative 1-17 | 10 to the matrix (χ inc ∼ 2000; χ int ∼ 2300; χ mat ∼ 100). The dipole error increases significantly in both directions from this minimum. These values indicate that (i) The spherical inclusion susceptibility is significantly larger when embedded in a matrix relative to its normalized susceptibility in vacuum (χ NP ∼ 500, see Section III A and Figure 3), and (ii) The interfacial MgO layer is much more polarizable than bulk MgO. The enhanced NP permittivity is largely due to an additional ionic polarizability component that is present in silver inclusions but absent in gas phase silver clusters: When embedded in MgO, the silver atoms develop nontrivial positive Born Effective Charges [28], and their reorganization in the presence of an electric field contributes to the overall polarization of the system. In contrast, neutral, spherical silver clusters only have an electronic component to their polarization in the gas phase.

D. Rod-and Disk-Like Inclusions
Finally, the dielectric ellipsoid model was applied to the NC model with particles embedded in an MgO matrix. Two orientations were considered for each inclusion, with the principal axis aligned either parallel (ε ) or perpendicular (ε ⊥ ) to the applied field. The ellipsoidal geometry of the inclusions in the continuum model was generated by finding the best-fit ellipsoid to the Bader zero flux surface of the inclusion from the CPMD-optimized NC electronic structure, where the ellipsoid volume was constrained to be equal to the Bader volume. Using the NP permittivity values derived in Section III B, we calculate the effective permittivities of the rod and disk NC materials using FEM as discussed in Section III C. The dielectric permittivities thus obtained are listed in Table I, as well as their CPMDcalculated counterparts. As implemented, the ellipsoidal model with a highly-polarizable 1 Å interface reproduces CPMD results for small rod-like and disks-like inclusions relatively well. Initial calculations on larger rods and disks, including those studied previously by this group [37], suggest that the model is more suitable to rod-like inclusions than disks, although further investigations in this respect are necessary.

IV. Conclusions
We derived a finite element model for the low-frequency dielectric response of nanocomposites of silver nanoparticles, described as polarizable dielectric spheres and ellipsoids, embedded in a dielectric matrix. When fit to a suitable training set, the ellipsoidal model faithfully reproduces the polarizability of gas phase particles, including the anisotropy due to shape effects. Extension to nanocomposite systems demonstrated the importance of a highlypolarizable interfacial layer in reconstructing the high effective permittivities calculated using CPMD. After fitting on spherical inclusions, the interfacial layer model was found to reproduce the effective permittivities of rod-and disk-like inclusions with reasonable fidelity (errors within 10%). Further work should be done to evaluate the potential of other non-ellipsoidal representations for particles embedded in matrix. For example, the flexibility of FEM allows modeling of other shapes such as rounded prisms, for which analytical models like effective medium approaches do not yet exist. For small, molecular-scale inclusions, the physical interpretation of the proposed interface model encountered some ambiguity as the definition of the interfacial shape and thickness caused discrete jumps in the number of enclosed atoms. We are currently refining the model to tie the observable properties of matrix atoms to those of the FEM model layers, and to create a systematic mapping of nanoparticle shapes onto FEM inclusion models.

V. Appendix A: Polarizability of Dielectric Ellipsoids in Vacuum
In addition to spheres, ellipsoids more generally also have a uniform polarization in a uniform external electric field. For a general ellipsoid, this uniform polarization density P j for an external field E ext applied along principal axis j is given in SI units as [41] P j = χε 0 1 + N j χ E ext (21) where N j is the depolarization factor along the axis j, given by and a i is the inclusion semi-axis length in direction i. Note that the depolarization factor obeys the relationship N x + N y +N z = 1, so that in the case of a sphere, N x = N y = N z = 1 3 . Multiplying this polarization density by the volume of the ellipsoid 4 3 πa x a y a z , and substituting χ = ε r − 1 gives an expression for the dipole moment at a given external field in terms of the ellipsoid dielectric constant (assuming isotropic permittivity).
Since polarizability α is defined by α = d E , the polarizability of the ellipsoidal inclusion along principal axis j is determined to be α = 4 3 πa x a y a z ε 0 (ε r − 1) 1 + N j (ε r − 1) One can verify that plugging in N = 1 3 yields the expression of the polarizability of a sphere given above in Equation 20.

Ellipsoid and External Field Not Aligned.
If the external field is not oriented along one of the ellipsoid axes, the polarization induced can be calculated as a superposition of the polarizations induced by the components of the external field that do lie along each ellipsoid axis. That is, if we call the ellipsoid axes j , then The three components P j give the polarization density P of the ellipsoid in the basis defined by its principal axes. Transforming these different components of polarization back to the global coordinate system can then be used to find the polarizability of the inclusion along the direction of applied field (written in the global basis). The ellipsoid principal axes a x , a y , and a z , can be written as normalized column vectorsâ x ,â y , andâ z in the global coordinate system. Together, these vectors form an orthonormal basis in real space. The matrix that changes from the ellipsoidal to the global coordinate system is A = â x â y â z (26) Thus, to write the electric field in the global coordinate system E from the electric field in the ellipsoidal coordinate system E , one would write And to write the electric field in the ellipsoidal coordinate system, one would write Writing the equation above for the polarization of the ellipsoid as a vector equation in the ellipsoidal coordinate system, we have where operations on the right-hand side are element-wise multiplication of the two vector quantities. The polarization in the global coordinate system is P = AP = P x P y P z T Assuming the external field was applied along an axis of the global coordinate system, say z, then the polarizability of the inclusion in the j th direction can be calculated as α jz = 4 3 πa x a y a z P j E z (32)