We demonstrate a method to compute the dielectric spectra of fluids in molecular dynamics (MD) by directly applying electric fields to the simulation. We obtain spectra from MD simulations with low magnitude electric fields (≈0.01 V/Å) in agreement with spectra from the fluctuation–dissipation method for water and acetonitrile. We examine this method’s trade-off between noise at low field magnitudes and the nonlinearity of the response at higher field magnitudes. We then apply the Booth equation to describe the nonlinear response of both fluids at low frequency (0.1 GHz) and high field magnitude (up to 0.5 V/Å). We develop a model of the frequency-dependent nonlinear response by combining the Booth description of the static nonlinear dielectric response of fluids with the frequency-dependent linear dielectric response of the Debye model. We find good agreement between our model and the MD simulations of the nonlinear dielectric response for both acetonitrile and water.

## I. INTRODUCTION

Dielectric spectroscopy is a sensitive probe of the behavior of collections of molecules and individual species in solution. Experimental broadband dielectric spectroscopy can measure solution spectra across a broad frequency range from 10^{3} to 10^{16} Hz.^{1} These spectra have peaks caused by molecular-level processes ranging from ion pairing and simple rotations of species in solution to disruptions of hydrogen bond networks and distortion of ionic clouds.^{2–5} Molecular dynamics (MD) simulations provide a unique opportunity to resolve the molecular origin of these peaks, but the usual method for simulating dielectric spectra, the fluctuation–dissipation theorem,^{6–10} has numerical limitations at low frequencies.

The fluctuation–dissipation theorem relates the autocorrelation of polarization fluctuations to the frequency dependent dielectric susceptibility *χ*(*ω*),^{11,12}

where V is the volume, k_{b} is Boltzmann’s constant, T is the temperature, *ɛ*_{0} is the dielectric permittivity of free space, *ω* is the angular frequency (i.e., *ω* = 2*πf*), $P\u20d7$(0) is the initial polarization of the system, and $P\u20d7$(*t*) is the polarization at time *t*. The fluctuation–dissipation method allows for the calculation of the entire spectrum using a single simulation. However, this method requires simulation times on the order of hundreds of nanoseconds to obtain sufficient statistics^{13,14} and suffers from accuracy problems with decreasing frequency, as we demonstrate in this paper.

Here, we first describe our direct applied field approach for simulating the dielectric response and develop a model for the nonlinear, frequency-dependent dielectric response of a liquid of dipoles. We then compare the direct electric field method to the fluctuation–dissipation theorem for both water and acetonitrile, describing how the dielectric response varies with electric field magnitude. We demonstrate that at high electric field strengths, the response becomes nonlinear, saturating at low frequencies. Finally, we compare the results of our model with those from our molecular dynamics simulations as a function of frequency for both water and acetonitrile.

## II. METHODS

### A. Computational details

We model the real and imaginary dielectric response of water and acetonitrile with molecular dynamics using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS).^{15} Each solvent is simulated with 1500 molecules (4500 atoms and 9000 atoms for the water and acetonitrile simulations, respectively), using the extended single point charge (SPC/E) model for water^{16} and an all-atom force field for the acetonitrile.^{17} Electrostatic interactions are computed by the Particle-Mesh–Ewald method, collecting the electric dipole vector and total electric dipole. The Lennard–Jones interactions cut-off was set to 1.2 nm for all simulations. The system sizes have been kept identical for both direct electric field and fluctuation–dissipation theorem methods.

All equilibration simulations were run at 298.15 K and 1 atm for 2 ns under the isothermal–isobaric ensemble (NPT). This was sufficient time to minimize the volumetric fluctuations of the materials studied. The resulting density of water used in further simulations is 0.992 g/cm^{3}, and the acetonitrile reached a density of 0.779 g/cm^{3}. These results are in excellent agreement with previously reported computational values.^{6,17} We find the bulk dielectric constant for water to be 71.8, in good agreement with previously reported values for SPC/E water (69.9 from Ref. 6), which is less than the average experimental value of 78.3.^{18–21} For acetonitrile, we find the bulk dielectric constant of 22.0 (experimental value of 35.95^{22}), and a Debye relaxation time of 4.14 ps, in good agreement with the experimentally reported value of 3.51 ps.^{23} The equilibrated systems were then used for production runs using the canonical ensemble (NVT).

In the direct electric field method that we developed here, a sinusoidally oscillating external electric field was applied after the 2 ns of equilibration. Each simulation was run with an electric field of a different frequency, in the range of 0.1–1000 GHz. All production runs were done using a 2 fs time step, and systems ran for 50 ns with data collection performed at each step. The fluctuation–dissipation theorem simulations were also done with a 2 fs timestep, and the data for dipoles were collected every 4 fs. The dipole components of each time series are Fourier transformed via Fast Fourier Transformation (FFT) to obtain the dielectric spectra.

Note: Certain software is identified in this paper to foster understanding. Such identification does not imply recommendation or endorsement by the National Institute of Standards and Technology nor does it imply that the software identified is necessarily the best available for the purpose.

### B. Direct electric field method

We apply a time-varying electric field $E\u20d7(t)=E\u2061sin(\omega t)z\u0302$ along the *z*-direction of the simulation box and measure the corresponding induced polarization density *P*_{z}(*t*) (total simulation cell dipole moment per unit volume). For each electric field amplitude *E* and frequency *ω*, we fit *P*_{z}(*t*) to a series of harmonics,

where the coefficients *P*_{n} depend on *E* and *ω*, and here, we analyze up to the *n* = 5 harmonic order. From the first harmonic *n* = 1, we extract the complex susceptibility,

where *ɛ*_{0} is the permittivity of free space. This frequency-dependent response is nonlinear with *E* in general and reduces to the linear susceptibility accessed by fluctuation–dissipation calculations as *E* → 0.

### C. Nonlinear response model

To explore the physics underlying the simulated polarization in different solvents, we develop a model for the nonlinear frequency-dependent response of an idealized dipolar fluid. The static limit of the nonlinear response is well-developed for solvation models,^{24}, especially for those applied to electrochemical systems that introduce significant electric fields at the interface.^{25,26}

The nonlinear response can be derived from a competition between the internal free energy function $A(P\u20d7)$, in terms of the polarization density $P\u20d7$, and the interaction energy with the external field *E* in classical density functional theory within a local polarization density approximation.^{24,27,28} Going forward, we drop the vectors for a response along a single direction for simplicity. Specifically, from the partition function of a liquid with number density *N* of dipoles *p*_{0} each in an external field *E*, the free energy function *A*(*P*) can be shown to be implicitly given by

and *P* = *Np*_{0}*L*(*βp*_{0}*cE*), where *β* = 1/(*k*_{B}*T*), *L*(*x*) ≡ coth(*x*) − 1/*x* is the Langevin function, and *c* is a factor that corrects the external field to the local field experienced by the dipole (accounting both for Kirkwood dipole correlations^{29,30} and Clausius–Mossotti cavity fields^{31}).

We extend this nonlinear response model from the static limit to polarization dynamics,

where the dipole scattering rate with damping coefficient *γ* on the left side is driven by a net force combining the external field and the restoring force due to the internal free energy *A*(*P*). We neglect inertial effects with a term proportional to *∂*^{2}*P*/*∂t*^{2} above,^{32} and it is shown in Fig. S3 of the supplementary material that this only impacts the response at frequencies much higher than the Debye peak that we focus on below.

From Eq. (4), we can solve for the restoring force in terms of the polarization as

Next, we define the relaxation time in the low-field limit, *τ* ≡ *γχ*_{0}, where $\chi 0\u2261limE\u21920(\u2202P/\u2202E)=\beta Np02c/3$ is the linear susceptibility [using *L*(*x*) ≈ *x*/3 for *x* ≪ 1]. Using these relations, we can rearrange Eq. (5) as

a nonlinear first-order differential equation in *P*(*t*) with external forcing *E*(*t*).

We solve Eq. (7) numerically to extract the nonlinear frequency-dependent polarization response of the dipole model fluid at general field strengths and frequencies to compare with molecular dynamics simulations. Before analyzing these numerical results, note that Eq. (7) reproduces two well-known analytical limits by construction. First, in the static (low-frequency) limit, *∂P*(*t*)/*∂t* is negligible, and Eq. (7) reduces to *P* = *Np*_{0}*L*(*βp*_{0}*cE*), the standard Langevin saturation response or Booth equation.^{29} Next, in the low-field limit, *L*^{−1}(*x*) ≈ 3*x*, reducing Eq. (7) to *τ∂P*/*∂t* + *P* = *χ*_{0}*E*. This can be solved readily in the frequency domain to yield the Debye model susceptibility, *χ*(*ω*) ≡ *P*/*E* = *χ*_{0}/(1 − *iωτ*).

## III. RESULTS AND DISCUSSION

### A. Linear response

Figure 1(a) compares the calculated frequency-dependent susceptibility from the fluctuation–dissipation theorem and the direct electric field approaches. Across the entire spectrum, the two methods are in very good agreement, indicating that the response for water is sufficiently linear at an applied field strength of 0.01 V/Å. Simulations of water perturbed by a 0.01 V/Å static electric field have been shown by others to have an approximately linear response^{7,33} as well. Fitting the data to the Debye model, the Debye relaxation times (direct electric field method: 10.3 ps; fluctuation–dissipation theory method: 9.7 ps, standard deviation of 0.1 ps) are in good agreement with simulation results from sufficiently long runs (10.72 ps, Cole–Cole fit).^{6} We note that the convergence challenges of fluctuation–dissipation theorem have resulted in a wide range of water relaxation times in the literature.^{34}

Despite the level of agreement in the resulting spectra, the variance differs between the approaches. Figure 1(a) depicts the average signal for both methods, with variance calculated from five separate simulations of 100 ns (FDT) and five 10 ns simulation intervals (from a 50 ns simulation) at each frequency (direct electric field).

The fluctuation–dissipation variance for the imaginary susceptibility peaks at around 10 GHz. These data are collected fromautocorrelations and are constrained to be zero at ω= 0 Eq. (1). Depending on the length of the run, very little to no autocorrelation data are collected at such low frequencies. This causes the variance displayed in Fig. 1(a) to be low, while not representing all sources of the uncertainty.

Above ∼10 GHz, the variances in the direct electric field method decrease with frequency, although this specific value may depend on the solvent. Therefore, to produce a similar signal to noise for each frequency, one could scale the individual simulation times by $1\omega $.

Variable length runs, coupled with a modified collection frequency of the polarization density, have been used with the fluctuation–dissipation theorem to produce data more evenly across a broad frequency range.^{35}

To compare the convergence between the direct electric field method and the fluctuation–dissipation theorem, Fig. 1(c) plots the average and variance of the low frequency (f = 1 GHz) dielectric constant from the two methods as a function of simulation time. The variances in the fluctuation–dissipation method simulations are larger than those from the direct electric field method, and they converge slower with simulation time.

One way to further reduce the variance in the spectrum collected with the direct electric field method is to increase the magnitude of the applied electric field so that the signal to noise is larger. However, if the field is too large, the water response will be nonlinear.

To explore the effects of the electric field magnitude on the polarization response, we have simulated the dielectric spectrum of water with applied field magnitudes ranging from 0.0001 to 0.01 V/Å. The resulting real and imaginary dielectric spectra are shown in Figs. 2(a) and 2(b), respectively.

With decreasing electric field magnitude, the real and imaginary components of the susceptibility become noisier. For electric field magnitudes of greater than 0.0001 V/Å, the noise appears to have the same frequency dependence noted for 0.01 V/Å. For 0.0001 V/Å, the signal is so close to the noise floor that this dependence is no longer clear, and the response of the fluid to the field is almost completely obscured by natural fluctuations in the fluid. When no field is applied, the polarization exhibits the inherent statistical noise of the water box. This noise fluctuation can be used to produce the dielectric spectrum via the fluctuation–dissipation theorem [Eq. (1)].

### B. Low frequency nonlinear response

We evaluate the low frequency nonlinear dielectric response of water and acetonitrile to develop guidance for choosing a perturbing electric field magnitude.

We first evaluate the nonlinear dielectric response of water at a frequency well below that of the Debye peak, 0.1 GHz. Figure 3(a) shows the polarization of a box of water molecules simulated with molecular dynamics in response to an applied sinusoidal electric field at a field strength of 0.1 V/Å. This applied field strength results in a nonlinear polarization of the water, and we include harmonics beyond the first order to improve the fitting of the response.

Figure 3(b) plots the response up to the fifth order as a function of electric field magnitude. At the lowest field strengths, noise obscures the contributions of *P*_{3} and *P*_{5}. When the field magnitude is sufficiently low, the polarization fitting captures the response from the natural fluctuations in the fluid rather than a response from an external field. These contributions will create an effective floor to the fit values for the direct electric field method, as seen for the lowest field strength fifth order harmonics in Fig. 3(b).

At slightly larger field strengths, the contributions of *P*_{3} and *P*_{5} appear to be linear on the log–log axes. This is expected when a Taylor series of the form *P*_{n} = *c*_{n}*E*^{n} + *c*_{n+2}*E*^{n+2} + ⋯ is valid for the polarization response, leading to *P*_{n} ∝ *E*^{n} for small enough *E*. With increasing field strengths, terms beyond the leading order become important, and the higher-order responses also exhibit saturation. This general behavior has been observed experimentally across a broad range of materials and electronic circuits.^{37}

Figure 3(b) illustrates how the magnitudes of the nonlinear components increase with E. However, the statistical noise is expected to increase with decreasing applied electric field magnitude. To balance between the statistical noise and the nonlinearity, one could select a value for *E* where the magnitude of the error is roughly equal to the size of the next order term. If this error is too large, the statistical noise will need to be reduced, e.g., by increasing the number of molecules in the simulation.

To predict the nonlinear behavior of the dielectric response from the first order harmonic fits, we model the low-frequency nonlinear response of water and acetonitrile with the Booth equation,^{29,30,36} following the approach of^{31}

where *ɛ*_{n} is the high frequency limit of the dielectric constant, *ɛ*_{b} is the low frequency limit of the bulk dielectric constant, and the fit parameter *α*_{E} = *βp*_{0}*c* is related to the molecular dipole moment and dipole correlation factor *c* discussed in Sec. II C. The extracted *α*_{E} corresponds to *c* = 4.0 and 0.89 for water and acetonitrile, respectively, in agreement with previous simulations.^{17,29–31} The nonlinearity of the response of water at low frequency, 0.1 GHz, agrees well with data from perturbation of SPC/E water by a *static* electric field^{36} (Fig. 4). The parameters for the Booth equation fit for water and acetonitrile are in Table III A.

The Booth equation parameters can provide an estimate for the nonlinear error in the response as a function of applied field strength (Table I).

### C. Frequency-dependent nonlinear response

We next evaluate the nonlinear dielectric response as a function of frequency and identify its molecular origins. Figure 5 plots the polarization of a model (Sec. II C) solvent and MD-simulated water and acetonitrile in response to a sinusoidal applied electric field. The polarization is plotted for frequencies below, at, and above the Debye peak at several electric field strengths. At a frequency below the Debye peak, the linear response is approximately sinusoidal, and at larger field strengths, the response saturates for all solvents. At this low frequency, the model solvent describes this saturation behavior well for both water and acetonitrile because the primary source of this saturation is the rotational saturation of the molecules as described by the Booth equation.

At the frequency of the Debye peak, the polarization response first sharply increases with time as the electric field strength increases until the response saturates, and then the response decreases smoothly as the applied electric field strength decreases. This behavior is seen for the model and the MD simulated solvents.

At a frequency above the Debye peak, the linear and nonlinear responses are out of phase from the applied field. The response does not appear to saturate even at large field strengths. However, the behavior of the model is not fully consistent with that of the molecular dynamics solvents. In both MD solvents, the polarization response is not fully phase shifted to match that of the model. This disagreement may be caused by both the acetonitrile and the water dielectric spectra containing features that are not fully Debye-like.^{4} A more complex model for the linear response of the fluid could improve the agreement between the model and the MD simulation results.

An additional disagreement between the model and MD simulation results can be found for the high field strength acetonitrile simulation. The zero crossing of the water polarization is roughly consistent across a large range of applied field strengths, whereas the acetonitrile polarization shifts with field strength. The acetonitrile force field includes flexibility at the molecular level, whereas the water force field is rigid. When the acetonitrile is simulated as a rigid molecule, we find that this shift reduces significantly, as seen in Fig. 6.

Finally, we evaluate the harmonic components of the nonlinear response of the model and compare it to that of the water and acetonitrile simulations (Fig. 7). We note the excellent agreement between the linear and nonlinear response of the model and both the water and acetonitrile.

## IV. CONCLUSIONS

In this work, we demonstrate a direct method for calculating the dielectric spectrum of solvents using an alternating current electric field. By using the direct electric field approach, we can focus on specific frequencies, more easily collect data at low frequencies, and investigate the nonlinearity of the dielectric response. We find that a field strength of 0.01 V/Å provides good agreement with the fluctuation–dissipation results for water and acetonitrile, with a primarily linear response at this field strength. At higher field strengths, the response becomes more nonlinear.

We develop a numerical model combining the frequency-dependence of the Debye response with the nonlinearity of the Booth equation. We apply this model and demonstrate that the nonlinear components of the water and acetonitrile responses are in good agreement with our predictions. Specifically, we find that the response becomes more linear at frequencies above the Debye peak, although the model predicts a fully Debye linear response, and both water and acetonitrile have deviations from that linear response. Our numerical model could be further developed to include other features of the linear response, such as resonances (e.g., bond vibrations), electronic polarization response, and more complex relaxations (e.g., Cole–Cole and Havriliak–Negami,^{38} using the framework of memory functions^{39}). This will also require significant advances in describing the nonlinear response associated with those spectral features and is expected to be highly molecule-specific.

This work opens new research possibilities for both simulating dielectric spectra with the direct field method and modeling the frequency dependence of the nonlinear response of polar fluids.

## SUPPLEMENTARY MATERIAL

See the supplementary material for Fig. S1 that plots the time-averaged raw polarization data without normalization. Figure S2 provides the raw coefficients of the polarization of water and acetonitrile from fitting the simulation data, without normalization. Figure S3 plots the normalized polarization coefficients of the model solvent with the inclusion of an inertial term, $M\u22022P(t)\u2202t2$ to the left side of Eq. (5) of the main text.

## ACKNOWLEDGMENTS

Michael Woodcox acknowledges support from the U.S. Department of Commerce, National Institute of Standards and Technology under the financial assistance Award No. 70NANB22H002. Ravishankar Sundararaman acknowledges support from the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award No. DE-SC0022247. We thank Alta Fang for her fluctuation–dissipation script and Abdenacer Idrissi for providing molecular dynamics force field parameters for acetonitrile used in this work.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Michael Woodcox**: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Avik Mahata**: Formal analysis (supporting); Investigation (supporting); Validation (equal). **Aaron Hagerstrom**: Conceptualization (equal); Formal analysis (supporting); Methodology (equal); Validation (equal). **Angela Stelson**: Conceptualization (equal); Methodology (equal); Validation (equal); Visualization (supporting). **Chris Muzny**: Conceptualization (equal); Project administration (equal); Supervision (equal); Validation (equal). **Ravishankar Sundararaman**: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Writing – review & editing (equal). **Kathleen Schwarz**: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are openly available on GitHub at https://github.com/usnistgov/DielectricSpectroscopyTheory.