Dielectric metal/metal oxide nanocomposites: modeling response properties at multiple scales

Nanocomposites with metallic inclusions show great promise as tunable functional materials, particularly for applications where high permittivities are desirable, such as charge-storage. These applications strain quantum mechanical computational approaches, as any representative sample of the material includes hundreds if not thousands of atoms. Many continuum methods offer some predictive power for matrix-inclusion composites, but cannot be directly applied to composites with small inclusions, for which quantum and interfacial effects dominate. Here, we develop an adjustable finite element approach to calculate the permittivities of composites consisting of a metal-oxide matrix with nanometer-scale silver inclusions, by introducing an interfacial layer in the model. The approach involves solving the Laplace equation with Dirichlet and Neumann boundary conditions. We demonstrate that such a continuum model, when appropriately informed using quantum mechanical results, can capture


Introduction
Materials with high permittivities are important research and development targets, as gate components in metal-oxide semiconductor transistors [1,2] and high energy density capacitors [3,4]. Industry demands ever-increasing performance in these devices, with capacitance increased most often by decreasing dielectric layer thicknesses. However, the thicknesses of traditional SiO 2 gate layers cannot be decreased below ∼1 nm [2], because increasing leakage currents negate capacitance improvements. Using thicker dielectrics with larger permittivities allows device footprints to continue to shrink while maintaining low leakage currents.
One way to increase the permittivity of a dielectric is to add highly polarizable inclusions [3,5]. As with other composites, the properties of the mixed material are potentially tunable. By adjusting the composition, size, and loading of inclusions, materials with desired properties such as static dielectric permittivities, breakdown voltages, and dielectric losses can be fabricated. This can be accomplished while maintaining favorable properties of the matrix, such as processability and interfacial fit with other device components.
Several theoretical methods have been developed to model and explain the properties of dielectric composites. Effective medium theories such as the Maxwell-Garnett and Bruggeman models approximate the composite as dielectric spheres in a homogeneous medium [15][16][17]. The approximations lead to describing the composite's static dielectric response through a mixing formula of the component permittivities. Similar analytical models have been developed to treat general ellipsoidal inclusions [18]. However, experimental comparisons suggest that interfacial effects between the inclusions and matrix materials can cause significant deviations from these model predictions [19][20][21]. Researchers have attempted to address this issue by including parameters for capturing the interface structure and polarizability in the effective medium formalism [17,22].
Quantum mechanical electronic structure methods provide a means of exploring the behavior of nanocomposites without a priori knowledge. However, these methods scale poorly with the number of electrons, 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 two earlier publications from our group [23,24], Car-Parrinello quantum mechanical molecular dynamics (quantum dynamics) simulations were used to study nanocomposite dielectrics by optimizing electronic and nuclear coordinates under an applied electric field. The method was applied to modeling the dielectric response of metal-oxide matrix materials and noble metal nanoparticle inclusions [23]. The addition of 8-atom silver nanoparticles to the magnesium oxide matrix increased the permittivity by up to 30% for atomic loadings of 3.7% Ag. Larger, anisotropic inclusions led to further increases in permittivity, at similar loadings. These increases were attributed to both the high polarizabilities of the inclusions, and the increased polarizabilities of matrix ions near the inclusions, due to greater electron delocalization-a quantum effect. While elegant, this method is computationally expensive.
Among classical continuum approaches, the finite element method (FEM) and the finite difference method have been previously considered for describing the polarization response of mixed materials. Previous works have provided significant insight into the dielectric behavior of mixtures including elongated and irregularly-shaped microscale inclusions, [25][26][27] three-phase mixtures [26], and random distributions of inclusions [28,29]. These models often include an explicit interfacial layer, to account for the effects mentioned above [26,27]. The results of these numerical calculations have prompted modifications to effective medium models, and helped elucidate the conditions under which those models lose their validity.
The applicability of FEMs to nanostructured materials has not been widely explored. In nanostructured systems, quantum effects due to small feature size and the electronic coupling between the inclusion and matrix can dominate the polarization response of the material [23,24]. Although not directly addressing quantum effects, the FEM can be parameterized for these interactions as well as the shape effects that become important for non-trivial loading.
In the sections that follow, we explore such a parameterized FEM model (a term we use interchangeably with continuum model), based on quantum-level computational results. The set of parameters necessary to capture inclusion quantum mechanical size, shape and interfacial effects is considered in detail. We assess the transferability of these parameters to a range of inclusion types for mixtures of MgO and nanoscale Ag inclusions. Along the way, we propose a definition of the continuum geometry based on Bader volumes calculated from quantum dynamics, define a procedure for fitting FEM parameters to quantum results, and validate the use of a continuum approach to the polarizability of isolated silver nanoparticles to motivate a similar approach for nanocomposites. Herein, the method devised by Hally and Paci [23] is used to generate quantum mechanical data for silver inclusions in MgO matrices. Based on this data, we parameterize a continuum model whose permittivity, ε, can be calculated via the FEM. The model is then used to investigate the effects of the shape and loading of inclusions on the permittivities of Ag/MgO nanocomposites [30].

Methods and systems
2.1. Nanocomposite systems 2.1.1. Models with 'spherical' Ag 8 inclusions. The methodology developed in this work is applied to a series of MgO nanocomposites with Ag inclusions. Using Car-Parrinello quantum mechanical molecular dynamics and the modern theory of polarization, Hally and Paci calculated the relative permittivity of bulk MgO to be ε r = 9.3 [23]. ε r = ε/ε 0 , where ε 0 is the permittivity of free space. The relative permittivity of an Ag 8 /MgO nanocomposite was also determined, and was found to be 11.9, for an inclusion atomic loading of 3.7% [23]. To explore the loading dependence of the permittivity, we examined two additional nanocomposites, comprised of MgO supercells of 64 and 512 atoms and Ag 8 inclusions. Additional details for these systems are provided in section S1.

Anisotropic inclusions.
Nonspherical inclusions can provide desirable material properties, and various nanoparticle synthesis techniques can be used to finely control particle shape. Thus, there is the potential to use rod-and disk-like inclusions to introduce an anisotropic dielectric response. With this in mind, we considered five composites with rod-and disk-like inclusions from [24] to investigate shape effects. All geometries were created by substitution of Ag at adjacent MgO lattice sites followed by optimization of the lattice and atomic coordinates. Hereafter, we denote these inclusions by the number of atoms substituted in each direction (2 × 2 × N for rods, and N × N × 2 for disks). These composites are described in more detail in section S1.2, and the smallest two anisotropic inclusions are shown in figure 1(d), alongside 'spherical' Ag 8 .

Isolated nanoparticles
To gain insight into the different contributions to the material's polarization response, we examine below a continuum model of the polarization of isolated nanoparticles. Quantum data used for the development of the model was generated as follows. Relaxed nanoparticle geometries were obtained using density functional theory (DFT) by annealing in SIESTA [31], using a PBE [32]/DZP approach and Troullier-Martins pseudopotentials [33,34] obtained from the SIESTA Database. Polarizabilities α ii along Cartesian directions i were then calculated by applying an electric field E i of 0.0001 a.u. (0.005 eV Å −1 ), and calculating the change in dipole moment µ i relative to the zero-field states: Unless otherwise noted, the reported polarizabilities are averaged over the diagonal elements of the polarizability tensor: The response was assumed to be linear due to the small size of the applied field [35], and interactions with image particles were minimized by using a cubic cell with 100 Å side lengths. The electron delocalization volumes were calculated as the volumes of the 0.001 a.u. isosurfaces of the optimized electronic charge densities, as in previous work [36].

Car-Parrinello molecular dynamics
Car-Parrinello quantum mechanical molecular dynamics were performed using plane wave basis sets in Quantum Espresso, version 6.6. After electronic minimization, both the ionic coordinates and cell size were relaxed using damped dynamics, with a force convergence threshold of 1.0 × e −4 a.u. Next, the vanishing field polarization was calculated before applying a homogeneous electric field of 0.001 a.u [37]. Finally, the ionic and electronic degrees of freedom were relaxed using damping, until the polarization had converged. The relative permittivity is calculated from where E is the magnitude of the applied electric field, ∆p is the change in simulation cell dipole moment relative to the vanishing field case, and Ω is the cell volume. Car-Parrinello molecular dynamics was also used to investigate structural coupling of interfacial ions to inclusions in nanocomposites. All structures were equilibrated to 400 K using a Nosé-Hoover thermostat for a time of 4 × 10 4 a.u. before the thermostat was removed and the NVE ensemble sampled for a time of another 8 × 10 4 a.u. Mean square displacements of ions were calculated for the NVE ensemble using the python wrapper of the analisi code, version 0.5.0 [38].

Continuum model and the FEM
In the continuum approach, a composite system is considered to be a union of matrix and inclusion materials, combined in the form of a representative volume element (RVE). Orthorhombic RVEs were generated (see figures 1(b) and (c)), using optimized cell lengths from quantum dynamics simulations. The geometries of the inclusions were approximated by ellipsoids, with the same aspect ratios and volumes as the quantum inclusions. The ratios and volumes were determined using Bader analysis, as explained in section 2.6. Symmetry allowed the simulation of just one-eighth of the full RVE (figure 1(c)).
There are multiple options in choosing a shape for the nanoparticle inclusions in the FEM model: nanoparticles can be described by meshing the full Bader volume, or geometric shapes such as prisms, cylinders or ellipsoids may be used. In the following pages, we chose to implement the theory using ellipsoidal shapes for several reasons: (1) there are simple analytical expressions for their polarizability in vacuum that facilitate analysis of isolated particles. (2) They are commonly used as models for anisotropic inclusions in effective medium theory [16]. (3) There are available procedures for fitting ellipsoids to arbitrary shapes, such as Bader surfaces [39]. (4) It is simple to scale and extrapolate ellipsoidal shapes for larger inclusions that are not simulated quantum mechanically. Using Bader surfaces themselves as inclusion shapes, in contrast, would require a quantum simulation of each inclusion to be simulated prior to performing the FEM.
Nevertheless, the FEM is flexible. Since changing inclusion shape does not affect the boundary conditions employed but only the complexity of the mesh used, arbitrary shapes could be used to model inclusions where appropriate. However, the fitting for inclusion and interfacial permittivity must be done separately for each type of shape used. Aside from ellipsoids, we discuss analysis done with rounded prism inclusions in section S9.
In order to compute the permittivity, a uniform potential was simulated across the RVEs, and the electrostatic Laplace equation was solved using the methodology outlined in section S2 as in previous studies [26][27][28]. The RVEs were created and meshed in Gmsh version 4.7.0 [40]. First order tetrahedral meshes were used, and they were optimized using the Netgen extension. Meshes were further refined on the surfaces of the inclusions, based on local curvature and until suitable convergence was obtained, as described in section S3. Finite element calculations were performed using the GetDP solver, version 3.3.1 [41].

Partitioning the dielectric response
The calculation of the permittivity of a composite via the FEM requires knowledge of the permittivity ε(r) at all points in the RVE. The bulk matrix response can be calculated via quantum dynamics or obtained by experiment. However, the response of the inclusions and that of the interface layers needs to be determined: the former may differ significantly from either bulk metal or nanoparticle-in-vacuum values, whereas the latter is completely unknown. We parameterize these values in a coarse-graining procedure discussed in sections 3.1 and 3.2.2, using quantum dynamics-calculated permittivities and the overall polarization, partitioned into inclusion and matrix contributions.
The response of a composite modeled using quantum dynamics 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 (see section S5.1 and [23] for a full description). Hally and Paci found that for the 216-atom Ag 8 nanocomposite under applied electric fields, the change in dipole moment relative to zero-field of the inclusion was 1.221 a.u., compared to a total change in the simulation cell dipole moment of 12.688 a.u. (calculated using equation (3) and the values reported for ε r , E, and Ω in [23]). Thus, the ratio of the inclusion dipole moment to the total cell dipole moment, which we call P r,quant , is 0.0926 for the 216-atom Ag 8 nanocomposite. We used this ratio as a point of comparison with the continuum model of the same 216-atom cell.
To obtain the corresponding dipole moment ratio for the continuum model, the polarization density P as defined by classical electrostatics can be integrated over the desired model component when the external field is applied. The projection of this density onto the applied field direction i takes the following form: Integrating equation (4) over the simulation cell yields the cell dipole moment. Similarly, the inclusion dipole moment is obtained by integrating the polarization density over the inclusion, Note that only the component of p inc along the direction of the applied field is considered. Figure S8 depicts the partitioning of dipole moments among different components of the composite in both the continuum and quantum models.
We call the ratio of the inclusion dipole moment p inc to the total cell moment in the continuum model P r,FEM . To directly compare the partitioning of the polarization between quantum and continuum models, we calculated the relative difference between P r,quant and P r,FEM , denoted ω, The continuum model with both a low permittivity error and the lowest ω was deemed the 'best-fit' model.

Volume loading definition
To draw parallels between quantum mechanical and FEM calculations of permittivity, it is necessary to rigorously define the inclusion volumes. The loading in quantum simulations is defined by atom volumes. Inclusion atoms replace matrix atoms in a one-to-one ratio. This suggests one option, which is to consider the percent atomic loading as a stand-in for volume loading [23,24]. While this is a sufficient definition when loading is used only as a system descriptor, we find that a more rigorous definition of volume is necessary when developing inclusion shape, size and loading parameters for FEM. A more precise description of electron delocalization volume in mixed materials requires the partition of electron density as done in Bader volume calculations [42]. Here, we define the inclusion volume as the summed Bader volumes of the inclusion atoms, calculated using the Henkelman Group code [43].
This total volume envelopes the full inclusion (see figure 2(a)). Furthermore, this treatment significantly improves the fit of a spherical inclusion model to the quantum results at low loadings, as shown in figure 2(b), for which different volume loadings of spheres with permittivities 5, 10, and 1000 times larger than the matrix permittivity (ε r = 9.3) were simulated using the FEM. Using a metallic sphere model does not significantly improve the fit to the quantum data relative to the most polarizable dielectric sphere shown at the volume loadings considered, and its curve would essentially overlay the ε inc /ε mat = 1000 curve. Higher loadings in quantum dynamics 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 3.2.2. Discrepancy between atomic and volume loadings of silver. The Bader volume of a Ag 8 inclusion is plotted in panel (a), for the 3.7% atomic loading/6.8% volume loading case. In panel (b), quantum dynamics-calculated permittivities are plotted using red diamonds and black triangles against the atomic and volume loading percent of silver, respectively. FEM-calculated permittivities are shown for different inclusion/matrix permittivity ratios. Several DFT-based studies have shown the static polarizability α of silver and gold nanoparticles to be directly proportional to the electron delocalization volume [35,36,44]. These analyses considered particles with up to several hundred atoms, and they are applicable to particles that include those with discrete energy levels. To confirm that this trend also holds for the responses considered in our study, we examined the polarizabilities of thirteen silver particles containing from four to 55 atoms, using the process described in section 2.2. For isolated silver nanoparticles of up to thirteen atoms, we found close agreement to the minimum energy structures and polarizabilites calculated by Pereiro and Baldomir [35]. Structures for larger nanoparticles were determined through a single annealing event and are expected to be local minima. The polarizability trends remain throughout the size range examined here, showing the expected linear volume dependence. The resulting relationship between volume and α is plotted in figure 3(a). The proportionality between particle polarizability and volume is similar to that of metallic spheres [45]. For a sphere of radius r, α is also directly proportional to its volume: [45][46][47]

Results and discussion
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 nanoparticle polarizability increases linearly with volume, the polarizability per atom was found to, in general, decrease as particle size increases, as shown in figure 3(b). This is due in part to the larger fraction of surface atoms in smaller particles, and the distinct contribution from surface and core atoms to properties such as α [36,44,48]. In addition to the downward trend, there are large fluctuations from particle to particle, 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 particle, we consider the electron delocalization volume. This volume depends on the electronic structure of the particle (see figure 3(a)), and 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(a) to equation (7), assuming spherical volumes. The fit assigned the particles a relative permittivity of 12.9. While dependent on the precise definition of particle volume, this value is a consistent characteristic of all particles surveyed when using our volume definition, yielding R 2 = 0.997. Note, the polarizabilities of the particles are less than those of analogous metallic spheres, which are represented by the gray line in the figure panel. In most cases, the particle volume was roughly spherical, becoming more spherical as the number of atoms increased. The change in electron density with the application of an electric field closely resembles the polarization charge induced on a sphere, or in a coarser approximation, the charge separation of a point dipole ( figure 3(c)).
This analysis considers only the electronic polarizability, neglecting the ionic contribution that is present when particles are dispersed in a matrix. Still, the results provide a clear justification for using the volume model for describing polarizability. The following sections evaluate the suitability of such a model for treating anisotropic gas-phase particles and particles embedded in host matrices to determine the transferability of fitted permittivities, and the importance of particle shapes and volumes in the polarization behavior of each.

Particle anisotropy.
Not all particles are well approximated as spheres. In the analysis of figure 3, we smoothed over anisotropy by considering the mean polarizability over the three principal axes. One way to capture the influence of non-sphericity is to fit an ellipsoidal model. As discussed in section 2.6, defining the volume of a particle is nontrivial. Moreover, a Bader volume cannot be calculated for a particle in vacuum. One way to proceed is to consider a bounding box that encloses the electron delocalization volume of the inclusion. The principal axes of the ellipsoid are defined parallel in directions and proportional in lengths to the box. In this way, this approach incorporates aspect ratio and particle orientation information. Boxes were calculated to enclose the 0.001 a.u. electron density isosurfaces for the thirteen particles discussed in section 3.1. The box sides were scaled so that the volumes of the ellipsoids were equal to those of the 0.001 a.u. isosurfaces.
The classical theory for a dielectric sphere can then be modified to treat ellipsoids oriented arbitrarily relative to an applied field, as described in section S6. This model accounts for most of the differences between the polarizability components calculated with quantum mechanics and the predictions of an analogous spherical model, as shown in figure 4. The permittivity of the ellipsoidal model was determined by fitting the Cartesian projections of the particle polarizabilities to the quantum mechanical data. A least squares fit of equation S(16) to the polarizability data was used to obtain the relative permittivities ε XX = 13.13, ε YY = 13.08, and ε ZZ = 11.78, giving an average of 12.7, which is very close to that of the spherical model (12.9). Nevertheless, the ellipsoidal model produces a significantly better fit to the quantum data, because shape and orientation are included in the model (see figure 4(b)). This indicates the importance of shape and orientation and provides justification for treating rod-and disklike inclusions with an ellipsoidal model (section 3.2.4).
Previous work on nanoparticle polarizability in vacuum used a classical mechanical cylindrical jellium model to account for shape effects, and showed a similar ability to reproduce quantum mechanical results [49]. The ellipsoidal model has the advantage that it is readily incorporated into FEM models for composite systems. To develop such a FEM model, we began by testing the ellipsoid model on two rod-and disk-like particles from [24], Ag 12 (2 × 2 × 3) and Ag 18 (3×3×2) (see figure 1(d)). Particles were re-optimized in vaccum using PBE/DZP using 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 3.1, for both the re-optimized particles and the geometries without re-optimization. Ellipsoidal approximations for each particle were initially obtained using the bounding box procedure described above. However, for Ag 12 , this produced an ellipsoid that was too wide. Using a least squares procedure [39], followed by scaling the principal axes to attain an ellipsoid volume equal to the 0.001 a.u. isosurface volume (see figure 4(b)) provided a better fit. The relative permittivity of the inclusions was set to 12.7, and the method outlined in section S6 was used to calculate polarizabilities. Figure 5 shows the diagonal elements of the polarizability tensor, calculated using quantum dynamics and the two classical models. The spherical model reproduces the average polarizabilities with reasonable accuracy. The ellipsoidal model reproduces the quantum results well, including capturing the subtleties associated with anisotropy.

Quantum modeling. Previous quantum mechanical investigations from our group of
MgO/Ag nanocomposites have elucidated several mechanisms by which inclusions enhance permittivity. First, Born effective charge analysis has shown that the large electronic polarizability of the embedded nanoparticles couples to both the electronic and ionic polarization of surrounding matrix atoms [23,24]. Second, Hally and Paci [23] explored the specific importance It is important to ascertain how a continuum model for these nanocomposites can capture these effects. The first phenomenon, that the polarizable nanoparticle generates an induced field which increases polarization in the surrounding matrix, is already qualitatively present in any continuum model with a dielectric ellipsoid with permittivity ε ell in a polarized dielectric matrix with permittivity ε mat , when ε ell > ε mat [50]. Nevertheless, there will be a quantitative difference in the screening of this induced field by the surrounding ions in the atomic model and by a continuum dielectric. In addition, figure 2 already shows that the continuum model of a dielectric sphere in matrix cannot reproduce the high polarization exhibited in quantum simulations. We believe this shortcoming is in large part explained by the greater electronic delocalization and coupling of the ionic dynamics of O and Ag at the interface. These quantum phenomena result in a layer of MgO atoms neighboring Ag that behave distinctly from bulk and are specifically more polarizable than bulk.
A distinct and more polarizable interfacial layer of MgO is supported by two additional analyses conducted in this work. First, Bader charge analysis (figure S1) shows a small transfer of electronic density from interfacial oxygen atoms to the inclusion, with a rapid return to bulk charges for subsequent layers. Second, we found significant increases in the mobilities of Mg and O ions neighboring Ag compared to bulk, furthering Hally's notion of ionic coupling. Mean squared displacements (MSDs) of ions over the course of Car-Parrinello molecular dynamics are shown in figure 6. The MSD for bare MgO shows comparable displacements of Mg and O ions, while nanocomposite MSD plots show that Ag has larger mobility within the unit cell. Furthermore, plots of the MSD of interfacial Mg and O ions, specifically those directly neighboring Ag atoms, show significantly larger displacement relative to other matrix ions. Thus, ionic mobility is enhanced at the interface, and this enhancement is localized.
This structural coupling manifests as a large enhancement of the ionic polarization of nanocomposites relative to MgO bulk. Table S1 shows that the enhancement of relaxed ion permittivity for nanocomposites relative to MgO is generally larger than the enhancement of electronic permittivity. This discrepancy is exaggerated as the percent loading of Ag increases, in which case the percent of interfacial Mg and O ions also increases, pointing to the importance of interfacial ionic coupling.

Continuum modeling of the interfacial layer.
Therefore, a continuum model must account for both the electronic coupling to the Ag inclusion and the increased ionic polarizability of the ions nearest the inclusion. To capture this coupling, we include an explicit interfacial layer, which is also ellipsoidal in shape, but with its axes all increased by a fixed thickness relative to the inclusion. This interface is assigned its own permittivity different than that of the bulk matrix by fitting the continuum model to quantum dynamics results as detailed in the following section. As expected, the best fit model includes a high-permittivity interface.

3.2.3.
Fitting procedure in coarse-graining the quantum response. Fitting of continuum to quantum dynamics permittivity was performed using a sweep over three parameters: the permittivities of the inclusion, that of the interfacial shell, and the interface thickness. The inclusion volume and matrix relative permittivity were fixed at the Bader value and 9.3, respectively. A 1 Å interface was selected to avoid encroachment on cell boundaries for high volume loadings (section S3.5). More importantly, this interface thickness also resulted in lower average errors than 2 or 3 Å interfaces, or no interface. Figure 7 depicts the parameter search space as well as the average permittivity errors τ for the three systems of figure 2 is the relative error. A region of this plot from the lower left to upper right shows relatively constant values of τ below 8%, which is within the error margins of the quantum dynamics method. This indicates that combinations with either (i) low inclusion and high interface permittivities or (ii) high inclusion and low interface permittivities, reproduced the quantum dynamics permittivity equally well. The induced dipole moment was used as an additional parameter in defining our FEM model. We calculated the optimal fit of the inclusion dipole by partitioning the polarization of the quantum dynamics model using the method described in section 2.5, and determined the relative error of the FEM inclusion dipole (ω) from equation (6). This procedure was performed for the 216-atom supercell using data from Hally and Paci [23] to obtain an ω for each tile in the heatmap of figure 7. A pictorial representation, as well as details of the scheme in the quantum dynamics and FEM models are provided in section S5 and figure S8.
Along the constant-τ curve (the region of black shading in figure 7(a)), the area of smallest discrepancy between the quantum dynamics-and FEM-calculated inclusion dipoles (ω ∼ 1.0%, indicated by the smallest, yellow circles) corresponds to high permittivities for both the inclusion and interface relative to the matrix (ε inc = 100; ε int = 300; ε mat = 9.3). The dipole error increases significantly in both directions from this minimum. These values indicate that (i) the inclusion permittivity is significantly larger when embedded in a matrix relative to its permittivity in vacuum (ε NP ∼ 13, see section 3.1 and figure 3), and (ii) the interfacial MgO layer is much more polarizable than bulk MgO. The enhanced nanoparticle permittivity is largely due to an additional ionic polarizability component that is present in silver inclusions but absent in gas phase silver particles: when embedded in MgO, the silver atoms develop nontrivial positive Born effective charges [23], and the reorganization of these silvers and the interfacial MgO ions in the presence of an electric field contributes to the overall polarization of the system. In contrast, spherical silver particles in the gas phase only have an electronic component to their polarization.
While it may seem strange for the interface permittivity to be larger than the inclusion permittivity, this is an artifact of representing discrete atomic layers as a continuum. With the selected thickness of 1.0 Å and inner radius of 3.27 Å, the interface contains 24 atoms and an atomic density of 0.134 atoms Å −3 . The inclusion, however, contains only eight silver atoms and a density of 0.055 atoms Å −3 . Thus the high interface permittivity is largely explained by high atomic density in the continuum interface region. This effect could be mitigated by enforcing an interface shape that is consistent with the matrix crystal system (cubic here) and by constraining the interface thickness to be integer multiples of the matrix lattice spacing.

Rod-and disk-like inclusions.
The ellipsoid model was applied to the series of rodand disk-like particles embedded in a MgO matrix from [24]. Two orientations were considered for each inclusion, with the long principal axes aligned either parallel (ε ∥ ) or perpendicular (ε ⊥ ) to the applied field. The continuum geometries of the 2 × 2 × 3 and 3 × 3 × 2 inclusions were generated using the bounding box method from section 3.1.2. Ellipsoids for the remaining inclusions were extrapolated using the procedure in section S10. Using the particle permittivity values derived in section 3.2.3, we calculate the permittivities of composites with rod-and disk-shaped inclusions using the FEM. The relative permittivities thus obtained, as well as their quantum-dynamics-calculated counterparts, are listed in table 1.
The ellipsoidal model displayed mixed performance in reproducing quantum dynamics results. The model correctly predicts the anisotropic response for rods and disks. However, while the error generally remains close to 10%, it is approximately double this for both the 2 × 2 × 4 Table 1. Relative permittivities for rod-and disk-like silver nanoparticles embedded in MgO matrices, calculated with quantum dynamics (QD) and the FEM. a Symbols ∥ and ⊥ indicate whether the inclusion has its long axis oriented parallel to the applied field and perpendicular to it, respectively. The field is oriented along z in each case, thus the ∥ notation indicates that the rod or disk long axis was along z, whereas the ⊥ notation indicates that the long axis was in the xy plane. τ is the permittivity error defined in equation (8). b Values from [24].
rod and the 4 × 4 × 2 disk. Analysis of the physical basis for discrepancies between the continuum and quantum models can both suggest improvements to the continuum model and provide insights difficult to obtain from quantum dynamics alone. These insights may inform design principles for high-permittivity composites. The continuum model performs best for the least anisotropic inclusions (such as the 2 × 2 × 3 rod and the 3 × 3 × 2 disk). The method can overestimate the permittivity of rodbased composites and underestimate that of disks. Several factors account for these discrepancies: (i) the anisotropy of inclusion polarization is only partially accounted for in a parameterization built on spherical inclusions; (ii) The ellipsoidal shape is an imperfect representation of rod-and disk-like inclusions, which do not exhibit sharp tips in the polarization direction, but may be better represented by cylindrical or rounded prism shapes; (iii) There is variability in the amount of padding between image inclusions in the different composites, restricted by volume loadings, that may not be fully captured by the interfacial region parameterized for the continuum model.
Notably, ab initio molecular dynamics simulations revealed differences in the ion mobilities of interfacial Mg and O ions between rod-and disk-like inclusions. Figure S7 shows that interfacial ion mobility is significantly increased around an Ag 18 disk relative to around an Ag 12 rod, which shows an interfacial ion MSD more similar to that of the Ag 8 nanocomposite from figure 6. This non-linear effect can help to account for the large ionic polarization in the studied disk-containing nanocomposites. In turn, it may also account for some of the shortcomings of a continuum model with a single parameter for interfacial permittivity, notably its overestimation of rod-like polarization and underestimation of disk-based nanocomposite polarization.
Two first-order forays into alternate parameterization options are discussed in the SI. Section S8 presents a parameterization that uses only low volume loading data to eliminate issues with interface and padding regions. The model significantly improves the fit for rodlike inclusions. The rounded-prism inclusion model discussed in section S9 underperforms the original ellipsoidal model.

Larger anisotropic inclusions.
We also calculated the permittivities of composites with unit cells that were too large to simulate quantum mechanically. We focused on anisotropic inclusions, with rods containing up to 56 silver atoms and aspect ratios of 14:2 and disks containing up to 392 atoms (also with aspect ratio 14:2). Volume loadings were varied by starting with a cell with two layers of MgO around inclusions and then sequentially adding layers of matrix. We approximated the dimensions of the inclusions using multiples of the dimensions of the 2 × 2 × 3 rod and 3 × 3 × 2 disk studied above (details in section S10) and used an ellipsoidal model with the optimal parameters from section 3.2.3. The dependence of relative permittivity on inclusion size and volume loading for these inclusions is presented in figure 8.
The volume ratio and anisotropy trends observed in quantum dynamics are conserved for the large inclusions in the FEM: the relative permittivity increases monotonically with the volume ratio. Small increases with axis ratios are seen for rod-like inclusions, when the inclusion is oriented parallel to the field, and small decreases are seen in perpendicular orientations (see figure 8(a)). The aspect ratio dependence of the permittivity is stronger for disks, as seen in figure 8(b), where each axis ratio represents a significantly larger number of atoms (and polarizable electrons) in the inclusion. It is worth noting the incipient saturation behavior for larger rods in parallel orientation to the field, at higher volume loadings. The requirement of a dielectric interface presents an upper limit to the volume loading that can be treated with the present FEM, and the higher loadings presented here for large rods approach that limit.

Conclusions
We propose a finite element model for the dielectric response of nanocomposites of silver nanoparticles, described as polarizable dielectric spheres and ellipsoids, embedded in a dielectric matrix. A methodology for fitting the FEM to a quantum mechanical polarization model is also presented. When fit to a suitable quantum mechanical training set, an ellipsoidal model faithfully reproduces the polarizability of gas phase particles, including the anisotropy due to shape effects. For nanocomposite systems, the FEM requires the use of an accurate definition for the electron delocalization volume, as well as the addition of an interfacial layer to the model. A model containing ellipsoidal inclusions and an interfacial layer of high permittivity was found to reproduce the permittivities of nanocomposites with rod-and disk-like inclusions with reasonable fidelity, even though the interfacial permittivity was fitted on a set of nanocomposites with spherical inclusions. The FEM allows the estimation of permittivities for systems that cannot be treated quantum mechanically, such as systems with large anisotropic inclusions, over a large range of inclusion loadings. Furthermore, analysis of the FEM data provides insight into local field effects and the coupling of polarization between components, that is complementary to that obtained from quantum simulations.
To improve the FEM parameter set, further investigations are needed into other, nonellipsoidal representations for anisotropic particles embedded in matrix. A more general mapping of the atomistic interface onto the interfacial layer of the FEM is also being examined for model improvement. Additionally, a broader set of systems, including other matrix and inclusion materials, are currently being considered, in order to generalize the FEM implementation for an array of currently relevant dielectric ceramics with highly-polarizable inclusions.

Data availability statement
The data that support the findings of this study are openly available at the following URL: https://github.com/Paci-Group/NanoComposite_FEM.