Supplemental Information: Transport mechanisms underlying ionic conductivity in nanoparticle-based single-ion electrolytes

Transport mechanisms underlying ionic conductivity in nanoparticle-based single-ion electrolytes Sanket Kadulkar,† Delia J. Milliron,† Thomas M. Truskett,∗,‡ and Venkat Ganesan∗,† †McKetta Department of Chemical Engineering, University of Texas at Austin, Austin, Texas 78712, USA ‡McKetta Department of Chemical Engineering and Department of Physics, University of Texas at Austin, Austin, Texas 78712, USA

The simulation cell is a cuboid with length L z in the direction perpendicular to the surfaces and length L = 13.46 σ in the two (periodically replicated) directions parallel to the flat surfaces. Based on the experimental study by Schaefer et al., 2 we fix the grafting density of anions to be equal to that of the polymer chains (0.0993 per σ 2 ); i.e., 18 polymer chains and 18 anions per surface. This polymer grafting density corresponds to ≈ 0.42 times the experimentally reported grafting density (190 poly(ethylene glycol) chains on 7 nm silica nanoparticle). 2 To maintain charge neutrality, the number of cations in the cell is equal to the number of anions (i.e., 36 ion pairs total). The number of oligomeric solvent chains in the system is varied depending on L z and the anion size, while the number of every other component in the system was fixed. When solvent chains are present, their number density is chosen to maintain a system volume fraction of 0.444 (Table S1). Such a value is equivalent to that of a fluid of monomers with number density of 0.85 σ −3 . Below a critical value of L z , no solvent is present, and the system volume fraction monotonically increases with increasing confinement due to the presence of tethered polymer chains and anions.  pair potential. Every graft bead is assumed to be charge neutral with no embedded dipoles.
The forces on the graft beads are fixed to be zero throughout the simulations to replicate nanoparticle-tethered polymer chains and anions. The graft beads are each bonded to a freely moving charge neutral bead with diameter σ with no dipole moment. 2 The charge neutral beads are then bonded to either the first bead of a polymer chain or an anion such that a uniform distribution of both polymer chains and anions is maintained across the flat surface. The polymer chains are randomly grown from the first bead attached to the charge neutral bead. Subsequently, the oligomeric solvent chains were grown randomly and dispersed throughout the simulation cell.
The simulations are executed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package. 3 All systems are initially energy minimized using the steepest descent algorithm followed by the conjugate gradient algorithm to remove heavy overlaps.
Short runs are then performed in the microcanonical (NVE-limit) ensemble during which the beads are restricted from moving more than 0.1 σ for 100 τ , where τ = σ 2 m/ , m is the monomer bead mass (assumed equal for all beads), and is the strength of the monomermonomer Lennard-Jones (LJ) attraction for the polymer described in Sec. S1.2 ( ≈ 5.603 × 10 −21 J). 1 This effectively removes softer particle overlaps and further relaxes the system. S-3

S1.2 Force Field Details
Most of the interaction potentials and the corresponding force field parameters are adapted from a previous study 6 in our group on Stockmayer polymer electrolytes. All beads interact through a repulsive Lennard-Jones potential E LJ : where r ij is the distance between the centers of two beads i and j. σ ij is the arithmetic mean of the bead diameters, σ i and σ j , determined as σ ij = (σ i + σ j )/2.
Further, the interaction of the beads with the nanoparticle is accounted through a repulsive LJ potential E wall : where r is the distance between the center of the bead i and the nearest wall (nanoparticle surface), and σ i is the diameter of the bead i.
All components bonded to each other (monomer-monomer, graft bead-charge neutral bead, charge neutral bead-anion, charge neutral bead-monomer) interact through the finitely extensible nonlinear elastic (FENE) potential: We use the spring constant k = 30 /σ 2 and cutoff radius R 0 = 1.5 × (σ i + σ j )/2. 7 In addition to the LJ potential, ions interact through the Coulomb potential E zz : where e is the fundamental proton charge, z is either +1 (cations) or −1 (anions), and ε 0 is Further, dipoles embedded in monomer beads interact through the following potential: where x ij is zero if the the monomers are directly bonded and unity otherwise.
In addition to the translational degrees of freedom, the presence of freely rotating dipoles in the monomers necessitates accounting for the associated rotational degrees of freedom and torques on the monomers. The monomers were simulated as finite-size particles with inertia of rotation I = 0.025 mσ 2 . 8 Sources of torque T on the monomer beads include ion-dipole interactions: and dipole-dipole interactions: To that end, we assume the contributions to the effective pair potential from the bare nanoparticles and the functionalizing components, respectively, to be additive, where U bare (r) is the pair interaction between bare silica nanoparticles and U func (r) is the effective pair interaction between non-interacting nanoparticles due to the solvent and the components functionalized on the nanoparticle surfaces.
A semiempirical pair potential proposed by Lee and Hua 9 is considered to represent the pair interactions of coarse-grained (bare) silica nanoparticles, where bare and α are nanoparticle size-dependent parameters. 9,10 Considering the nanoparticle size reported in the experimental study and the parametrization in our study, these S-6 parameters are calculated to be bare = 40.84 and α = 33.9.
We calculate the potential of mean force (PMF) between two non-interacting particles as a function of their distance to quantify the contribution of functionalized components (including solvent) to the effective pair potential between particles. Explicitly, we approximate the force between the spherical particles, F (r), based on the Derjaguin approximation for two spheres of equal radii R, 11 where W (D) planes is the interaction energy per unit area of the two functionalized planar surfaces at a separation distance D. The pair potential due to the functionalized components, U func (r), is expressed in terms of the mean force as

S2.2 Simulation Methodology
We perform molecular dynamics (MD) simulations to obtain the equilibrium configurations of nanoparticles using the LAMMPS 3 package. To carry out such simulations, we start from a random initial configuration of particles and equilibrate the system for 5 × 10 8 steps in the canonical ensemble with time step dt = 0.001 σ 2 m/ . The temperature considered is T = 373 K. 2 Temperature is maintained using a Nosé-Hoover thermostat. 4,5

S3 Simulations for the ionic conductivity
On-lattice kinetic Monte Carlo (kMC) simulations are based on a master equation 12 in which the probability distribution P (ω i , t) for being in state ω i at time t evolves as follows: where k ij is the transition (hopping) rate from state ω i to ω j . In the context of our system, the lattice sites with different transition rates or different affinities to host a cation are labeled with unique states. A transition or a move corresponds to the hop of a cation to a neighboring site (in the ±x, ±y, ±z direction). In this work, we implemented the direct kMC method, 13-15 one of the widely used rejection-free algorithms to solve the master equation (Eq. 13). Based on the analysis of lattice spacing effects in kMC simulations (Section S5.2), we choose a lattice spacing of 0.32 σ. We first map the nanoparticle centers onto the lattice sites. With these sites as center, we then assume the sites within a sphere of radius 0.5 σ NP to be occupied by the nanoparticles. The simulation cell consists of 1024 noninteracting tracer S-8 cations placed randomly on sites not occupied by the particle. We note that a large number of tracer cations are present in order to increase sampling in the estimation of the cation selfdiffusion coefficient (i.e., their concentration in the simulation is not reflective of a physical concentration of cations in the system). For our system, the nanoparticles are the only source of cations. Hence, N cat is proportional to the particle volume fraction. Based on the choice of coarse-grained parameters in our study, we assume N cat = 80 × N , to appropriately map the number of Li + for every particle as reported in the experimental system. 2 Every kMC simulation run lasted for 3 × 10 6 moves. From the resulting trajectories of the tracers, we extract the mobility of the cations using the Einstein-Smoluchowski relation: where (r(t) − r(0)) 2 is the mean squared displacement (MSD) of cations at time t. The cation diffusivity: where (r (t) − r (0)) 2 surface is the mean squared displacement (MSD) of cations in the directions parallel to the surface in the region, 0 ≤ z ≤ L surface z , at time t. The survival probability of cations P (t) denotes the averaged probability for the cations to remain in the surface region at time t. Figure S2 shows the cation dynamics in the non-overlapping surface region for all the cases considered in this study.

S4.2 Calculating D in the overlapping surface region
For the cases with L z < 2 L surface z , the surface regions associated with both the nanoparticles merge. In this scenario, the cation transport through the entire region is in the vicinity of anions. Hence, the local cation diffusivity is assumed to be uniform throughout the region between the two surfaces, and is calculated using the Einstein-Smoluchoski relation: As shown in Figure S3, the cations reach the diffusive regime for all the cases.
S-10 We denote the region occupied by the solvent molecules as "bulk", and we assume that the cation diffusivity in this region is independent of the L z value and the anion size. To that end, we calculate the D value in the bulk region for a system with large separation distance between the two flat surfaces (L z = 30 σ). For the low polarity case, due to the negligible cation presence in the bulk region, the methodology employed to calculate the local cation diffusivity in Section S4.1 cannot be extended to extract the local cation diffusivity in the bulk region. For the high polarity case, the local diffusivities in the bulk region, calculated using the similar methodology, were inconclusive due to the faster movement of cations in the bulk region resulting in considerable statistical noise. Alternatively, to facilitate the cation presence in the bulk region, the cations were initialized within the region, 7 σ ≤ z ≤ 23 σ, where z in this context denotes the distance from either surface. The cations were restricted to move outside this region throughout the dynamics. Additionally, to maintain a uniform distribution of cations throughout this region, an external biasing potential is applied on the cations. Explicitly, in addition to its interaction with the solvent chains, the cations experience an additional harmonic spring potential E bias : where bias = 2 and 0.25 for the low and high polarity cases, respectively. The above parameters were chosen to ensure uniform distribution of cations throughout the bulk region, as shown in Figure S4.
Restricting the movement of cations in the bulk region allowed for calculating reliable D values (using Eq. 16). The cation MSDs in the bulk region are shown in Figure S5.

S4.4 Assigning transition rates in the on-lattice kMC model
Lattice sites with similar local diffusivities and affinities to host a cation are denoted as unique states ω i (as defined in Eq. 13) for cations (tracers). To that end, the sites occupied by the nanoparticles are denoted as state ω 0 . Based on the distribution of cations and local cation diffusivities in the MD simulations, we define six unique states (ω i , i = 1,....,6) for the tracers to be present on the sites not occupied by nanoparticles ( Figure S6).
Explicitly, sites belonging to states ω 1 , ω 2 and ω 3 represent the "surface" sites. Sites S-13 with unique local diffusivities are assigned to states ω 2 and ω 3 as shown in Figure S6(c). The "bulk" sites were represented by sites corresponding to the remaining states (ω 4 , ω 5 and ω 6 ).
In contrast to the states representing the surface sites, the different states representing the bulk region offer similar local diffusivities but exhibit different affinities to host a cation.
Precisely, to capture the non-uniform cation distribution in the bulk region, observed for S-14 the cases with L z = 15 σ ( Figure S7), we define two sub-regions in the bulk region denoted by states ω 4 and ω 6 . Further, for the case with relatively smaller separations, the cation presence in the bulk region between the nanoparticles is slightly enhanced ( Figure S7) due to the relatively stronger electrostatic interactions with the anions tethered to either surface.
Hence, we define a unique state ω 5 as illustrated in Figures S6(b) and S6(c) representing the sites corresponding to the bulk region between two nanoparticles with relatively smaller separation distance.
The value of k ij is then fixed based on the condition in Eq. 18.

S5.1 Effective pair potential and nanoparticle morphologies
The results for the effective pair potential, U (r), between the functionalized nanoparticles are presented in Figure S8(a). We observe that the effective interaction between nanoparticles is repulsive for all separation distances due to the dominant contribution from the functionalized components (U func (r)) to the effective pair potential. From the resulting equilibrium morphologies at different nanoparticle volume fractions, we quantify the proximity between neighboring nanoparticles. To that end, we denote d to be the shortest distance between the surface of a particle and the surface of its nearest neighbor, averaged over all particles. From the results reported in Figure S8(b), it is evident that at low volume fractions, the nanoparticles are distant from each other due to the repulsive interactions between the nanoparticles ( Figure S8(a)). The d values are seen to decrease with increasing nanoparticle loading. Figure S8: (a) Effective pair interactions between functionalized nanoparticles along with the contribution of individual components. (b) Shortest distance between the surfaces of a particle and its nearest neighbor, averaged over all particles, at different nanoparticle loadings.

S5.3 Mean-squared displacement of tracer cations from kMC simulations
The cation diffusion coefficient in the multiparticle system was extracted from the slope of the long-time mean-squared displacement of tracer cations (Eq. 14). From the MSDs shown in Figure S11, it is evident that the tracers reach the diffusive regime in the kMC simulations for all the cases.

S5.4 Anion-cation coordination
To rationalize the lower cation diffusivities in the overlapping surface regions compared to those in the non-overlapping surface region for the low polarity case (Figure 2c in the main text), we characterized the ion pair electrostatic interactions in the overlapping surface regions (L z ≤ L surface z ). To that end, we report the average anion-cation coordination number (n) for two cases with different anion sizes. The cation size was fixed to be σ cat = 0.5 σ.
The cutoff distance for anion-cation coordination was decided based on the anion-cation radial distribution function ( Figure S10(a)). From the results displayed in Figure S10 Figure S12(a), we present results for the fraction of cations hosted on different site types. Explicitly, the surface sites are categorized into non-overlapping and overlapping surface sites. It is evident that for both the cases with different anion sizes, the cation presence on the overlapping surface sites monotonically increases with particle volume fraction at higher particle loadings (φ > 0.25).
This results in a significant influence of the cation transport through the overlapping surface region on the overall cation diffusivities. To further support this hypothesis, we repeat our S-23 calculations for a fictitious case with D values in the overlapping surface region being equal to those in the non-overlapping surface region. To that end, we modify the kMC inputs in Tables S2 and S4. Explicitly, we set k 22 and k 33 values to be equal to k 11 . From the results displayed in Figure S12 Figures S13(a) and S13(b) shows the cation distribution on different site types to be similar for all cation sizes.

S5.7 Cation-solvent association autocorrelation function
The metrics employed to understand the trend for cation diffusivities in the bulk region with increasing cation size is the cation-solvent association autocorrelation function C(t). The

S-24
association autocorrelation function is defined as: where h(t) takes a value of one if a cation-solvent association is present at time t and zero otherwise. 17 The autocorrelation function C(t) shows on average how long a cation-solvent pair remains associated. The cutoff distance for cation-solvent association is decided based on the cation-solvent radial distribution functions depicted in Figure S14(a). From the results displayed in Fig-ure S14(b), it is evident that for the high polarity case (µ = 2 µ EO ), the cation-solvent association function decays faster with increasing cation size. This indicates that the increased cation diffusivity in the bulk region with increasing cation size can be understood by the movement of cations in contact with the dipolar solvent beads. For systems with highly asymmetric sizes (σ cat = 0.25 σ), the cation strongly associates with solvent chains and cation transport by vehicular motion also seems to contribute considerably. S-25